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++) {
 23:     VecAssemblyEnd(vs->v[i]);
 24:   }
 25:   return 0;
 26: }

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

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

 49:   PetscFree(vs);
 50:   return 0;
 51: }

 53: /* supports nested blocks */
 54: static PetscErrorCode VecCopy_Nest(Vec x,Vec y)
 55: {
 56:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 57:   Vec_Nest       *by = (Vec_Nest*)y->data;
 58:   PetscInt       i;

 61:   VecNestCheckCompatible2(x,1,y,2);
 62:   for (i=0; i<bx->nb; i++) {
 63:     VecCopy(bx->v[i],by->v[i]);
 64:   }
 65:   return 0;
 66: }

 68: /* supports nested blocks */
 69: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
 70: {
 71:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 72:   Vec            Y;
 73:   Vec            *sub;
 74:   PetscInt       i;

 76:   PetscMalloc1(bx->nb,&sub);
 77:   for (i=0; i<bx->nb; i++) {
 78:     VecDuplicate(bx->v[i],&sub[i]);
 79:   }
 80:   VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
 81:   for (i=0; i<bx->nb; i++) {
 82:     VecDestroy(&sub[i]);
 83:   }
 84:   PetscFree(sub);
 85:   *y   = Y;
 86:   return 0;
 87: }

 89: /* supports nested blocks */
 90: static PetscErrorCode VecDot_Nest(Vec x,Vec y,PetscScalar *val)
 91: {
 92:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 93:   Vec_Nest       *by = (Vec_Nest*)y->data;
 94:   PetscInt       i,nr;
 95:   PetscScalar    x_dot_y,_val;

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

107: /* supports nested blocks */
108: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
109: {
110:   Vec_Nest       *bx = (Vec_Nest*)x->data;
111:   Vec_Nest       *by = (Vec_Nest*)y->data;
112:   PetscInt       i,nr;
113:   PetscScalar    x_dot_y,_val;

115:   nr   = bx->nb;
116:   _val = 0.0;
117:   for (i=0; i<nr; i++) {
118:     VecTDot(bx->v[i],by->v[i],&x_dot_y);
119:     _val = _val + x_dot_y;
120:   }
121:   *val = _val;
122:   return 0;
123: }

125: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
126: {
127:   Vec_Nest       *bx = (Vec_Nest*)x->data;
128:   Vec_Nest       *by = (Vec_Nest*)y->data;
129:   PetscInt       i,nr;
130:   PetscScalar    x_dot_y,_dp,_nm;
131:   PetscReal      norm2_y;

133:   nr  = bx->nb;
134:   _dp = 0.0;
135:   _nm = 0.0;
136:   for (i=0; i<nr; i++) {
137:     VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
138:     _dp += x_dot_y;
139:     _nm += norm2_y;
140:   }
141:   *dp = _dp;
142:   *nm = _nm;
143:   return 0;
144: }

146: static PetscErrorCode VecAXPY_Nest(Vec y,PetscScalar alpha,Vec x)
147: {
148:   Vec_Nest       *bx = (Vec_Nest*)x->data;
149:   Vec_Nest       *by = (Vec_Nest*)y->data;
150:   PetscInt       i,nr;

152:   nr = bx->nb;
153:   for (i=0; i<nr; i++) {
154:     VecAXPY(by->v[i],alpha,bx->v[i]);
155:   }
156:   return 0;
157: }

159: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
160: {
161:   Vec_Nest       *bx = (Vec_Nest*)x->data;
162:   Vec_Nest       *by = (Vec_Nest*)y->data;
163:   PetscInt       i,nr;

165:   nr = bx->nb;
166:   for (i=0; i<nr; i++) {
167:     VecAYPX(by->v[i],alpha,bx->v[i]);
168:   }
169:   return 0;
170: }

172: static PetscErrorCode VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
173: {
174:   Vec_Nest       *bx = (Vec_Nest*)x->data;
175:   Vec_Nest       *by = (Vec_Nest*)y->data;
176:   PetscInt       i,nr;

178:   nr = bx->nb;
179:   for (i=0; i<nr; i++) {
180:     VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
181:   }
182:   return 0;
183: }

185: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
186: {
187:   Vec_Nest       *bx = (Vec_Nest*)x->data;
188:   Vec_Nest       *by = (Vec_Nest*)y->data;
189:   Vec_Nest       *bz = (Vec_Nest*)z->data;
190:   PetscInt       i,nr;

192:   nr = bx->nb;
193:   for (i=0; i<nr; i++) {
194:     VecAXPBYPCZ(bz->v[i],alpha,beta,gamma,bx->v[i],by->v[i]);
195:   }
196:   return 0;
197: }

199: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
200: {
201:   Vec_Nest       *bx = (Vec_Nest*)x->data;
202:   PetscInt       i,nr;

204:   nr = bx->nb;
205:   for (i=0; i<nr; i++) {
206:     VecScale(bx->v[i],alpha);
207:   }
208:   return 0;
209: }

211: static PetscErrorCode VecPointwiseMult_Nest(Vec w,Vec x,Vec y)
212: {
213:   Vec_Nest       *bx = (Vec_Nest*)x->data;
214:   Vec_Nest       *by = (Vec_Nest*)y->data;
215:   Vec_Nest       *bw = (Vec_Nest*)w->data;
216:   PetscInt       i,nr;

218:   VecNestCheckCompatible3(w,1,x,2,y,3);
219:   nr = bx->nb;
220:   for (i=0; i<nr; i++) {
221:     VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
222:   }
223:   return 0;
224: }

226: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
227: {
228:   Vec_Nest       *bx = (Vec_Nest*)x->data;
229:   Vec_Nest       *by = (Vec_Nest*)y->data;
230:   Vec_Nest       *bw = (Vec_Nest*)w->data;
231:   PetscInt       i,nr;

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

235:   nr = bx->nb;
236:   for (i=0; i<nr; i++) {
237:     VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
238:   }
239:   return 0;
240: }

242: static PetscErrorCode VecReciprocal_Nest(Vec x)
243: {
244:   Vec_Nest       *bx = (Vec_Nest*)x->data;
245:   PetscInt       i,nr;

247:   nr = bx->nb;
248:   for (i=0; i<nr; i++) {
249:     VecReciprocal(bx->v[i]);
250:   }
251:   return 0;
252: }

254: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
255: {
256:   Vec_Nest       *bx = (Vec_Nest*)xin->data;
257:   PetscInt       i,nr;
258:   PetscReal      z_i;
259:   PetscReal      _z;

261:   nr = bx->nb;
262:   _z = 0.0;

264:   if (type == NORM_2) {
265:     PetscScalar dot;
266:     VecDot(xin,xin,&dot);
267:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
268:   } else if (type == NORM_1) {
269:     for (i=0; i<nr; i++) {
270:       VecNorm(bx->v[i],type,&z_i);
271:       _z = _z + z_i;
272:     }
273:   } else if (type == NORM_INFINITY) {
274:     for (i=0; i<nr; i++) {
275:       VecNorm(bx->v[i],type,&z_i);
276:       if (z_i > _z) _z = z_i;
277:     }
278:   }

280:   *z = _z;
281:   return 0;
282: }

284: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
285: {
286:   PetscInt       v;

288:   for (v=0; v<nv; v++) {
289:     /* Do axpy on each vector,v */
290:     VecAXPY(y,alpha[v],x[v]);
291:   }
292:   return 0;
293: }

295: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
296: {
297:   PetscInt       j;

299:   for (j=0; j<nv; j++) {
300:     VecDot(x,y[j],&val[j]);
301:   }
302:   return 0;
303: }

305: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
306: {
307:   PetscInt       j;

309:   for (j=0; j<nv; j++) {
310:     VecTDot(x,y[j],&val[j]);
311:   }
312:   return 0;
313: }

315: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
316: {
317:   Vec_Nest       *bx = (Vec_Nest*)x->data;
318:   PetscInt       i,nr;

320:   nr = bx->nb;
321:   for (i=0; i<nr; i++) {
322:     VecSet(bx->v[i],alpha);
323:   }
324:   return 0;
325: }

327: static PetscErrorCode VecConjugate_Nest(Vec x)
328: {
329:   Vec_Nest       *bx = (Vec_Nest*)x->data;
330:   PetscInt       j,nr;

332:   nr = bx->nb;
333:   for (j=0; j<nr; j++) {
334:     VecConjugate(bx->v[j]);
335:   }
336:   return 0;
337: }

339: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
340: {
341:   Vec_Nest       *bx = (Vec_Nest*)x->data;
342:   Vec_Nest       *by = (Vec_Nest*)y->data;
343:   PetscInt       i,nr;

345:   VecNestCheckCompatible2(x,1,y,2);
346:   nr = bx->nb;
347:   for (i=0; i<nr; i++) {
348:     VecSwap(bx->v[i],by->v[i]);
349:   }
350:   return 0;
351: }

353: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
354: {
355:   Vec_Nest       *bx = (Vec_Nest*)x->data;
356:   Vec_Nest       *by = (Vec_Nest*)y->data;
357:   Vec_Nest       *bw = (Vec_Nest*)w->data;
358:   PetscInt       i,nr;

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

362:   nr = bx->nb;
363:   for (i=0; i<nr; i++) {
364:     VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
365:   }
366:   return 0;
367: }

369: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
370: {
371:   Vec_Nest       *bx;
372:   PetscInt       i,nr;
373:   PetscBool      isnest;
374:   PetscInt       L;
375:   PetscInt       _entry_loc;
376:   PetscReal      _entry_val;

378:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
379:   if (!isnest) {
380:     /* Not nest */
381:     VecMax(x,&_entry_loc,&_entry_val);
382:     if (_entry_val > *max) {
383:       *max = _entry_val;
384:       if (p) *p = _entry_loc + *cnt;
385:     }
386:     VecGetSize(x,&L);
387:     *cnt = *cnt + L;
388:     return 0;
389:   }

391:   /* Otherwise we have a nest */
392:   bx = (Vec_Nest*)x->data;
393:   nr = bx->nb;

395:   /* now descend recursively */
396:   for (i=0; i<nr; i++) {
397:     VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
398:   }
399:   return 0;
400: }

402: /* supports nested blocks */
403: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
404: {
405:   PetscInt       cnt;

407:   cnt  = 0;
408:   if (p) *p = 0;
409:   *max = PETSC_MIN_REAL;
410:   VecMax_Nest_Recursive(x,&cnt,p,max);
411:   return 0;
412: }

414: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
415: {
416:   Vec_Nest       *bx = (Vec_Nest*)x->data;
417:   PetscInt       i,nr,L,_entry_loc;
418:   PetscBool      isnest;
419:   PetscReal      _entry_val;

421:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
422:   if (!isnest) {
423:     /* Not nest */
424:     VecMin(x,&_entry_loc,&_entry_val);
425:     if (_entry_val < *min) {
426:       *min = _entry_val;
427:       if (p) *p = _entry_loc + *cnt;
428:     }
429:     VecGetSize(x,&L);
430:     *cnt = *cnt + L;
431:     return 0;
432:   }

434:   /* Otherwise we have a nest */
435:   nr = bx->nb;

437:   /* now descend recursively */
438:   for (i=0; i<nr; i++) {
439:     VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
440:   }
441:   return 0;
442: }

444: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
445: {
446:   PetscInt       cnt;

448:   cnt  = 0;
449:   if (p) *p = 0;
450:   *min = PETSC_MAX_REAL;
451:   VecMin_Nest_Recursive(x,&cnt,p,min);
452:   return 0;
453: }

455: /* supports nested blocks */
456: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
457: {
458:   Vec_Nest       *bx = (Vec_Nest*)x->data;
459:   PetscBool      isascii;
460:   PetscInt       i;

462:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
463:   if (isascii) {
464:     PetscViewerASCIIPushTab(viewer);
465:     PetscViewerASCIIPrintf(viewer,"VecNest, rows=%" PetscInt_FMT ",  structure: \n",bx->nb);
466:     for (i=0; i<bx->nb; i++) {
467:       VecType  type;
468:       char     name[256] = "",prefix[256] = "";
469:       PetscInt NR;

471:       VecGetSize(bx->v[i],&NR);
472:       VecGetType(bx->v[i],&type);
473:       if (((PetscObject)bx->v[i])->name) PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);
474:       if (((PetscObject)bx->v[i])->prefix) PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);

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

478:       PetscViewerASCIIPushTab(viewer);             /* push1 */
479:       VecView(bx->v[i],viewer);
480:       PetscViewerASCIIPopTab(viewer);              /* pop1 */
481:     }
482:     PetscViewerASCIIPopTab(viewer);                /* pop0 */
483:   }
484:   return 0;
485: }

487: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
488: {
489:   Vec_Nest       *bx;
490:   PetscInt       size,i,nr;
491:   PetscBool      isnest;

493:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
494:   if (!isnest) {
495:     /* Not nest */
496:     if (globalsize) VecGetSize(x,&size);
497:     else VecGetLocalSize(x,&size);
498:     *L = *L + size;
499:     return 0;
500:   }

502:   /* Otherwise we have a nest */
503:   bx = (Vec_Nest*)x->data;
504:   nr = bx->nb;

506:   /* now descend recursively */
507:   for (i=0; i<nr; i++) {
508:     VecSize_Nest_Recursive(bx->v[i],globalsize,L);
509:   }
510:   return 0;
511: }

513: /* Returns the sum of the global size of all the consituent vectors in the nest */
514: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
515: {
516:   *N = x->map->N;
517:   return 0;
518: }

520: /* Returns the sum of the local size of all the consituent vectors in the nest */
521: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
522: {
523:   *n = x->map->n;
524:   return 0;
525: }

527: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
528: {
529:   Vec_Nest       *bx = (Vec_Nest*)x->data;
530:   Vec_Nest       *by = (Vec_Nest*)y->data;
531:   PetscInt       i,nr;
532:   PetscReal      local_max,m;

534:   VecNestCheckCompatible2(x,1,y,2);
535:   nr = bx->nb;
536:   m  = 0.0;
537:   for (i=0; i<nr; i++) {
538:     VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
539:     if (local_max > m) m = local_max;
540:   }
541:   *max = m;
542:   return 0;
543: }

545: static PetscErrorCode  VecGetSubVector_Nest(Vec X,IS is,Vec *x)
546: {
547:   Vec_Nest       *bx = (Vec_Nest*)X->data;
548:   PetscInt       i;

550:   *x = NULL;
551:   for (i=0; i<bx->nb; i++) {
552:     PetscBool issame = PETSC_FALSE;
553:     ISEqual(is,bx->is[i],&issame);
554:     if (issame) {
555:       *x   = bx->v[i];
556:       PetscObjectReference((PetscObject)(*x));
557:       break;
558:     }
559:   }
561:   return 0;
562: }

564: static PetscErrorCode  VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
565: {
566:   VecDestroy(x);
567:   return 0;
568: }

570: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
571: {
572:   Vec_Nest       *bx = (Vec_Nest*)X->data;
573:   PetscInt       i,m,rstart,rend;

575:   VecGetOwnershipRange(X,&rstart,&rend);
576:   VecGetLocalSize(X,&m);
577:   PetscMalloc1(m,x);
578:   for (i=0; i<bx->nb; i++) {
579:     Vec               subvec = bx->v[i];
580:     IS                isy    = bx->is[i];
581:     PetscInt          j,sm;
582:     const PetscInt    *ixy;
583:     const PetscScalar *y;
584:     VecGetLocalSize(subvec,&sm);
585:     VecGetArrayRead(subvec,&y);
586:     ISGetIndices(isy,&ixy);
587:     for (j=0; j<sm; j++) {
588:       PetscInt ix = ixy[j];
590:       (*x)[ix-rstart] = y[j];
591:     }
592:     ISRestoreIndices(isy,&ixy);
593:     VecRestoreArrayRead(subvec,&y);
594:   }
595:   return 0;
596: }

598: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
599: {
600:   Vec_Nest       *bx = (Vec_Nest*)X->data;
601:   PetscInt       i,m,rstart,rend;

603:   VecGetOwnershipRange(X,&rstart,&rend);
604:   VecGetLocalSize(X,&m);
605:   for (i=0; i<bx->nb; i++) {
606:     Vec            subvec = bx->v[i];
607:     IS             isy    = bx->is[i];
608:     PetscInt       j,sm;
609:     const PetscInt *ixy;
610:     PetscScalar    *y;
611:     VecGetLocalSize(subvec,&sm);
612:     VecGetArray(subvec,&y);
613:     ISGetIndices(isy,&ixy);
614:     for (j=0; j<sm; j++) {
615:       PetscInt ix = ixy[j];
617:       y[j] = (*x)[ix-rstart];
618:     }
619:     ISRestoreIndices(isy,&ixy);
620:     VecRestoreArray(subvec,&y);
621:   }
622:   PetscFree(*x);
623:   return 0;
624: }

626: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X,const PetscScalar **x)
627: {
628:   PetscFree(*x);
629:   return 0;
630: }

632: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
633: {
635:   return 0;
636: }

638: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
639: {
640:   ops->duplicate               = VecDuplicate_Nest;
641:   ops->duplicatevecs           = VecDuplicateVecs_Default;
642:   ops->destroyvecs             = VecDestroyVecs_Default;
643:   ops->dot                     = VecDot_Nest;
644:   ops->mdot                    = VecMDot_Nest;
645:   ops->norm                    = VecNorm_Nest;
646:   ops->tdot                    = VecTDot_Nest;
647:   ops->mtdot                   = VecMTDot_Nest;
648:   ops->scale                   = VecScale_Nest;
649:   ops->copy                    = VecCopy_Nest;
650:   ops->set                     = VecSet_Nest;
651:   ops->swap                    = VecSwap_Nest;
652:   ops->axpy                    = VecAXPY_Nest;
653:   ops->axpby                   = VecAXPBY_Nest;
654:   ops->maxpy                   = VecMAXPY_Nest;
655:   ops->aypx                    = VecAYPX_Nest;
656:   ops->waxpy                   = VecWAXPY_Nest;
657:   ops->axpbypcz                = NULL;
658:   ops->pointwisemult           = VecPointwiseMult_Nest;
659:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
660:   ops->setvalues               = NULL;
661:   ops->assemblybegin           = VecAssemblyBegin_Nest;
662:   ops->assemblyend             = VecAssemblyEnd_Nest;
663:   ops->getarray                = VecGetArray_Nest;
664:   ops->getsize                 = VecGetSize_Nest;
665:   ops->getlocalsize            = VecGetLocalSize_Nest;
666:   ops->restorearray            = VecRestoreArray_Nest;
667:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
668:   ops->max                     = VecMax_Nest;
669:   ops->min                     = VecMin_Nest;
670:   ops->setrandom               = NULL;
671:   ops->setoption               = NULL;
672:   ops->setvaluesblocked        = NULL;
673:   ops->destroy                 = VecDestroy_Nest;
674:   ops->view                    = VecView_Nest;
675:   ops->placearray              = NULL;
676:   ops->replacearray            = NULL;
677:   ops->dot_local               = VecDot_Nest;
678:   ops->tdot_local              = VecTDot_Nest;
679:   ops->norm_local              = VecNorm_Nest;
680:   ops->mdot_local              = VecMDot_Nest;
681:   ops->mtdot_local             = VecMTDot_Nest;
682:   ops->load                    = NULL;
683:   ops->reciprocal              = VecReciprocal_Nest;
684:   ops->conjugate               = VecConjugate_Nest;
685:   ops->setlocaltoglobalmapping = NULL;
686:   ops->setvalueslocal          = NULL;
687:   ops->resetarray              = NULL;
688:   ops->setfromoptions          = NULL;
689:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
690:   ops->load                    = NULL;
691:   ops->pointwisemax            = NULL;
692:   ops->pointwisemaxabs         = NULL;
693:   ops->pointwisemin            = NULL;
694:   ops->getvalues               = NULL;
695:   ops->sqrt                    = NULL;
696:   ops->abs                     = NULL;
697:   ops->exp                     = NULL;
698:   ops->shift                   = NULL;
699:   ops->create                  = NULL;
700:   ops->stridegather            = NULL;
701:   ops->stridescatter           = NULL;
702:   ops->dotnorm2                = VecDotNorm2_Nest;
703:   ops->getsubvector            = VecGetSubVector_Nest;
704:   ops->restoresubvector        = VecRestoreSubVector_Nest;
705:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
706:   ops->concatenate             = VecConcatenate_Nest;
707:   return 0;
708: }

710: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
711: {
712:   Vec_Nest *b = (Vec_Nest*)x->data;
713:   PetscInt i;
714:   PetscInt row;

716:   if (!m) return 0;
717:   for (i=0; i<m; i++) {
718:     row = idxm[i];
720:     vec[i] = b->v[row];
721:   }
722:   return 0;
723: }

725: PetscErrorCode  VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
726: {
727:   PetscObjectStateIncrease((PetscObject)X);
728:   VecNestGetSubVecs_Private(X,1,&idxm,sx);
729:   return 0;
730: }

732: /*@
733:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

735:  Not collective

737:  Input Parameters:
738: +  X  - nest vector
739: -  idxm - index of the vector within the nest

741:  Output Parameter:
742: .  sx - vector at index idxm within the nest

744:  Notes:

746:  Level: developer

748: .seealso: VecNestGetSize(), VecNestGetSubVecs()
749: @*/
750: PetscErrorCode  VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
751: {
752:   PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
753:   return 0;
754: }

756: PetscErrorCode  VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
757: {
758:   Vec_Nest *b = (Vec_Nest*)X->data;

760:   PetscObjectStateIncrease((PetscObject)X);
761:   if (N)  *N  = b->nb;
762:   if (sx) *sx = b->v;
763:   return 0;
764: }

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

769:  Not collective

771:  Input Parameter:
772: .  X  - nest vector

774:  Output Parameters:
775: +  N - number of nested vecs
776: -  sx - array of vectors

778:  Notes:
779:  The user should not free the array sx.

781:  Fortran Notes:
782:  The caller must allocate the array to hold the subvectors.

784:  Level: developer

786: .seealso: VecNestGetSize(), VecNestGetSubVec()
787: @*/
788: PetscErrorCode  VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
789: {
790:   PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
791:   return 0;
792: }

794: static PetscErrorCode  VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
795: {
796:   Vec_Nest       *bx = (Vec_Nest*)X->data;
797:   PetscInt       i,offset=0,n=0,bs;
798:   IS             is;
799:   PetscBool      issame = PETSC_FALSE;
800:   PetscInt       N=0;

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

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

809:   /* check if we need to update the IS for the block */
810:   offset = X->map->rstart;
811:   for (i=0; i<idxm; i++) {
812:     n=0;
813:     VecGetLocalSize(bx->v[i],&n);
814:     offset += n;
815:   }

817:   /* get the local size and block size */
818:   VecGetLocalSize(x,&n);
819:   VecGetBlockSize(x,&bs);

821:   /* create the new IS */
822:   ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
823:   ISSetBlockSize(is,bs);

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

828:   if (!issame) {
829:     /* The IS of given vector has a different layout compared to the existing block vector.
830:      Destroy the existing reference and update the IS. */
831:     ISDestroy(&bx->is[idxm]);
832:     ISDuplicate(is,&bx->is[idxm]);
833:     ISCopy(is,bx->is[idxm]);

835:     offset += n;
836:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
837:     for (i=idxm+1; i<bx->nb; i++) {
838:       /* get the local size and block size */
839:       VecGetLocalSize(bx->v[i],&n);
840:       VecGetBlockSize(bx->v[i],&bs);

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

847:       offset += n;
848:     }

850:     n=0;
851:     VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
852:     VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
853:     PetscLayoutSetSize(X->map,N);
854:     PetscLayoutSetLocalSize(X->map,n);
855:   }

857:   ISDestroy(&is);
858:   return 0;
859: }

861: PetscErrorCode  VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
862: {
863:   PetscObjectStateIncrease((PetscObject)X);
864:   VecNestSetSubVec_Private(X,idxm,sx);
865:   return 0;
866: }

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

871:    Not collective

873:    Input Parameters:
874: +  X  - nest vector
875: .  idxm - index of the vector within the nest vector
876: -  sx - vector at index idxm within the nest vector

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

881:    Level: developer

883: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
884: @*/
885: PetscErrorCode  VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
886: {
887:   PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
888:   return 0;
889: }

891: PetscErrorCode  VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
892: {
893:   PetscInt       i;

895:   PetscObjectStateIncrease((PetscObject)X);
896:   for (i=0; i<N; i++) {
897:     VecNestSetSubVec_Private(X,idxm[i],sx[i]);
898:   }
899:   return 0;
900: }

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

905:    Not collective

907:    Input Parameters:
908: +  X  - nest vector
909: .  N - number of component vecs in sx
910: .  idxm - indices of component vecs that are to be replaced
911: -  sx - array of vectors

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

917:    Level: developer

919: .seealso: VecNestGetSize(), VecNestGetSubVec()
920: @*/
921: PetscErrorCode  VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
922: {
923:   PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
924:   return 0;
925: }

927: PetscErrorCode  VecNestGetSize_Nest(Vec X,PetscInt *N)
928: {
929:   Vec_Nest *b = (Vec_Nest*)X->data;

931:   *N = b->nb;
932:   return 0;
933: }

935: /*@
936:  VecNestGetSize - Returns the size of the nest vector.

938:  Not collective

940:  Input Parameter:
941: .  X  - nest vector

943:  Output Parameter:
944: .  N - number of nested vecs

946:  Notes:

948:  Level: developer

950: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
951: @*/
952: PetscErrorCode  VecNestGetSize(Vec X,PetscInt *N)
953: {
956:   PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
957:   return 0;
958: }

960: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
961: {
962:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
963:   PetscInt       i;

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

967:   ctx->nb = nb;

970:   /* Create space */
971:   PetscMalloc1(ctx->nb,&ctx->v);
972:   for (i=0; i<ctx->nb; i++) {
973:     ctx->v[i] = x[i];
974:     PetscObjectReference((PetscObject)x[i]);
975:     /* Do not allocate memory for internal sub blocks */
976:   }

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

980:   ctx->setup_called = PETSC_TRUE;
981:   return 0;
982: }

984: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
985: {
986:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
987:   PetscInt       i,offset,m,n,M,N;

989:   if (is) {                     /* Do some consistency checks and reference the is */
990:     offset = V->map->rstart;
991:     for (i=0; i<ctx->nb; i++) {
992:       ISGetSize(is[i],&M);
993:       VecGetSize(ctx->v[i],&N);
995:       ISGetLocalSize(is[i],&m);
996:       VecGetLocalSize(ctx->v[i],&n);
998:       if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
999:         PetscInt  start;
1000:         PetscBool contiguous;
1001:         ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1004:       }
1005:       PetscObjectReference((PetscObject)is[i]);
1006:       ctx->is[i] = is[i];
1007:       offset += n;
1008:     }
1009:   } else {                      /* Create a contiguous ISStride for each entry */
1010:     offset = V->map->rstart;
1011:     for (i=0; i<ctx->nb; i++) {
1012:       PetscInt bs;
1013:       VecGetLocalSize(ctx->v[i],&n);
1014:       VecGetBlockSize(ctx->v[i],&bs);
1015:       ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1016:       ISSetBlockSize(ctx->is[i],bs);
1017:       offset += n;
1018:     }
1019:   }
1020:   return 0;
1021: }

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

1026:    Collective on Vec

1028:    Input Parameters:
1029: +  comm - Communicator for the new Vec
1030: .  nb - number of nested blocks
1031: .  is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1032: -  x - array of nb sub-vectors

1034:    Output Parameter:
1035: .  Y - new vector

1037:    Level: advanced

1039: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1040: @*/
1041: PetscErrorCode  VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1042: {
1043:   Vec            V;
1044:   Vec_Nest       *s;
1045:   PetscInt       n,N;

1047:   VecCreate(comm,&V);
1048:   PetscObjectChangeTypeName((PetscObject)V,VECNEST);

1050:   /* allocate and set pointer for implememtation data */
1051:   PetscNew(&s);
1052:   V->data          = (void*)s;
1053:   s->setup_called  = PETSC_FALSE;
1054:   s->nb            = -1;
1055:   s->v             = NULL;

1057:   VecSetUp_Nest_Private(V,nb,x);

1059:   n = N = 0;
1060:   VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1061:   VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1062:   PetscLayoutSetSize(V->map,N);
1063:   PetscLayoutSetLocalSize(V->map,n);
1064:   PetscLayoutSetBlockSize(V->map,1);
1065:   PetscLayoutSetUp(V->map);

1067:   VecSetUp_NestIS_Private(V,nb,is);

1069:   VecNestSetOps_Private(V->ops);
1070:   V->petscnative = PETSC_FALSE;

1072:   /* expose block api's */
1073:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1074:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1075:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1076:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1077:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);

1079:   *Y = V;
1080:   return 0;
1081: }

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

1086:   Level: intermediate

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

1092: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1093: M*/