Actual source code: vecnest.c


  2: #include <../src/vec/vec/impls/nest/vecnestimpl.h>

  4: /* check all blocks are filled */
  5: static PetscErrorCode VecAssemblyBegin_Nest(Vec v)
  6: {
  7:   Vec_Nest *vs = (Vec_Nest *)v->data;
  8:   PetscInt  i;

 10:   PetscFunctionBegin;
 11:   for (i = 0; i < vs->nb; i++) {
 12:     PetscCheck(vs->v[i], PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Nest  vector cannot contain NULL blocks");
 13:     PetscCall(VecAssemblyBegin(vs->v[i]));
 14:   }
 15:   PetscFunctionReturn(PETSC_SUCCESS);
 16: }

 18: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
 19: {
 20:   Vec_Nest *vs = (Vec_Nest *)v->data;
 21:   PetscInt  i;

 23:   PetscFunctionBegin;
 24:   for (i = 0; i < vs->nb; i++) PetscCall(VecAssemblyEnd(vs->v[i]));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: static PetscErrorCode VecDestroy_Nest(Vec v)
 29: {
 30:   Vec_Nest *vs = (Vec_Nest *)v->data;
 31:   PetscInt  i;

 33:   PetscFunctionBegin;
 34:   if (vs->v) {
 35:     for (i = 0; i < vs->nb; i++) PetscCall(VecDestroy(&vs->v[i]));
 36:     PetscCall(PetscFree(vs->v));
 37:   }
 38:   for (i = 0; i < vs->nb; i++) PetscCall(ISDestroy(&vs->is[i]));
 39:   PetscCall(PetscFree(vs->is));
 40:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVec_C", NULL));
 41:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecs_C", NULL));
 42:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVec_C", NULL));
 43:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL));
 44:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL));

 46:   PetscCall(PetscFree(vs));
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
 51: {
 52:   Vec_Nest *bx = (Vec_Nest *)x->data;
 53:   Vec_Nest *by = (Vec_Nest *)y->data;
 54:   PetscInt  i;

 56:   PetscFunctionBegin;
 57:   VecNestCheckCompatible2(x, 1, y, 2);
 58:   for (i = 0; i < bx->nb; i++) PetscCall(VecCopy(bx->v[i], by->v[i]));
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
 63: {
 64:   Vec_Nest *bx = (Vec_Nest *)x->data;
 65:   Vec       Y;
 66:   Vec      *sub;
 67:   PetscInt  i;

 69:   PetscFunctionBegin;
 70:   PetscCall(PetscMalloc1(bx->nb, &sub));
 71:   for (i = 0; i < bx->nb; i++) PetscCall(VecDuplicate(bx->v[i], &sub[i]));
 72:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y));
 73:   for (i = 0; i < bx->nb; i++) PetscCall(VecDestroy(&sub[i]));
 74:   PetscCall(PetscFree(sub));
 75:   *y = Y;
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

 79: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
 80: {
 81:   Vec_Nest   *bx = (Vec_Nest *)x->data;
 82:   Vec_Nest   *by = (Vec_Nest *)y->data;
 83:   PetscInt    i, nr;
 84:   PetscScalar x_dot_y, _val;

 86:   PetscFunctionBegin;
 87:   VecNestCheckCompatible2(x, 1, y, 2);
 88:   nr   = bx->nb;
 89:   _val = 0.0;
 90:   for (i = 0; i < nr; i++) {
 91:     PetscCall(VecDot(bx->v[i], by->v[i], &x_dot_y));
 92:     _val = _val + x_dot_y;
 93:   }
 94:   *val = _val;
 95:   PetscFunctionReturn(PETSC_SUCCESS);
 96: }

 98: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
 99: {
100:   Vec_Nest   *bx = (Vec_Nest *)x->data;
101:   Vec_Nest   *by = (Vec_Nest *)y->data;
102:   PetscInt    i, nr;
103:   PetscScalar x_dot_y, _val;

105:   PetscFunctionBegin;
106:   VecNestCheckCompatible2(x, 1, y, 2);
107:   nr   = bx->nb;
108:   _val = 0.0;
109:   for (i = 0; i < nr; i++) {
110:     PetscCall(VecTDot(bx->v[i], by->v[i], &x_dot_y));
111:     _val = _val + x_dot_y;
112:   }
113:   *val = _val;
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
118: {
119:   Vec_Nest   *bx = (Vec_Nest *)x->data;
120:   Vec_Nest   *by = (Vec_Nest *)y->data;
121:   PetscInt    i, nr;
122:   PetscScalar x_dot_y, _dp, _nm;
123:   PetscReal   norm2_y;

125:   PetscFunctionBegin;
126:   VecNestCheckCompatible2(x, 1, y, 2);
127:   nr  = bx->nb;
128:   _dp = 0.0;
129:   _nm = 0.0;
130:   for (i = 0; i < nr; i++) {
131:     PetscCall(VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y));
132:     _dp += x_dot_y;
133:     _nm += norm2_y;
134:   }
135:   *dp = _dp;
136:   *nm = _nm;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
141: {
142:   Vec_Nest *bx = (Vec_Nest *)x->data;
143:   Vec_Nest *by = (Vec_Nest *)y->data;
144:   PetscInt  i, nr;

146:   PetscFunctionBegin;
147:   VecNestCheckCompatible2(y, 1, x, 3);
148:   nr = bx->nb;
149:   for (i = 0; i < nr; i++) PetscCall(VecAXPY(by->v[i], alpha, bx->v[i]));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
154: {
155:   Vec_Nest *bx = (Vec_Nest *)x->data;
156:   Vec_Nest *by = (Vec_Nest *)y->data;
157:   PetscInt  i, nr;

159:   PetscFunctionBegin;
160:   VecNestCheckCompatible2(y, 1, x, 3);
161:   nr = bx->nb;
162:   for (i = 0; i < nr; i++) PetscCall(VecAYPX(by->v[i], alpha, bx->v[i]));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
167: {
168:   Vec_Nest *bx = (Vec_Nest *)x->data;
169:   Vec_Nest *by = (Vec_Nest *)y->data;
170:   PetscInt  i, nr;

172:   PetscFunctionBegin;
173:   VecNestCheckCompatible2(y, 1, x, 4);
174:   nr = bx->nb;
175:   for (i = 0; i < nr; i++) PetscCall(VecAXPBY(by->v[i], alpha, beta, bx->v[i]));
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
180: {
181:   Vec_Nest *bx = (Vec_Nest *)x->data;
182:   Vec_Nest *by = (Vec_Nest *)y->data;
183:   Vec_Nest *bz = (Vec_Nest *)z->data;
184:   PetscInt  i, nr;

186:   PetscFunctionBegin;
187:   VecNestCheckCompatible3(z, 1, x, 5, y, 6);
188:   nr = bx->nb;
189:   for (i = 0; i < nr; i++) PetscCall(VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]));
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
194: {
195:   Vec_Nest *bx = (Vec_Nest *)x->data;
196:   PetscInt  i, nr;

198:   PetscFunctionBegin;
199:   nr = bx->nb;
200:   for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
205: {
206:   Vec_Nest *bx = (Vec_Nest *)x->data;
207:   Vec_Nest *by = (Vec_Nest *)y->data;
208:   Vec_Nest *bw = (Vec_Nest *)w->data;
209:   PetscInt  i, nr;

211:   PetscFunctionBegin;
212:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);
213:   nr = bx->nb;
214:   for (i = 0; i < nr; i++) PetscCall(VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]));
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
219: {
220:   Vec_Nest *bx = (Vec_Nest *)x->data;
221:   Vec_Nest *by = (Vec_Nest *)y->data;
222:   Vec_Nest *bw = (Vec_Nest *)w->data;
223:   PetscInt  i, nr;

225:   PetscFunctionBegin;
226:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);
227:   nr = bx->nb;
228:   for (i = 0; i < nr; i++) PetscCall(VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]));
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: static PetscErrorCode VecReciprocal_Nest(Vec x)
233: {
234:   Vec_Nest *bx = (Vec_Nest *)x->data;
235:   PetscInt  i, nr;

237:   PetscFunctionBegin;
238:   nr = bx->nb;
239:   for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
244: {
245:   Vec_Nest *bx = (Vec_Nest *)xin->data;
246:   PetscInt  i, nr;
247:   PetscReal z_i;
248:   PetscReal _z;

250:   PetscFunctionBegin;
251:   nr = bx->nb;
252:   _z = 0.0;

254:   if (type == NORM_2) {
255:     PetscScalar dot;
256:     PetscCall(VecDot(xin, xin, &dot));
257:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
258:   } else if (type == NORM_1) {
259:     for (i = 0; i < nr; i++) {
260:       PetscCall(VecNorm(bx->v[i], type, &z_i));
261:       _z = _z + z_i;
262:     }
263:   } else if (type == NORM_INFINITY) {
264:     for (i = 0; i < nr; i++) {
265:       PetscCall(VecNorm(bx->v[i], type, &z_i));
266:       if (z_i > _z) _z = z_i;
267:     }
268:   }

270:   *z = _z;
271:   PetscFunctionReturn(PETSC_SUCCESS);
272: }

274: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
275: {
276:   PetscFunctionBegin;
277:   /* TODO: implement proper MAXPY. For now, do axpy on each vector */
278:   for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
283: {
284:   PetscFunctionBegin;
285:   /* TODO: implement proper MDOT. For now, do dot on each vector */
286:   for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
287:   PetscFunctionReturn(PETSC_SUCCESS);
288: }

290: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
291: {
292:   PetscFunctionBegin;
293:   /* TODO: implement proper MTDOT. For now, do tdot on each vector */
294:   for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
299: {
300:   Vec_Nest *bx = (Vec_Nest *)x->data;
301:   PetscInt  i, nr;

303:   PetscFunctionBegin;
304:   nr = bx->nb;
305:   for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }

309: static PetscErrorCode VecConjugate_Nest(Vec x)
310: {
311:   Vec_Nest *bx = (Vec_Nest *)x->data;
312:   PetscInt  j, nr;

314:   PetscFunctionBegin;
315:   nr = bx->nb;
316:   for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
321: {
322:   Vec_Nest *bx = (Vec_Nest *)x->data;
323:   Vec_Nest *by = (Vec_Nest *)y->data;
324:   PetscInt  i, nr;

326:   PetscFunctionBegin;
327:   VecNestCheckCompatible2(x, 1, y, 2);
328:   nr = bx->nb;
329:   for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
330:   PetscFunctionReturn(PETSC_SUCCESS);
331: }

333: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
334: {
335:   Vec_Nest *bx = (Vec_Nest *)x->data;
336:   Vec_Nest *by = (Vec_Nest *)y->data;
337:   Vec_Nest *bw = (Vec_Nest *)w->data;
338:   PetscInt  i, nr;

340:   PetscFunctionBegin;
341:   VecNestCheckCompatible3(w, 1, x, 3, y, 4);
342:   nr = bx->nb;
343:   for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
344:   PetscFunctionReturn(PETSC_SUCCESS);
345: }

347: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
348: {
349:   PetscInt  i, lp = -1, lw = -1;
350:   PetscReal lv;
351:   Vec_Nest *bx = (Vec_Nest *)x->data;

353:   PetscFunctionBegin;
354:   *v = PETSC_MAX_REAL;
355:   for (i = 0; i < bx->nb; i++) {
356:     PetscInt tp;
357:     PetscCall(VecMin(bx->v[i], &tp, &lv));
358:     if (lv < *v) {
359:       *v = lv;
360:       lw = i;
361:       lp = tp;
362:     }
363:   }
364:   if (p && lw > -1) {
365:     PetscInt        st, en;
366:     const PetscInt *idxs;

368:     *p = -1;
369:     PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
370:     if (lp >= st && lp < en) {
371:       PetscCall(ISGetIndices(bx->is[lw], &idxs));
372:       *p = idxs[lp - st];
373:       PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
374:     }
375:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
376:   }
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

380: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
381: {
382:   PetscInt  i, lp = -1, lw = -1;
383:   PetscReal lv;
384:   Vec_Nest *bx = (Vec_Nest *)x->data;

386:   PetscFunctionBegin;
387:   *v = PETSC_MIN_REAL;
388:   for (i = 0; i < bx->nb; i++) {
389:     PetscInt tp;
390:     PetscCall(VecMax(bx->v[i], &tp, &lv));
391:     if (lv > *v) {
392:       *v = lv;
393:       lw = i;
394:       lp = tp;
395:     }
396:   }
397:   if (p && lw > -1) {
398:     PetscInt        st, en;
399:     const PetscInt *idxs;

401:     *p = -1;
402:     PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
403:     if (lp >= st && lp < en) {
404:       PetscCall(ISGetIndices(bx->is[lw], &idxs));
405:       *p = idxs[lp - st];
406:       PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
407:     }
408:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
409:   }
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
414: {
415:   Vec_Nest *bx = (Vec_Nest *)x->data;
416:   PetscBool isascii;
417:   PetscInt  i;

419:   PetscFunctionBegin;
420:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
421:   if (isascii) {
422:     PetscCall(PetscViewerASCIIPushTab(viewer));
423:     PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ",  structure: \n", bx->nb));
424:     for (i = 0; i < bx->nb; i++) {
425:       VecType  type;
426:       char     name[256] = "", prefix[256] = "";
427:       PetscInt NR;

429:       PetscCall(VecGetSize(bx->v[i], &NR));
430:       PetscCall(VecGetType(bx->v[i], &type));
431:       if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
432:       if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));

434:       PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n", i, name, prefix, type, NR));

436:       PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
437:       PetscCall(VecView(bx->v[i], viewer));
438:       PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
439:     }
440:     PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
441:   }
442:   PetscFunctionReturn(PETSC_SUCCESS);
443: }

445: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
446: {
447:   Vec_Nest *bx;
448:   PetscInt  size, i, nr;
449:   PetscBool isnest;

451:   PetscFunctionBegin;
452:   PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
453:   if (!isnest) {
454:     /* Not nest */
455:     if (globalsize) PetscCall(VecGetSize(x, &size));
456:     else PetscCall(VecGetLocalSize(x, &size));
457:     *L = *L + size;
458:     PetscFunctionReturn(PETSC_SUCCESS);
459:   }

461:   /* Otherwise we have a nest */
462:   bx = (Vec_Nest *)x->data;
463:   nr = bx->nb;

465:   /* now descend recursively */
466:   for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: /* Returns the sum of the global size of all the constituent vectors in the nest */
471: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
472: {
473:   PetscFunctionBegin;
474:   *N = x->map->N;
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

478: /* Returns the sum of the local size of all the constituent vectors in the nest */
479: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
480: {
481:   PetscFunctionBegin;
482:   *n = x->map->n;
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
487: {
488:   Vec_Nest *bx = (Vec_Nest *)x->data;
489:   Vec_Nest *by = (Vec_Nest *)y->data;
490:   PetscInt  i, nr;
491:   PetscReal local_max, m;

493:   PetscFunctionBegin;
494:   VecNestCheckCompatible2(x, 1, y, 2);
495:   nr = bx->nb;
496:   m  = 0.0;
497:   for (i = 0; i < nr; i++) {
498:     PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
499:     if (local_max > m) m = local_max;
500:   }
501:   *max = m;
502:   PetscFunctionReturn(PETSC_SUCCESS);
503: }

505: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
506: {
507:   Vec_Nest *bx = (Vec_Nest *)X->data;
508:   PetscInt  i;

510:   PetscFunctionBegin;
511:   *x = NULL;
512:   for (i = 0; i < bx->nb; i++) {
513:     PetscBool issame = PETSC_FALSE;
514:     PetscCall(ISEqual(is, bx->is[i], &issame));
515:     if (issame) {
516:       *x = bx->v[i];
517:       PetscCall(PetscObjectReference((PetscObject)(*x)));
518:       break;
519:     }
520:   }
521:   PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
522:   PetscFunctionReturn(PETSC_SUCCESS);
523: }

525: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
526: {
527:   PetscFunctionBegin;
528:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
529:   PetscCall(VecDestroy(x));
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
534: {
535:   Vec_Nest *bx = (Vec_Nest *)X->data;
536:   PetscInt  i, m, rstart, rend;

538:   PetscFunctionBegin;
539:   PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
540:   PetscCall(VecGetLocalSize(X, &m));
541:   PetscCall(PetscMalloc1(m, x));
542:   for (i = 0; i < bx->nb; i++) {
543:     Vec                subvec = bx->v[i];
544:     IS                 isy    = bx->is[i];
545:     PetscInt           j, sm;
546:     const PetscInt    *ixy;
547:     const PetscScalar *y;
548:     PetscCall(VecGetLocalSize(subvec, &sm));
549:     PetscCall(VecGetArrayRead(subvec, &y));
550:     PetscCall(ISGetIndices(isy, &ixy));
551:     for (j = 0; j < sm; j++) {
552:       PetscInt ix = ixy[j];
553:       PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
554:       (*x)[ix - rstart] = y[j];
555:     }
556:     PetscCall(ISRestoreIndices(isy, &ixy));
557:     PetscCall(VecRestoreArrayRead(subvec, &y));
558:   }
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
563: {
564:   Vec_Nest *bx = (Vec_Nest *)X->data;
565:   PetscInt  i, m, rstart, rend;

567:   PetscFunctionBegin;
568:   PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
569:   PetscCall(VecGetLocalSize(X, &m));
570:   for (i = 0; i < bx->nb; i++) {
571:     Vec             subvec = bx->v[i];
572:     IS              isy    = bx->is[i];
573:     PetscInt        j, sm;
574:     const PetscInt *ixy;
575:     PetscScalar    *y;
576:     PetscCall(VecGetLocalSize(subvec, &sm));
577:     PetscCall(VecGetArray(subvec, &y));
578:     PetscCall(ISGetIndices(isy, &ixy));
579:     for (j = 0; j < sm; j++) {
580:       PetscInt ix = ixy[j];
581:       PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
582:       y[j] = (*x)[ix - rstart];
583:     }
584:     PetscCall(ISRestoreIndices(isy, &ixy));
585:     PetscCall(VecRestoreArray(subvec, &y));
586:   }
587:   PetscCall(PetscFree(*x));
588:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
593: {
594:   PetscFunctionBegin;
595:   PetscCall(PetscFree(*x));
596:   PetscFunctionReturn(PETSC_SUCCESS);
597: }

599: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
600: {
601:   PetscFunctionBegin;
602:   PetscCheck(nx <= 0, PetscObjectComm((PetscObject)(*X)), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
603:   PetscFunctionReturn(PETSC_SUCCESS);
604: }

606: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
607: {
608:   Vec      *ww;
609:   IS       *wis;
610:   Vec_Nest *bv = (Vec_Nest *)v->data;
611:   PetscInt  i;

613:   PetscFunctionBegin;
614:   PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
615:   for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
616:   for (i = 0; i < bv->nb; i++) {
617:     PetscInt off;

619:     PetscCall(VecGetOwnershipRange(v, &off, NULL));
620:     PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
621:     PetscCall(ISShift(wis[i], -off, wis[i]));
622:   }
623:   PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
624:   for (i = 0; i < bv->nb; i++) {
625:     PetscCall(VecDestroy(&ww[i]));
626:     PetscCall(ISDestroy(&wis[i]));
627:   }
628:   PetscCall(PetscFree2(ww, wis));
629:   PetscFunctionReturn(PETSC_SUCCESS);
630: }

632: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
633: {
634:   Vec_Nest *bv = (Vec_Nest *)v->data;
635:   Vec_Nest *bw = (Vec_Nest *)w->data;
636:   PetscInt  i;

638:   PetscFunctionBegin;
639:   PetscCheckSameType(v, 1, w, 2);
640:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
641:   for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
646: {
647:   Vec_Nest *bv = (Vec_Nest *)v->data;
648:   Vec_Nest *bw = (Vec_Nest *)w->data;
649:   PetscInt  i;

651:   PetscFunctionBegin;
652:   PetscCheckSameType(v, 1, w, 2);
653:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
654:   for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
655:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
656:   PetscFunctionReturn(PETSC_SUCCESS);
657: }

659: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
660: {
661:   Vec_Nest *bv = (Vec_Nest *)v->data;
662:   Vec_Nest *bw = (Vec_Nest *)w->data;
663:   PetscInt  i;

665:   PetscFunctionBegin;
666:   PetscCheckSameType(v, 1, w, 2);
667:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
668:   for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
669:   PetscFunctionReturn(PETSC_SUCCESS);
670: }

672: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
673: {
674:   Vec_Nest *bv = (Vec_Nest *)v->data;
675:   Vec_Nest *bw = (Vec_Nest *)w->data;
676:   PetscInt  i;

678:   PetscFunctionBegin;
679:   PetscCheckSameType(v, 1, w, 2);
680:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
681:   for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
682:   PetscFunctionReturn(PETSC_SUCCESS);
683: }

685: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
686: {
687:   Vec_Nest *bv = (Vec_Nest *)v->data;
688:   PetscInt  i;

690:   PetscFunctionBegin;
691:   for (i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }

695: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
696: {
697:   PetscFunctionBegin;
698:   ops->duplicate               = VecDuplicate_Nest;
699:   ops->duplicatevecs           = VecDuplicateVecs_Default;
700:   ops->destroyvecs             = VecDestroyVecs_Default;
701:   ops->dot                     = VecDot_Nest;
702:   ops->mdot                    = VecMDot_Nest;
703:   ops->norm                    = VecNorm_Nest;
704:   ops->tdot                    = VecTDot_Nest;
705:   ops->mtdot                   = VecMTDot_Nest;
706:   ops->scale                   = VecScale_Nest;
707:   ops->copy                    = VecCopy_Nest;
708:   ops->set                     = VecSet_Nest;
709:   ops->swap                    = VecSwap_Nest;
710:   ops->axpy                    = VecAXPY_Nest;
711:   ops->axpby                   = VecAXPBY_Nest;
712:   ops->maxpy                   = VecMAXPY_Nest;
713:   ops->aypx                    = VecAYPX_Nest;
714:   ops->waxpy                   = VecWAXPY_Nest;
715:   ops->axpbypcz                = NULL;
716:   ops->pointwisemult           = VecPointwiseMult_Nest;
717:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
718:   ops->setvalues               = NULL;
719:   ops->assemblybegin           = VecAssemblyBegin_Nest;
720:   ops->assemblyend             = VecAssemblyEnd_Nest;
721:   ops->getarray                = VecGetArray_Nest;
722:   ops->getsize                 = VecGetSize_Nest;
723:   ops->getlocalsize            = VecGetLocalSize_Nest;
724:   ops->restorearray            = VecRestoreArray_Nest;
725:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
726:   ops->max                     = VecMax_Nest;
727:   ops->min                     = VecMin_Nest;
728:   ops->setrandom               = NULL;
729:   ops->setoption               = NULL;
730:   ops->setvaluesblocked        = NULL;
731:   ops->destroy                 = VecDestroy_Nest;
732:   ops->view                    = VecView_Nest;
733:   ops->placearray              = NULL;
734:   ops->replacearray            = NULL;
735:   ops->dot_local               = VecDot_Nest;
736:   ops->tdot_local              = VecTDot_Nest;
737:   ops->norm_local              = VecNorm_Nest;
738:   ops->mdot_local              = VecMDot_Nest;
739:   ops->mtdot_local             = VecMTDot_Nest;
740:   ops->load                    = NULL;
741:   ops->reciprocal              = VecReciprocal_Nest;
742:   ops->conjugate               = VecConjugate_Nest;
743:   ops->setlocaltoglobalmapping = NULL;
744:   ops->setvalueslocal          = NULL;
745:   ops->resetarray              = NULL;
746:   ops->setfromoptions          = NULL;
747:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
748:   ops->load                    = NULL;
749:   ops->pointwisemax            = NULL;
750:   ops->pointwisemaxabs         = NULL;
751:   ops->pointwisemin            = NULL;
752:   ops->getvalues               = NULL;
753:   ops->sqrt                    = NULL;
754:   ops->abs                     = NULL;
755:   ops->exp                     = NULL;
756:   ops->shift                   = NULL;
757:   ops->create                  = NULL;
758:   ops->stridegather            = NULL;
759:   ops->stridescatter           = NULL;
760:   ops->dotnorm2                = VecDotNorm2_Nest;
761:   ops->getsubvector            = VecGetSubVector_Nest;
762:   ops->restoresubvector        = VecRestoreSubVector_Nest;
763:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
764:   ops->concatenate             = VecConcatenate_Nest;
765:   ops->createlocalvector       = VecCreateLocalVector_Nest;
766:   ops->getlocalvector          = VecGetLocalVector_Nest;
767:   ops->getlocalvectorread      = VecGetLocalVectorRead_Nest;
768:   ops->restorelocalvector      = VecRestoreLocalVector_Nest;
769:   ops->restorelocalvectorread  = VecRestoreLocalVectorRead_Nest;
770:   ops->setrandom               = VecSetRandom_Nest;
771:   PetscFunctionReturn(PETSC_SUCCESS);
772: }

774: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
775: {
776:   Vec_Nest *b = (Vec_Nest *)x->data;
777:   PetscInt  i;
778:   PetscInt  row;

780:   PetscFunctionBegin;
781:   if (!m) PetscFunctionReturn(PETSC_SUCCESS);
782:   for (i = 0; i < m; i++) {
783:     row = idxm[i];
784:     PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
785:     vec[i] = b->v[row];
786:   }
787:   PetscFunctionReturn(PETSC_SUCCESS);
788: }

790: PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
791: {
792:   PetscFunctionBegin;
793:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
794:   PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
795:   PetscFunctionReturn(PETSC_SUCCESS);
796: }

798: /*@
799:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

801:  Not Collective

803:  Input Parameters:
804: +  X  - nest vector
805: -  idxm - index of the vector within the nest

807:  Output Parameter:
808: .  sx - vector at index `idxm` within the nest

810:  Level: developer

812: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
813: @*/
814: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
815: {
816:   PetscFunctionBegin;
817:   PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
818:   PetscFunctionReturn(PETSC_SUCCESS);
819: }

821: PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
822: {
823:   Vec_Nest *b = (Vec_Nest *)X->data;

825:   PetscFunctionBegin;
826:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
827:   if (N) *N = b->nb;
828:   if (sx) *sx = b->v;
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: /*@C
833:  VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.

835:  Not Collective

837:  Input Parameter:
838: .  X  - nest vector

840:  Output Parameters:
841: +  N - number of nested vecs
842: -  sx - array of vectors

844:  Level: developer

846:  Note:
847:  The user should not free the array `sx`.

849:  Fortran Note:
850:  The caller must allocate the array to hold the subvectors.

852: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
853: @*/
854: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec **sx)
855: {
856:   PetscFunctionBegin;
857:   PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
858:   PetscFunctionReturn(PETSC_SUCCESS);
859: }

861: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
862: {
863:   Vec_Nest *bx = (Vec_Nest *)X->data;
864:   PetscInt  i, offset = 0, n = 0, bs;
865:   IS        is;
866:   PetscBool issame = PETSC_FALSE;
867:   PetscInt  N      = 0;

869:   PetscFunctionBegin;
870:   /* check if idxm < bx->nb */
871:   PetscCheck(idxm < bx->nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, idxm, bx->nb);

873:   PetscCall(VecDestroy(&bx->v[idxm]));      /* destroy the existing vector */
874:   PetscCall(VecDuplicate(x, &bx->v[idxm])); /* duplicate the layout of given vector */
875:   PetscCall(VecCopy(x, bx->v[idxm]));       /* copy the contents of the given vector */

877:   /* check if we need to update the IS for the block */
878:   offset = X->map->rstart;
879:   for (i = 0; i < idxm; i++) {
880:     n = 0;
881:     PetscCall(VecGetLocalSize(bx->v[i], &n));
882:     offset += n;
883:   }

885:   /* get the local size and block size */
886:   PetscCall(VecGetLocalSize(x, &n));
887:   PetscCall(VecGetBlockSize(x, &bs));

889:   /* create the new IS */
890:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
891:   PetscCall(ISSetBlockSize(is, bs));

893:   /* check if they are equal */
894:   PetscCall(ISEqual(is, bx->is[idxm], &issame));

896:   if (!issame) {
897:     /* The IS of given vector has a different layout compared to the existing block vector.
898:      Destroy the existing reference and update the IS. */
899:     PetscCall(ISDestroy(&bx->is[idxm]));
900:     PetscCall(ISDuplicate(is, &bx->is[idxm]));
901:     PetscCall(ISCopy(is, bx->is[idxm]));

903:     offset += n;
904:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
905:     for (i = idxm + 1; i < bx->nb; i++) {
906:       /* get the local size and block size */
907:       PetscCall(VecGetLocalSize(bx->v[i], &n));
908:       PetscCall(VecGetBlockSize(bx->v[i], &bs));

910:       /* destroy the old and create the new IS */
911:       PetscCall(ISDestroy(&bx->is[i]));
912:       PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
913:       PetscCall(ISSetBlockSize(bx->is[i], bs));

915:       offset += n;
916:     }

918:     n = 0;
919:     PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
920:     PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
921:     PetscCall(PetscLayoutSetSize(X->map, N));
922:     PetscCall(PetscLayoutSetLocalSize(X->map, n));
923:   }

925:   PetscCall(ISDestroy(&is));
926:   PetscFunctionReturn(PETSC_SUCCESS);
927: }

929: PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
930: {
931:   PetscFunctionBegin;
932:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
933:   PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
934:   PetscFunctionReturn(PETSC_SUCCESS);
935: }

937: /*@
938:    VecNestSetSubVec - Set a single component vector in a nest vector at specified index.

940:    Not Collective

942:    Input Parameters:
943: +  X  - nest vector
944: .  idxm - index of the vector within the nest vector
945: -  sx - vector at index `idxm` within the nest vector

947:    Level: developer

949:    Note:
950:    The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.

952: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
953: @*/
954: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
955: {
956:   PetscFunctionBegin;
957:   PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
958:   PetscFunctionReturn(PETSC_SUCCESS);
959: }

961: PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
962: {
963:   PetscInt i;

965:   PetscFunctionBegin;
966:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
967:   for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
968:   PetscFunctionReturn(PETSC_SUCCESS);
969: }

971: /*@C
972:    VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.

974:    Not Collective

976:    Input Parameters:
977: +  X  - nest vector
978: .  N - number of component vecs in `sx`
979: .  idxm - indices of component vectors that are to be replaced
980: -  sx - array of vectors

982:    Level: developer

984:    Note:
985:    The components in the vector array `sx` do not have to be of the same size as corresponding
986:    components in `X`. The user can also free the array `sx` after the call.

988: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
989: @*/
990: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
991: {
992:   PetscFunctionBegin;
993:   PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
994:   PetscFunctionReturn(PETSC_SUCCESS);
995: }

997: PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
998: {
999:   Vec_Nest *b = (Vec_Nest *)X->data;

1001:   PetscFunctionBegin;
1002:   *N = b->nb;
1003:   PetscFunctionReturn(PETSC_SUCCESS);
1004: }

1006: /*@
1007:  VecNestGetSize - Returns the size of the nest vector.

1009:  Not Collective

1011:  Input Parameter:
1012: .  X  - nest vector

1014:  Output Parameter:
1015: .  N - number of nested vecs

1017:  Level: developer

1019: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1020: @*/
1021: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1022: {
1023:   PetscFunctionBegin;
1026:   PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1027:   PetscFunctionReturn(PETSC_SUCCESS);
1028: }

1030: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1031: {
1032:   Vec_Nest *ctx = (Vec_Nest *)V->data;
1033:   PetscInt  i;

1035:   PetscFunctionBegin;
1036:   if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);

1038:   ctx->nb = nb;
1039:   PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");

1041:   /* Create space */
1042:   PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1043:   for (i = 0; i < ctx->nb; i++) {
1044:     ctx->v[i] = x[i];
1045:     PetscCall(PetscObjectReference((PetscObject)x[i]));
1046:     /* Do not allocate memory for internal sub blocks */
1047:   }

1049:   PetscCall(PetscMalloc1(ctx->nb, &ctx->is));

1051:   ctx->setup_called = PETSC_TRUE;
1052:   PetscFunctionReturn(PETSC_SUCCESS);
1053: }

1055: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1056: {
1057:   Vec_Nest *ctx = (Vec_Nest *)V->data;
1058:   PetscInt  i, offset, m, n, M, N;

1060:   PetscFunctionBegin;
1061:   if (is) { /* Do some consistency checks and reference the is */
1062:     offset = V->map->rstart;
1063:     for (i = 0; i < ctx->nb; i++) {
1064:       PetscCall(ISGetSize(is[i], &M));
1065:       PetscCall(VecGetSize(ctx->v[i], &N));
1066:       PetscCheck(M == N, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of size %" PetscInt_FMT " is not compatible with Vec of size %" PetscInt_FMT, i, M, N);
1067:       PetscCall(ISGetLocalSize(is[i], &m));
1068:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
1069:       PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of local size %" PetscInt_FMT " is not compatible with Vec of local size %" PetscInt_FMT, i, m, n);
1070:       PetscCall(PetscObjectReference((PetscObject)is[i]));
1071:       ctx->is[i] = is[i];
1072:       offset += n;
1073:     }
1074:   } else { /* Create a contiguous ISStride for each entry */
1075:     offset = V->map->rstart;
1076:     for (i = 0; i < ctx->nb; i++) {
1077:       PetscInt bs;
1078:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
1079:       PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1080:       PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1081:       PetscCall(ISSetBlockSize(ctx->is[i], bs));
1082:       offset += n;
1083:     }
1084:   }
1085:   PetscFunctionReturn(PETSC_SUCCESS);
1086: }

1088: /*MC
1089:   VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.

1091:   Level: intermediate

1093:   Notes:
1094:   This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1095:   It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.

1097: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1098: M*/

1100: /*@C
1101:    VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately

1103:    Collective

1105:    Input Parameters:
1106: +  comm - Communicator for the new `Vec`
1107: .  nb - number of nested blocks
1108: .  is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1109: -  x - array of `nb` sub-vectors

1111:    Output Parameter:
1112: .  Y - new vector

1114:    Level: advanced

1116: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`, `VECNEST`
1117: @*/
1118: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1119: {
1120:   Vec       V;
1121:   Vec_Nest *s;
1122:   PetscInt  n, N;

1124:   PetscFunctionBegin;
1125:   PetscCall(VecCreate(comm, &V));
1126:   PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));

1128:   /* allocate and set pointer for implementation data */
1129:   PetscCall(PetscNew(&s));
1130:   V->data         = (void *)s;
1131:   s->setup_called = PETSC_FALSE;
1132:   s->nb           = -1;
1133:   s->v            = NULL;

1135:   PetscCall(VecSetUp_Nest_Private(V, nb, x));

1137:   n = N = 0;
1138:   PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1139:   PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1140:   PetscCall(PetscLayoutSetSize(V->map, N));
1141:   PetscCall(PetscLayoutSetLocalSize(V->map, n));
1142:   PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1143:   PetscCall(PetscLayoutSetUp(V->map));

1145:   PetscCall(VecSetUp_NestIS_Private(V, nb, is));

1147:   PetscCall(VecNestSetOps_Private(V->ops));
1148:   V->petscnative = PETSC_FALSE;

1150:   /* expose block api's */
1151:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1152:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1153:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1154:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1155:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));

1157:   *Y = V;
1158:   PetscFunctionReturn(PETSC_SUCCESS);
1159: }