Actual source code: vecnest.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

693: static PetscErrorCode VecErrorWeightedNorms_Nest(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
694: {
695:   Vec_Nest *bu = (Vec_Nest *)U->data;
696:   Vec_Nest *by = (Vec_Nest *)Y->data;
697:   Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
698:   Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
699:   Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;

701:   PetscFunctionBegin;
702:   VecNestCheckCompatible2(U, 1, Y, 2);
703:   if (E) VecNestCheckCompatible2(U, 1, E, 3);
704:   if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
705:   if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
706:   *norm      = 0.0;
707:   *norma     = 0.0;
708:   *normr     = 0.0;
709:   *norm_loc  = 0;
710:   *norma_loc = 0;
711:   *normr_loc = 0;
712:   for (PetscInt i = 0; i < bu->nb; i++) {
713:     PetscReal n, na, nr;
714:     PetscInt  n_loc, na_loc, nr_loc;

716:     PetscCall(VecErrorWeightedNorms(bu->v[i], by->v[i], be ? be->v[i] : NULL, wnormtype, atol, ba ? ba->v[i] : NULL, rtol, br ? br->v[i] : NULL, ignore_max, &n, &n_loc, &na, &na_loc, &nr, &nr_loc));
717:     if (wnormtype == NORM_INFINITY) {
718:       *norm  = PetscMax(*norm, n);
719:       *norma = PetscMax(*norma, na);
720:       *normr = PetscMax(*normr, nr);
721:     } else {
722:       *norm += PetscSqr(n);
723:       *norma += PetscSqr(na);
724:       *normr += PetscSqr(nr);
725:     }
726:     *norm_loc += n_loc;
727:     *norma_loc += na_loc;
728:     *normr_loc += nr_loc;
729:   }
730:   if (wnormtype == NORM_2) {
731:     *norm  = PetscSqrtReal(*norm);
732:     *norma = PetscSqrtReal(*norma);
733:     *normr = PetscSqrtReal(*normr);
734:   }
735:   PetscFunctionReturn(PETSC_SUCCESS);
736: }

738: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
739: {
740:   PetscFunctionBegin;
741:   ops->duplicate               = VecDuplicate_Nest;
742:   ops->duplicatevecs           = VecDuplicateVecs_Default;
743:   ops->destroyvecs             = VecDestroyVecs_Default;
744:   ops->dot                     = VecDot_Nest;
745:   ops->mdot                    = VecMDot_Nest;
746:   ops->norm                    = VecNorm_Nest;
747:   ops->tdot                    = VecTDot_Nest;
748:   ops->mtdot                   = VecMTDot_Nest;
749:   ops->scale                   = VecScale_Nest;
750:   ops->copy                    = VecCopy_Nest;
751:   ops->set                     = VecSet_Nest;
752:   ops->swap                    = VecSwap_Nest;
753:   ops->axpy                    = VecAXPY_Nest;
754:   ops->axpby                   = VecAXPBY_Nest;
755:   ops->maxpy                   = VecMAXPY_Nest;
756:   ops->aypx                    = VecAYPX_Nest;
757:   ops->waxpy                   = VecWAXPY_Nest;
758:   ops->axpbypcz                = NULL;
759:   ops->pointwisemult           = VecPointwiseMult_Nest;
760:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
761:   ops->setvalues               = NULL;
762:   ops->assemblybegin           = VecAssemblyBegin_Nest;
763:   ops->assemblyend             = VecAssemblyEnd_Nest;
764:   ops->getarray                = VecGetArray_Nest;
765:   ops->getsize                 = VecGetSize_Nest;
766:   ops->getlocalsize            = VecGetLocalSize_Nest;
767:   ops->restorearray            = VecRestoreArray_Nest;
768:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
769:   ops->max                     = VecMax_Nest;
770:   ops->min                     = VecMin_Nest;
771:   ops->setrandom               = NULL;
772:   ops->setoption               = NULL;
773:   ops->setvaluesblocked        = NULL;
774:   ops->destroy                 = VecDestroy_Nest;
775:   ops->view                    = VecView_Nest;
776:   ops->placearray              = NULL;
777:   ops->replacearray            = NULL;
778:   ops->dot_local               = VecDot_Nest;
779:   ops->tdot_local              = VecTDot_Nest;
780:   ops->norm_local              = VecNorm_Nest;
781:   ops->mdot_local              = VecMDot_Nest;
782:   ops->mtdot_local             = VecMTDot_Nest;
783:   ops->load                    = NULL;
784:   ops->reciprocal              = VecReciprocal_Nest;
785:   ops->conjugate               = VecConjugate_Nest;
786:   ops->setlocaltoglobalmapping = NULL;
787:   ops->setvalueslocal          = NULL;
788:   ops->resetarray              = NULL;
789:   ops->setfromoptions          = NULL;
790:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
791:   ops->load                    = NULL;
792:   ops->pointwisemax            = NULL;
793:   ops->pointwisemaxabs         = NULL;
794:   ops->pointwisemin            = NULL;
795:   ops->getvalues               = NULL;
796:   ops->sqrt                    = NULL;
797:   ops->abs                     = NULL;
798:   ops->exp                     = NULL;
799:   ops->shift                   = NULL;
800:   ops->create                  = NULL;
801:   ops->stridegather            = NULL;
802:   ops->stridescatter           = NULL;
803:   ops->dotnorm2                = VecDotNorm2_Nest;
804:   ops->getsubvector            = VecGetSubVector_Nest;
805:   ops->restoresubvector        = VecRestoreSubVector_Nest;
806:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
807:   ops->concatenate             = VecConcatenate_Nest;
808:   ops->createlocalvector       = VecCreateLocalVector_Nest;
809:   ops->getlocalvector          = VecGetLocalVector_Nest;
810:   ops->getlocalvectorread      = VecGetLocalVectorRead_Nest;
811:   ops->restorelocalvector      = VecRestoreLocalVector_Nest;
812:   ops->restorelocalvectorread  = VecRestoreLocalVectorRead_Nest;
813:   ops->setrandom               = VecSetRandom_Nest;
814:   ops->errorwnorm              = VecErrorWeightedNorms_Nest;
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }

818: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
819: {
820:   Vec_Nest *b = (Vec_Nest *)x->data;
821:   PetscInt  i;
822:   PetscInt  row;

824:   PetscFunctionBegin;
825:   if (!m) PetscFunctionReturn(PETSC_SUCCESS);
826:   for (i = 0; i < m; i++) {
827:     row = idxm[i];
828:     PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
829:     vec[i] = b->v[row];
830:   }
831:   PetscFunctionReturn(PETSC_SUCCESS);
832: }

834: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
835: {
836:   PetscFunctionBegin;
837:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
838:   PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
839:   PetscFunctionReturn(PETSC_SUCCESS);
840: }

842: /*@
843:   VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

845:   Not Collective

847:   Input Parameters:
848: + X    - nest vector
849: - idxm - index of the vector within the nest

851:   Output Parameter:
852: . sx - vector at index `idxm` within the nest

854:   Level: developer

856: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
857: @*/
858: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
859: {
860:   PetscFunctionBegin;
861:   PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
862:   PetscFunctionReturn(PETSC_SUCCESS);
863: }

865: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
866: {
867:   Vec_Nest *b = (Vec_Nest *)X->data;

869:   PetscFunctionBegin;
870:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
871:   if (N) *N = b->nb;
872:   if (sx) *sx = b->v;
873:   PetscFunctionReturn(PETSC_SUCCESS);
874: }

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

879:   Not Collective

881:   Input Parameter:
882: . X - nest vector

884:   Output Parameters:
885: + N  - number of nested vecs
886: - sx - array of vectors

888:   Level: developer

890:   Note:
891:   The user should not free the array `sx`.

893:   Fortran Notes:
894:   The caller must allocate the array to hold the subvectors.

896: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
897: @*/
898: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec **sx)
899: {
900:   PetscFunctionBegin;
901:   PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
902:   PetscFunctionReturn(PETSC_SUCCESS);
903: }

905: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
906: {
907:   Vec_Nest *bx = (Vec_Nest *)X->data;
908:   PetscInt  i, offset = 0, n = 0, bs;
909:   IS        is;
910:   PetscBool issame = PETSC_FALSE;
911:   PetscInt  N      = 0;

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

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

921:   /* check if we need to update the IS for the block */
922:   offset = X->map->rstart;
923:   for (i = 0; i < idxm; i++) {
924:     n = 0;
925:     PetscCall(VecGetLocalSize(bx->v[i], &n));
926:     offset += n;
927:   }

929:   /* get the local size and block size */
930:   PetscCall(VecGetLocalSize(x, &n));
931:   PetscCall(VecGetBlockSize(x, &bs));

933:   /* create the new IS */
934:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
935:   PetscCall(ISSetBlockSize(is, bs));

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

940:   if (!issame) {
941:     /* The IS of given vector has a different layout compared to the existing block vector.
942:      Destroy the existing reference and update the IS. */
943:     PetscCall(ISDestroy(&bx->is[idxm]));
944:     PetscCall(ISDuplicate(is, &bx->is[idxm]));
945:     PetscCall(ISCopy(is, bx->is[idxm]));

947:     offset += n;
948:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
949:     for (i = idxm + 1; i < bx->nb; i++) {
950:       /* get the local size and block size */
951:       PetscCall(VecGetLocalSize(bx->v[i], &n));
952:       PetscCall(VecGetBlockSize(bx->v[i], &bs));

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

959:       offset += n;
960:     }

962:     n = 0;
963:     PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
964:     PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
965:     PetscCall(PetscLayoutSetSize(X->map, N));
966:     PetscCall(PetscLayoutSetLocalSize(X->map, n));
967:   }

969:   PetscCall(ISDestroy(&is));
970:   PetscFunctionReturn(PETSC_SUCCESS);
971: }

973: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
974: {
975:   PetscFunctionBegin;
976:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
977:   PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
978:   PetscFunctionReturn(PETSC_SUCCESS);
979: }

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

984:   Not Collective

986:   Input Parameters:
987: + X    - nest vector
988: . idxm - index of the vector within the nest vector
989: - sx   - vector at index `idxm` within the nest vector

991:   Level: developer

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

996: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
997: @*/
998: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
999: {
1000:   PetscFunctionBegin;
1001:   PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

1005: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1006: {
1007:   PetscInt i;

1009:   PetscFunctionBegin;
1010:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
1011:   for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1012:   PetscFunctionReturn(PETSC_SUCCESS);
1013: }

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

1018:   Not Collective

1020:   Input Parameters:
1021: + X    - nest vector
1022: . N    - number of component vecs in `sx`
1023: . idxm - indices of component vectors that are to be replaced
1024: - sx   - array of vectors

1026:   Level: developer

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

1032: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1033: @*/
1034: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1035: {
1036:   PetscFunctionBegin;
1037:   PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

1041: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1042: {
1043:   Vec_Nest *b = (Vec_Nest *)X->data;

1045:   PetscFunctionBegin;
1046:   *N = b->nb;
1047:   PetscFunctionReturn(PETSC_SUCCESS);
1048: }

1050: /*@
1051:   VecNestGetSize - Returns the size of the nest vector.

1053:   Not Collective

1055:   Input Parameter:
1056: . X - nest vector

1058:   Output Parameter:
1059: . N - number of nested vecs

1061:   Level: developer

1063: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1064: @*/
1065: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1066: {
1067:   PetscFunctionBegin;
1069:   PetscAssertPointer(N, 2);
1070:   PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1071:   PetscFunctionReturn(PETSC_SUCCESS);
1072: }

1074: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1075: {
1076:   Vec_Nest *ctx = (Vec_Nest *)V->data;
1077:   PetscInt  i;

1079:   PetscFunctionBegin;
1080:   if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);

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

1085:   /* Create space */
1086:   PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1087:   for (i = 0; i < ctx->nb; i++) {
1088:     ctx->v[i] = x[i];
1089:     PetscCall(PetscObjectReference((PetscObject)x[i]));
1090:     /* Do not allocate memory for internal sub blocks */
1091:   }

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

1095:   ctx->setup_called = PETSC_TRUE;
1096:   PetscFunctionReturn(PETSC_SUCCESS);
1097: }

1099: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1100: {
1101:   Vec_Nest *ctx = (Vec_Nest *)V->data;
1102:   PetscInt  i, offset, m, n, M, N;

1104:   PetscFunctionBegin;
1105:   (void)nb;
1106:   if (is) { /* Do some consistency checks and reference the is */
1107:     offset = V->map->rstart;
1108:     for (i = 0; i < ctx->nb; i++) {
1109:       PetscCall(ISGetSize(is[i], &M));
1110:       PetscCall(VecGetSize(ctx->v[i], &N));
1111:       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);
1112:       PetscCall(ISGetLocalSize(is[i], &m));
1113:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
1114:       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);
1115:       PetscCall(PetscObjectReference((PetscObject)is[i]));
1116:       ctx->is[i] = is[i];
1117:       offset += n;
1118:     }
1119:   } else { /* Create a contiguous ISStride for each entry */
1120:     offset = V->map->rstart;
1121:     for (i = 0; i < ctx->nb; i++) {
1122:       PetscInt bs;
1123:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
1124:       PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1125:       PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1126:       PetscCall(ISSetBlockSize(ctx->is[i], bs));
1127:       offset += n;
1128:     }
1129:   }
1130:   PetscFunctionReturn(PETSC_SUCCESS);
1131: }

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

1136:   Level: intermediate

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

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

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

1148:   Collective

1150:   Input Parameters:
1151: + comm - Communicator for the new `Vec`
1152: . nb   - number of nested blocks
1153: . is   - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1154: - x    - array of `nb` sub-vectors

1156:   Output Parameter:
1157: . Y - new vector

1159:   Level: advanced

1161: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1162: @*/
1163: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1164: {
1165:   Vec       V;
1166:   Vec_Nest *s;
1167:   PetscInt  n, N;

1169:   PetscFunctionBegin;
1170:   PetscCall(VecCreate(comm, &V));
1171:   PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));

1173:   /* allocate and set pointer for implementation data */
1174:   PetscCall(PetscNew(&s));
1175:   V->data         = (void *)s;
1176:   s->setup_called = PETSC_FALSE;
1177:   s->nb           = -1;
1178:   s->v            = NULL;

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

1182:   n = N = 0;
1183:   PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1184:   PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1185:   PetscCall(PetscLayoutSetSize(V->map, N));
1186:   PetscCall(PetscLayoutSetLocalSize(V->map, n));
1187:   PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1188:   PetscCall(PetscLayoutSetUp(V->map));

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

1192:   PetscCall(VecNestSetOps_Private(V->ops));
1193:   V->petscnative = PETSC_FALSE;

1195:   /* expose block api's */
1196:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1197:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1198:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1199:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1200:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));

1202:   *Y = V;
1203:   PetscFunctionReturn(PETSC_SUCCESS);
1204: }