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:   for (i = 0; i < vs->nb; i++) {
 12:     VecAssemblyBegin(vs->v[i]);
 13:   }
 14:   return 0;
 15: }

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

 22:   for (i = 0; i < vs->nb; i++) VecAssemblyEnd(vs->v[i]);
 23:   return 0;
 24: }

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

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

 43:   PetscFree(vs);
 44:   return 0;
 45: }

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

 54:   VecNestCheckCompatible2(x, 1, y, 2);
 55:   for (i = 0; i < bx->nb; i++) VecCopy(bx->v[i], by->v[i]);
 56:   return 0;
 57: }

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

 66:   PetscMalloc1(bx->nb, &sub);
 67:   for (i = 0; i < bx->nb; i++) VecDuplicate(bx->v[i], &sub[i]);
 68:   VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y);
 69:   for (i = 0; i < bx->nb; i++) VecDestroy(&sub[i]);
 70:   PetscFree(sub);
 71:   *y = Y;
 72:   return 0;
 73: }

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

 82:   nr   = bx->nb;
 83:   _val = 0.0;
 84:   for (i = 0; i < nr; i++) {
 85:     VecDot(bx->v[i], by->v[i], &x_dot_y);
 86:     _val = _val + x_dot_y;
 87:   }
 88:   *val = _val;
 89:   return 0;
 90: }

 92: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
 93: {
 94:   Vec_Nest   *bx = (Vec_Nest *)x->data;
 95:   Vec_Nest   *by = (Vec_Nest *)y->data;
 96:   PetscInt    i, nr;
 97:   PetscScalar x_dot_y, _val;

 99:   nr   = bx->nb;
100:   _val = 0.0;
101:   for (i = 0; i < nr; i++) {
102:     VecTDot(bx->v[i], by->v[i], &x_dot_y);
103:     _val = _val + x_dot_y;
104:   }
105:   *val = _val;
106:   return 0;
107: }

109: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
110: {
111:   Vec_Nest   *bx = (Vec_Nest *)x->data;
112:   Vec_Nest   *by = (Vec_Nest *)y->data;
113:   PetscInt    i, nr;
114:   PetscScalar x_dot_y, _dp, _nm;
115:   PetscReal   norm2_y;

117:   nr  = bx->nb;
118:   _dp = 0.0;
119:   _nm = 0.0;
120:   for (i = 0; i < nr; i++) {
121:     VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y);
122:     _dp += x_dot_y;
123:     _nm += norm2_y;
124:   }
125:   *dp = _dp;
126:   *nm = _nm;
127:   return 0;
128: }

130: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
131: {
132:   Vec_Nest *bx = (Vec_Nest *)x->data;
133:   Vec_Nest *by = (Vec_Nest *)y->data;
134:   PetscInt  i, nr;

136:   nr = bx->nb;
137:   for (i = 0; i < nr; i++) VecAXPY(by->v[i], alpha, bx->v[i]);
138:   return 0;
139: }

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

147:   nr = bx->nb;
148:   for (i = 0; i < nr; i++) VecAYPX(by->v[i], alpha, bx->v[i]);
149:   return 0;
150: }

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

158:   nr = bx->nb;
159:   for (i = 0; i < nr; i++) VecAXPBY(by->v[i], alpha, beta, bx->v[i]);
160:   return 0;
161: }

163: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
164: {
165:   Vec_Nest *bx = (Vec_Nest *)x->data;
166:   Vec_Nest *by = (Vec_Nest *)y->data;
167:   Vec_Nest *bz = (Vec_Nest *)z->data;
168:   PetscInt  i, nr;

170:   nr = bx->nb;
171:   for (i = 0; i < nr; i++) VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]);
172:   return 0;
173: }

175: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
176: {
177:   Vec_Nest *bx = (Vec_Nest *)x->data;
178:   PetscInt  i, nr;

180:   nr = bx->nb;
181:   for (i = 0; i < nr; i++) VecScale(bx->v[i], alpha);
182:   return 0;
183: }

185: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
186: {
187:   Vec_Nest *bx = (Vec_Nest *)x->data;
188:   Vec_Nest *by = (Vec_Nest *)y->data;
189:   Vec_Nest *bw = (Vec_Nest *)w->data;
190:   PetscInt  i, nr;

192:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);
193:   nr = bx->nb;
194:   for (i = 0; i < nr; i++) VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]);
195:   return 0;
196: }

198: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
199: {
200:   Vec_Nest *bx = (Vec_Nest *)x->data;
201:   Vec_Nest *by = (Vec_Nest *)y->data;
202:   Vec_Nest *bw = (Vec_Nest *)w->data;
203:   PetscInt  i, nr;

205:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);

207:   nr = bx->nb;
208:   for (i = 0; i < nr; i++) VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]);
209:   return 0;
210: }

212: static PetscErrorCode VecReciprocal_Nest(Vec x)
213: {
214:   Vec_Nest *bx = (Vec_Nest *)x->data;
215:   PetscInt  i, nr;

217:   nr = bx->nb;
218:   for (i = 0; i < nr; i++) VecReciprocal(bx->v[i]);
219:   return 0;
220: }

222: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
223: {
224:   Vec_Nest *bx = (Vec_Nest *)xin->data;
225:   PetscInt  i, nr;
226:   PetscReal z_i;
227:   PetscReal _z;

229:   nr = bx->nb;
230:   _z = 0.0;

232:   if (type == NORM_2) {
233:     PetscScalar dot;
234:     VecDot(xin, xin, &dot);
235:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
236:   } else if (type == NORM_1) {
237:     for (i = 0; i < nr; i++) {
238:       VecNorm(bx->v[i], type, &z_i);
239:       _z = _z + z_i;
240:     }
241:   } else if (type == NORM_INFINITY) {
242:     for (i = 0; i < nr; i++) {
243:       VecNorm(bx->v[i], type, &z_i);
244:       if (z_i > _z) _z = z_i;
245:     }
246:   }

248:   *z = _z;
249:   return 0;
250: }

252: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
253: {
254:   PetscInt v;

256:   for (v = 0; v < nv; v++) {
257:     /* Do axpy on each vector,v */
258:     VecAXPY(y, alpha[v], x[v]);
259:   }
260:   return 0;
261: }

263: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
264: {
265:   PetscInt j;

267:   for (j = 0; j < nv; j++) VecDot(x, y[j], &val[j]);
268:   return 0;
269: }

271: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
272: {
273:   PetscInt j;

275:   for (j = 0; j < nv; j++) VecTDot(x, y[j], &val[j]);
276:   return 0;
277: }

279: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
280: {
281:   Vec_Nest *bx = (Vec_Nest *)x->data;
282:   PetscInt  i, nr;

284:   nr = bx->nb;
285:   for (i = 0; i < nr; i++) VecSet(bx->v[i], alpha);
286:   return 0;
287: }

289: static PetscErrorCode VecConjugate_Nest(Vec x)
290: {
291:   Vec_Nest *bx = (Vec_Nest *)x->data;
292:   PetscInt  j, nr;

294:   nr = bx->nb;
295:   for (j = 0; j < nr; j++) VecConjugate(bx->v[j]);
296:   return 0;
297: }

299: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
300: {
301:   Vec_Nest *bx = (Vec_Nest *)x->data;
302:   Vec_Nest *by = (Vec_Nest *)y->data;
303:   PetscInt  i, nr;

305:   VecNestCheckCompatible2(x, 1, y, 2);
306:   nr = bx->nb;
307:   for (i = 0; i < nr; i++) VecSwap(bx->v[i], by->v[i]);
308:   return 0;
309: }

311: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
312: {
313:   Vec_Nest *bx = (Vec_Nest *)x->data;
314:   Vec_Nest *by = (Vec_Nest *)y->data;
315:   Vec_Nest *bw = (Vec_Nest *)w->data;
316:   PetscInt  i, nr;

318:   VecNestCheckCompatible3(w, 1, x, 3, y, 4);

320:   nr = bx->nb;
321:   for (i = 0; i < nr; i++) VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]);
322:   return 0;
323: }

325: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
326: {
327:   PetscInt  i, lp = -1, lw = -1;
328:   PetscReal lv;
329:   Vec_Nest *bx = (Vec_Nest *)x->data;

331:   *v = PETSC_MAX_REAL;
332:   for (i = 0; i < bx->nb; i++) {
333:     PetscInt tp;
334:     VecMin(bx->v[i], &tp, &lv);
335:     if (lv < *v) {
336:       *v = lv;
337:       lw = i;
338:       lp = tp;
339:     }
340:   }
341:   if (p && lw > -1) {
342:     PetscInt        st, en;
343:     const PetscInt *idxs;

345:     *p = -1;
346:     VecGetOwnershipRange(bx->v[lw], &st, &en);
347:     if (lp >= st && lp < en) {
348:       ISGetIndices(bx->is[lw], &idxs);
349:       *p = idxs[lp - st];
350:       ISRestoreIndices(bx->is[lw], &idxs);
351:     }
352:     MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x));
353:   }
354:   return 0;
355: }

357: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
358: {
359:   PetscInt  i, lp = -1, lw = -1;
360:   PetscReal lv;
361:   Vec_Nest *bx = (Vec_Nest *)x->data;

363:   *v = PETSC_MIN_REAL;
364:   for (i = 0; i < bx->nb; i++) {
365:     PetscInt tp;
366:     VecMax(bx->v[i], &tp, &lv);
367:     if (lv > *v) {
368:       *v = lv;
369:       lw = i;
370:       lp = tp;
371:     }
372:   }
373:   if (p && lw > -1) {
374:     PetscInt        st, en;
375:     const PetscInt *idxs;

377:     *p = -1;
378:     VecGetOwnershipRange(bx->v[lw], &st, &en);
379:     if (lp >= st && lp < en) {
380:       ISGetIndices(bx->is[lw], &idxs);
381:       *p = idxs[lp - st];
382:       ISRestoreIndices(bx->is[lw], &idxs);
383:     }
384:     MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x));
385:   }
386:   return 0;
387: }

389: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
390: {
391:   Vec_Nest *bx = (Vec_Nest *)x->data;
392:   PetscBool isascii;
393:   PetscInt  i;

395:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
396:   if (isascii) {
397:     PetscViewerASCIIPushTab(viewer);
398:     PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ",  structure: \n", bx->nb);
399:     for (i = 0; i < bx->nb; i++) {
400:       VecType  type;
401:       char     name[256] = "", prefix[256] = "";
402:       PetscInt NR;

404:       VecGetSize(bx->v[i], &NR);
405:       VecGetType(bx->v[i], &type);
406:       if (((PetscObject)bx->v[i])->name) PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name);
407:       if (((PetscObject)bx->v[i])->prefix) PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix);

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

411:       PetscViewerASCIIPushTab(viewer); /* push1 */
412:       VecView(bx->v[i], viewer);
413:       PetscViewerASCIIPopTab(viewer); /* pop1 */
414:     }
415:     PetscViewerASCIIPopTab(viewer); /* pop0 */
416:   }
417:   return 0;
418: }

420: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
421: {
422:   Vec_Nest *bx;
423:   PetscInt  size, i, nr;
424:   PetscBool isnest;

426:   PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest);
427:   if (!isnest) {
428:     /* Not nest */
429:     if (globalsize) VecGetSize(x, &size);
430:     else VecGetLocalSize(x, &size);
431:     *L = *L + size;
432:     return 0;
433:   }

435:   /* Otherwise we have a nest */
436:   bx = (Vec_Nest *)x->data;
437:   nr = bx->nb;

439:   /* now descend recursively */
440:   for (i = 0; i < nr; i++) VecSize_Nest_Recursive(bx->v[i], globalsize, L);
441:   return 0;
442: }

444: /* Returns the sum of the global size of all the constituent vectors in the nest */
445: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
446: {
447:   *N = x->map->N;
448:   return 0;
449: }

451: /* Returns the sum of the local size of all the constituent vectors in the nest */
452: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
453: {
454:   *n = x->map->n;
455:   return 0;
456: }

458: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
459: {
460:   Vec_Nest *bx = (Vec_Nest *)x->data;
461:   Vec_Nest *by = (Vec_Nest *)y->data;
462:   PetscInt  i, nr;
463:   PetscReal local_max, m;

465:   VecNestCheckCompatible2(x, 1, y, 2);
466:   nr = bx->nb;
467:   m  = 0.0;
468:   for (i = 0; i < nr; i++) {
469:     VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max);
470:     if (local_max > m) m = local_max;
471:   }
472:   *max = m;
473:   return 0;
474: }

476: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
477: {
478:   Vec_Nest *bx = (Vec_Nest *)X->data;
479:   PetscInt  i;

481:   *x = NULL;
482:   for (i = 0; i < bx->nb; i++) {
483:     PetscBool issame = PETSC_FALSE;
484:     ISEqual(is, bx->is[i], &issame);
485:     if (issame) {
486:       *x = bx->v[i];
487:       PetscObjectReference((PetscObject)(*x));
488:       break;
489:     }
490:   }
492:   return 0;
493: }

495: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
496: {
497:   PetscObjectStateIncrease((PetscObject)X);
498:   VecDestroy(x);
499:   return 0;
500: }

502: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
503: {
504:   Vec_Nest *bx = (Vec_Nest *)X->data;
505:   PetscInt  i, m, rstart, rend;

507:   VecGetOwnershipRange(X, &rstart, &rend);
508:   VecGetLocalSize(X, &m);
509:   PetscMalloc1(m, x);
510:   for (i = 0; i < bx->nb; i++) {
511:     Vec                subvec = bx->v[i];
512:     IS                 isy    = bx->is[i];
513:     PetscInt           j, sm;
514:     const PetscInt    *ixy;
515:     const PetscScalar *y;
516:     VecGetLocalSize(subvec, &sm);
517:     VecGetArrayRead(subvec, &y);
518:     ISGetIndices(isy, &ixy);
519:     for (j = 0; j < sm; j++) {
520:       PetscInt ix = ixy[j];
522:       (*x)[ix - rstart] = y[j];
523:     }
524:     ISRestoreIndices(isy, &ixy);
525:     VecRestoreArrayRead(subvec, &y);
526:   }
527:   return 0;
528: }

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

535:   VecGetOwnershipRange(X, &rstart, &rend);
536:   VecGetLocalSize(X, &m);
537:   for (i = 0; i < bx->nb; i++) {
538:     Vec             subvec = bx->v[i];
539:     IS              isy    = bx->is[i];
540:     PetscInt        j, sm;
541:     const PetscInt *ixy;
542:     PetscScalar    *y;
543:     VecGetLocalSize(subvec, &sm);
544:     VecGetArray(subvec, &y);
545:     ISGetIndices(isy, &ixy);
546:     for (j = 0; j < sm; j++) {
547:       PetscInt ix = ixy[j];
549:       y[j] = (*x)[ix - rstart];
550:     }
551:     ISRestoreIndices(isy, &ixy);
552:     VecRestoreArray(subvec, &y);
553:   }
554:   PetscFree(*x);
555:   PetscObjectStateIncrease((PetscObject)X);
556:   return 0;
557: }

559: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
560: {
561:   PetscFree(*x);
562:   return 0;
563: }

565: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
566: {
568:   return 0;
569: }

571: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
572: {
573:   Vec      *ww;
574:   IS       *wis;
575:   Vec_Nest *bv = (Vec_Nest *)v->data;
576:   PetscInt  i;

578:   PetscMalloc2(bv->nb, &ww, bv->nb, &wis);
579:   for (i = 0; i < bv->nb; i++) VecCreateLocalVector(bv->v[i], &ww[i]);
580:   for (i = 0; i < bv->nb; i++) {
581:     PetscInt off;

583:     VecGetOwnershipRange(v, &off, NULL);
584:     ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]);
585:     ISShift(wis[i], -off, wis[i]);
586:   }
587:   VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w);
588:   for (i = 0; i < bv->nb; i++) {
589:     VecDestroy(&ww[i]);
590:     ISDestroy(&wis[i]);
591:   }
592:   PetscFree2(ww, wis);
593:   return 0;
594: }

596: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
597: {
598:   Vec_Nest *bv = (Vec_Nest *)v->data;
599:   Vec_Nest *bw = (Vec_Nest *)w->data;
600:   PetscInt  i;

604:   for (i = 0; i < bv->nb; i++) VecGetLocalVector(bv->v[i], bw->v[i]);
605:   return 0;
606: }

608: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
609: {
610:   Vec_Nest *bv = (Vec_Nest *)v->data;
611:   Vec_Nest *bw = (Vec_Nest *)w->data;
612:   PetscInt  i;

616:   for (i = 0; i < bv->nb; i++) VecRestoreLocalVector(bv->v[i], bw->v[i]);
617:   PetscObjectStateIncrease((PetscObject)v);
618:   return 0;
619: }

621: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
622: {
623:   Vec_Nest *bv = (Vec_Nest *)v->data;
624:   Vec_Nest *bw = (Vec_Nest *)w->data;
625:   PetscInt  i;

629:   for (i = 0; i < bv->nb; i++) VecGetLocalVectorRead(bv->v[i], bw->v[i]);
630:   return 0;
631: }

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

641:   for (i = 0; i < bv->nb; i++) VecRestoreLocalVectorRead(bv->v[i], bw->v[i]);
642:   return 0;
643: }

645: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
646: {
647:   Vec_Nest *bv = (Vec_Nest *)v->data;
648:   PetscInt  i;

650:   for (i = 0; i < bv->nb; i++) VecSetRandom(bv->v[i], r);
651:   return 0;
652: }

654: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
655: {
656:   ops->duplicate               = VecDuplicate_Nest;
657:   ops->duplicatevecs           = VecDuplicateVecs_Default;
658:   ops->destroyvecs             = VecDestroyVecs_Default;
659:   ops->dot                     = VecDot_Nest;
660:   ops->mdot                    = VecMDot_Nest;
661:   ops->norm                    = VecNorm_Nest;
662:   ops->tdot                    = VecTDot_Nest;
663:   ops->mtdot                   = VecMTDot_Nest;
664:   ops->scale                   = VecScale_Nest;
665:   ops->copy                    = VecCopy_Nest;
666:   ops->set                     = VecSet_Nest;
667:   ops->swap                    = VecSwap_Nest;
668:   ops->axpy                    = VecAXPY_Nest;
669:   ops->axpby                   = VecAXPBY_Nest;
670:   ops->maxpy                   = VecMAXPY_Nest;
671:   ops->aypx                    = VecAYPX_Nest;
672:   ops->waxpy                   = VecWAXPY_Nest;
673:   ops->axpbypcz                = NULL;
674:   ops->pointwisemult           = VecPointwiseMult_Nest;
675:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
676:   ops->setvalues               = NULL;
677:   ops->assemblybegin           = VecAssemblyBegin_Nest;
678:   ops->assemblyend             = VecAssemblyEnd_Nest;
679:   ops->getarray                = VecGetArray_Nest;
680:   ops->getsize                 = VecGetSize_Nest;
681:   ops->getlocalsize            = VecGetLocalSize_Nest;
682:   ops->restorearray            = VecRestoreArray_Nest;
683:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
684:   ops->max                     = VecMax_Nest;
685:   ops->min                     = VecMin_Nest;
686:   ops->setrandom               = NULL;
687:   ops->setoption               = NULL;
688:   ops->setvaluesblocked        = NULL;
689:   ops->destroy                 = VecDestroy_Nest;
690:   ops->view                    = VecView_Nest;
691:   ops->placearray              = NULL;
692:   ops->replacearray            = NULL;
693:   ops->dot_local               = VecDot_Nest;
694:   ops->tdot_local              = VecTDot_Nest;
695:   ops->norm_local              = VecNorm_Nest;
696:   ops->mdot_local              = VecMDot_Nest;
697:   ops->mtdot_local             = VecMTDot_Nest;
698:   ops->load                    = NULL;
699:   ops->reciprocal              = VecReciprocal_Nest;
700:   ops->conjugate               = VecConjugate_Nest;
701:   ops->setlocaltoglobalmapping = NULL;
702:   ops->setvalueslocal          = NULL;
703:   ops->resetarray              = NULL;
704:   ops->setfromoptions          = NULL;
705:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
706:   ops->load                    = NULL;
707:   ops->pointwisemax            = NULL;
708:   ops->pointwisemaxabs         = NULL;
709:   ops->pointwisemin            = NULL;
710:   ops->getvalues               = NULL;
711:   ops->sqrt                    = NULL;
712:   ops->abs                     = NULL;
713:   ops->exp                     = NULL;
714:   ops->shift                   = NULL;
715:   ops->create                  = NULL;
716:   ops->stridegather            = NULL;
717:   ops->stridescatter           = NULL;
718:   ops->dotnorm2                = VecDotNorm2_Nest;
719:   ops->getsubvector            = VecGetSubVector_Nest;
720:   ops->restoresubvector        = VecRestoreSubVector_Nest;
721:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
722:   ops->concatenate             = VecConcatenate_Nest;
723:   ops->createlocalvector       = VecCreateLocalVector_Nest;
724:   ops->getlocalvector          = VecGetLocalVector_Nest;
725:   ops->getlocalvectorread      = VecGetLocalVectorRead_Nest;
726:   ops->restorelocalvector      = VecRestoreLocalVector_Nest;
727:   ops->restorelocalvectorread  = VecRestoreLocalVectorRead_Nest;
728:   ops->setrandom               = VecSetRandom_Nest;
729:   return 0;
730: }

732: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
733: {
734:   Vec_Nest *b = (Vec_Nest *)x->data;
735:   PetscInt  i;
736:   PetscInt  row;

738:   if (!m) return 0;
739:   for (i = 0; i < m; i++) {
740:     row = idxm[i];
742:     vec[i] = b->v[row];
743:   }
744:   return 0;
745: }

747: PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
748: {
749:   PetscObjectStateIncrease((PetscObject)X);
750:   VecNestGetSubVecs_Private(X, 1, &idxm, sx);
751:   return 0;
752: }

754: /*@
755:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

757:  Not collective

759:  Input Parameters:
760: +  X  - nest vector
761: -  idxm - index of the vector within the nest

763:  Output Parameter:
764: .  sx - vector at index idxm within the nest

766:  Notes:

768:  Level: developer

770: .seealso: `VecNestGetSize()`, `VecNestGetSubVecs()`
771: @*/
772: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
773: {
774:   PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
775:   return 0;
776: }

778: PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
779: {
780:   Vec_Nest *b = (Vec_Nest *)X->data;

782:   PetscObjectStateIncrease((PetscObject)X);
783:   if (N) *N = b->nb;
784:   if (sx) *sx = b->v;
785:   return 0;
786: }

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

791:  Not collective

793:  Input Parameter:
794: .  X  - nest vector

796:  Output Parameters:
797: +  N - number of nested vecs
798: -  sx - array of vectors

800:  Notes:
801:  The user should not free the array sx.

803:  Fortran Notes:
804:  The caller must allocate the array to hold the subvectors.

806:  Level: developer

808: .seealso: `VecNestGetSize()`, `VecNestGetSubVec()`
809: @*/
810: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec **sx)
811: {
812:   PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
813:   return 0;
814: }

816: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
817: {
818:   Vec_Nest *bx = (Vec_Nest *)X->data;
819:   PetscInt  i, offset = 0, n = 0, bs;
820:   IS        is;
821:   PetscBool issame = PETSC_FALSE;
822:   PetscInt  N      = 0;

824:   /* check if idxm < bx->nb */

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

831:   /* check if we need to update the IS for the block */
832:   offset = X->map->rstart;
833:   for (i = 0; i < idxm; i++) {
834:     n = 0;
835:     VecGetLocalSize(bx->v[i], &n);
836:     offset += n;
837:   }

839:   /* get the local size and block size */
840:   VecGetLocalSize(x, &n);
841:   VecGetBlockSize(x, &bs);

843:   /* create the new IS */
844:   ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is);
845:   ISSetBlockSize(is, bs);

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

850:   if (!issame) {
851:     /* The IS of given vector has a different layout compared to the existing block vector.
852:      Destroy the existing reference and update the IS. */
853:     ISDestroy(&bx->is[idxm]);
854:     ISDuplicate(is, &bx->is[idxm]);
855:     ISCopy(is, bx->is[idxm]);

857:     offset += n;
858:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
859:     for (i = idxm + 1; i < bx->nb; i++) {
860:       /* get the local size and block size */
861:       VecGetLocalSize(bx->v[i], &n);
862:       VecGetBlockSize(bx->v[i], &bs);

864:       /* destroy the old and create the new IS */
865:       ISDestroy(&bx->is[i]);
866:       ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]);
867:       ISSetBlockSize(bx->is[i], bs);

869:       offset += n;
870:     }

872:     n = 0;
873:     VecSize_Nest_Recursive(X, PETSC_TRUE, &N);
874:     VecSize_Nest_Recursive(X, PETSC_FALSE, &n);
875:     PetscLayoutSetSize(X->map, N);
876:     PetscLayoutSetLocalSize(X->map, n);
877:   }

879:   ISDestroy(&is);
880:   return 0;
881: }

883: PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
884: {
885:   PetscObjectStateIncrease((PetscObject)X);
886:   VecNestSetSubVec_Private(X, idxm, sx);
887:   return 0;
888: }

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

893:    Not collective

895:    Input Parameters:
896: +  X  - nest vector
897: .  idxm - index of the vector within the nest vector
898: -  sx - vector at index idxm within the nest vector

900:    Notes:
901:    The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.

903:    Level: developer

905: .seealso: `VecNestSetSubVecs()`, `VecNestGetSubVec()`
906: @*/
907: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
908: {
909:   PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
910:   return 0;
911: }

913: PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
914: {
915:   PetscInt i;

917:   PetscObjectStateIncrease((PetscObject)X);
918:   for (i = 0; i < N; i++) VecNestSetSubVec_Private(X, idxm[i], sx[i]);
919:   return 0;
920: }

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

925:    Not collective

927:    Input Parameters:
928: +  X  - nest vector
929: .  N - number of component vecs in sx
930: .  idxm - indices of component vecs that are to be replaced
931: -  sx - array of vectors

933:    Notes:
934:    The components in the vector array sx do not have to be of the same size as corresponding
935:    components in X. The user can also free the array "sx" after the call.

937:    Level: developer

939: .seealso: `VecNestGetSize()`, `VecNestGetSubVec()`
940: @*/
941: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
942: {
943:   PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
944:   return 0;
945: }

947: PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
948: {
949:   Vec_Nest *b = (Vec_Nest *)X->data;

951:   *N = b->nb;
952:   return 0;
953: }

955: /*@
956:  VecNestGetSize - Returns the size of the nest vector.

958:  Not collective

960:  Input Parameter:
961: .  X  - nest vector

963:  Output Parameter:
964: .  N - number of nested vecs

966:  Notes:

968:  Level: developer

970: .seealso: `VecNestGetSubVec()`, `VecNestGetSubVecs()`
971: @*/
972: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
973: {
976:   PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
977:   return 0;
978: }

980: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
981: {
982:   Vec_Nest *ctx = (Vec_Nest *)V->data;
983:   PetscInt  i;

985:   if (ctx->setup_called) return 0;

987:   ctx->nb = nb;

990:   /* Create space */
991:   PetscMalloc1(ctx->nb, &ctx->v);
992:   for (i = 0; i < ctx->nb; i++) {
993:     ctx->v[i] = x[i];
994:     PetscObjectReference((PetscObject)x[i]);
995:     /* Do not allocate memory for internal sub blocks */
996:   }

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

1000:   ctx->setup_called = PETSC_TRUE;
1001:   return 0;
1002: }

1004: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1005: {
1006:   Vec_Nest *ctx = (Vec_Nest *)V->data;
1007:   PetscInt  i, offset, m, n, M, N;

1009:   if (is) { /* Do some consistency checks and reference the is */
1010:     offset = V->map->rstart;
1011:     for (i = 0; i < ctx->nb; i++) {
1012:       ISGetSize(is[i], &M);
1013:       VecGetSize(ctx->v[i], &N);
1015:       ISGetLocalSize(is[i], &m);
1016:       VecGetLocalSize(ctx->v[i], &n);
1018:       PetscObjectReference((PetscObject)is[i]);
1019:       ctx->is[i] = is[i];
1020:       offset += n;
1021:     }
1022:   } else { /* Create a contiguous ISStride for each entry */
1023:     offset = V->map->rstart;
1024:     for (i = 0; i < ctx->nb; i++) {
1025:       PetscInt bs;
1026:       VecGetLocalSize(ctx->v[i], &n);
1027:       VecGetBlockSize(ctx->v[i], &bs);
1028:       ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]);
1029:       ISSetBlockSize(ctx->is[i], bs);
1030:       offset += n;
1031:     }
1032:   }
1033:   return 0;
1034: }

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

1039:    Collective

1041:    Input Parameters:
1042: +  comm - Communicator for the new `Vec`
1043: .  nb - number of nested blocks
1044: .  is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1045: -  x - array of nb sub-vectors

1047:    Output Parameter:
1048: .  Y - new vector

1050:    Level: advanced

1052: .seealso: `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`, `VECNEST`
1053: @*/
1054: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1055: {
1056:   Vec       V;
1057:   Vec_Nest *s;
1058:   PetscInt  n, N;

1060:   VecCreate(comm, &V);
1061:   PetscObjectChangeTypeName((PetscObject)V, VECNEST);

1063:   /* allocate and set pointer for implementation data */
1064:   PetscNew(&s);
1065:   V->data         = (void *)s;
1066:   s->setup_called = PETSC_FALSE;
1067:   s->nb           = -1;
1068:   s->v            = NULL;

1070:   VecSetUp_Nest_Private(V, nb, x);

1072:   n = N = 0;
1073:   VecSize_Nest_Recursive(V, PETSC_TRUE, &N);
1074:   VecSize_Nest_Recursive(V, PETSC_FALSE, &n);
1075:   PetscLayoutSetSize(V->map, N);
1076:   PetscLayoutSetLocalSize(V->map, n);
1077:   PetscLayoutSetBlockSize(V->map, 1);
1078:   PetscLayoutSetUp(V->map);

1080:   VecSetUp_NestIS_Private(V, nb, is);

1082:   VecNestSetOps_Private(V->ops);
1083:   V->petscnative = PETSC_FALSE;

1085:   /* expose block api's */
1086:   PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest);
1087:   PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest);
1088:   PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest);
1089:   PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest);
1090:   PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest);

1092:   *Y = V;
1093:   return 0;
1094: }

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

1099:   Level: intermediate

1101:   Notes:
1102:   This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1103:   It is usually used with MATNEST and DMComposite via DMSetVecType().

1105: .seealso: `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1106: M*/