Actual source code: vecnest.c

petsc-master 2017-02-20
Report Typos and Errors

  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;

 12:   for (i=0; i<vs->nb; i++) {
 13:     if (!vs->v[i]) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_SUP,"Nest  vector cannot contain NULL blocks");
 14:     VecAssemblyBegin(vs->v[i]);
 15:   }
 16:   return(0);
 17: }

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

 26:   for (i=0; i<vs->nb; i++) {
 27:     VecAssemblyEnd(vs->v[i]);
 28:   }
 29:   return(0);
 30: }

 32: static PetscErrorCode VecDestroy_Nest(Vec v)
 33: {
 34:   Vec_Nest       *vs = (Vec_Nest*)v->data;
 35:   PetscInt       i;

 39:   if (vs->v) {
 40:     for (i=0; i<vs->nb; i++) {
 41:       VecDestroy(&vs->v[i]);
 42:     }
 43:     PetscFree(vs->v);
 44:   }
 45:   for (i=0; i<vs->nb; i++) {
 46:     ISDestroy(&vs->is[i]);
 47:   }
 48:   PetscFree(vs->is);
 49:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 50:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 51:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 52:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 53:   PetscObjectComposeFunction((PetscObject)v,"",NULL);

 55:   PetscFree(vs);
 56:   return(0);
 57: }

 59: /* supports nested blocks */
 60: static PetscErrorCode VecCopy_Nest(Vec x,Vec y)
 61: {
 62:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 63:   Vec_Nest       *by = (Vec_Nest*)y->data;
 64:   PetscInt       i;

 68:   VecNestCheckCompatible2(x,1,y,2);
 69:   for (i=0; i<bx->nb; i++) {
 70:     VecCopy(bx->v[i],by->v[i]);
 71:   }
 72:   return(0);
 73: }

 75: /* supports nested blocks */
 76: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
 77: {
 78:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 79:   Vec            Y;
 80:   Vec            *sub;
 81:   PetscInt       i;

 85:   PetscMalloc1(bx->nb,&sub);
 86:   for (i=0; i<bx->nb; i++) {
 87:     VecDuplicate(bx->v[i],&sub[i]);
 88:   }
 89:   VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
 90:   for (i=0; i<bx->nb; i++) {
 91:     VecDestroy(&sub[i]);
 92:   }
 93:   PetscFree(sub);
 94:   *y   = Y;
 95:   return(0);
 96: }

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

108:   nr   = bx->nb;
109:   _val = 0.0;
110:   for (i=0; i<nr; i++) {
111:     VecDot(bx->v[i],by->v[i],&x_dot_y);
112:     _val = _val + x_dot_y;
113:   }
114:   *val = _val;
115:   return(0);
116: }

118: /* supports nested blocks */
119: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
120: {
121:   Vec_Nest       *bx = (Vec_Nest*)x->data;
122:   Vec_Nest       *by = (Vec_Nest*)y->data;
123:   PetscInt       i,nr;
124:   PetscScalar    x_dot_y,_val;

128:   nr   = bx->nb;
129:   _val = 0.0;
130:   for (i=0; i<nr; i++) {
131:     VecTDot(bx->v[i],by->v[i],&x_dot_y);
132:     _val = _val + x_dot_y;
133:   }
134:   *val = _val;
135:   return(0);
136: }

138: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
139: {
140:   Vec_Nest       *bx = (Vec_Nest*)x->data;
141:   Vec_Nest       *by = (Vec_Nest*)y->data;
142:   PetscInt       i,nr;
143:   PetscScalar    x_dot_y,_dp,_nm;
144:   PetscReal      norm2_y;

148:   nr  = bx->nb;
149:   _dp = 0.0;
150:   _nm = 0.0;
151:   for (i=0; i<nr; i++) {
152:     VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
153:     _dp += x_dot_y;
154:     _nm += norm2_y;
155:   }
156:   *dp = _dp;
157:   *nm = _nm;
158:   return(0);
159: }

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

169:   nr = bx->nb;
170:   for (i=0; i<nr; i++) {
171:     VecAXPY(by->v[i],alpha,bx->v[i]);
172:   }
173:   return(0);
174: }

176: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
177: {
178:   Vec_Nest       *bx = (Vec_Nest*)x->data;
179:   Vec_Nest       *by = (Vec_Nest*)y->data;
180:   PetscInt       i,nr;

184:   nr = bx->nb;
185:   for (i=0; i<nr; i++) {
186:     VecAYPX(by->v[i],alpha,bx->v[i]);
187:   }
188:   return(0);
189: }

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

199:   nr = bx->nb;
200:   for (i=0; i<nr; i++) {
201:     VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
202:   }
203:   return(0);
204: }

206: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
207: {
208:   Vec_Nest       *bx = (Vec_Nest*)x->data;
209:   PetscInt       i,nr;

213:   nr = bx->nb;
214:   for (i=0; i<nr; i++) {
215:     VecScale(bx->v[i],alpha);
216:   }
217:   return(0);
218: }

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

229:   VecNestCheckCompatible3(w,1,x,3,y,4);
230:   nr = bx->nb;
231:   for (i=0; i<nr; i++) {
232:     VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
233:   }
234:   return(0);
235: }

237: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
238: {
239:   Vec_Nest       *bx = (Vec_Nest*)x->data;
240:   Vec_Nest       *by = (Vec_Nest*)y->data;
241:   Vec_Nest       *bw = (Vec_Nest*)w->data;
242:   PetscInt       i,nr;

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

248:   nr = bx->nb;
249:   for (i=0; i<nr; i++) {
250:     VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
251:   }
252:   return(0);
253: }

255: static PetscErrorCode VecReciprocal_Nest(Vec x)
256: {
257:   Vec_Nest       *bx = (Vec_Nest*)x->data;
258:   PetscInt       i,nr;

262:   nr = bx->nb;
263:   for (i=0; i<nr; i++) {
264:     VecReciprocal(bx->v[i]);
265:   }
266:   return(0);
267: }

269: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
270: {
271:   Vec_Nest       *bx = (Vec_Nest*)xin->data;
272:   PetscInt       i,nr;
273:   PetscReal      z_i;
274:   PetscReal      _z;

278:   nr = bx->nb;
279:   _z = 0.0;

281:   if (type == NORM_2) {
282:     PetscScalar dot;
283:     VecDot(xin,xin,&dot);
284:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
285:   } else if (type == NORM_1) {
286:     for (i=0; i<nr; i++) {
287:       VecNorm(bx->v[i],type,&z_i);
288:       _z = _z + z_i;
289:     }
290:   } else if (type == NORM_INFINITY) {
291:     for (i=0; i<nr; i++) {
292:       VecNorm(bx->v[i],type,&z_i);
293:       if (z_i > _z) _z = z_i;
294:     }
295:   }

297:   *z = _z;
298:   return(0);
299: }

301: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
302: {
303:   PetscInt       v;

307:   for (v=0; v<nv; v++) {
308:     /* Do axpy on each vector,v */
309:     VecAXPY(y,alpha[v],x[v]);
310:   }
311:   return(0);
312: }

314: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
315: {
316:   PetscInt       j;

320:   for (j=0; j<nv; j++) {
321:     VecDot(x,y[j],&val[j]);
322:   }
323:   return(0);
324: }

326: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
327: {
328:   PetscInt       j;

332:   for (j=0; j<nv; j++) {
333:     VecTDot(x,y[j],&val[j]);
334:   }
335:   return(0);
336: }

338: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
339: {
340:   Vec_Nest       *bx = (Vec_Nest*)x->data;
341:   PetscInt       i,nr;

345:   nr = bx->nb;
346:   for (i=0; i<nr; i++) {
347:     VecSet(bx->v[i],alpha);
348:   }
349:   return(0);
350: }

352: static PetscErrorCode VecConjugate_Nest(Vec x)
353: {
354:   Vec_Nest       *bx = (Vec_Nest*)x->data;
355:   PetscInt       j,nr;

359:   nr = bx->nb;
360:   for (j=0; j<nr; j++) {
361:     VecConjugate(bx->v[j]);
362:   }
363:   return(0);
364: }

366: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
367: {
368:   Vec_Nest       *bx = (Vec_Nest*)x->data;
369:   Vec_Nest       *by = (Vec_Nest*)y->data;
370:   PetscInt       i,nr;

374:   VecNestCheckCompatible2(x,1,y,2);
375:   nr = bx->nb;
376:   for (i=0; i<nr; i++) {
377:     VecSwap(bx->v[i],by->v[i]);
378:   }
379:   return(0);
380: }

382: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
383: {
384:   Vec_Nest       *bx = (Vec_Nest*)x->data;
385:   Vec_Nest       *by = (Vec_Nest*)y->data;
386:   Vec_Nest       *bw = (Vec_Nest*)w->data;
387:   PetscInt       i,nr;

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

393:   nr = bx->nb;
394:   for (i=0; i<nr; i++) {
395:     VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
396:   }
397:   return(0);
398: }

400: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
401: {
402:   Vec_Nest       *bx;
403:   PetscInt       i,nr;
404:   PetscBool      isnest;
405:   PetscInt       L;
406:   PetscInt       _entry_loc;
407:   PetscReal      _entry_val;

411:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
412:   if (!isnest) {
413:     /* Not nest */
414:     VecMax(x,&_entry_loc,&_entry_val);
415:     if (_entry_val > *max) {
416:       *max = _entry_val;
417:       *p   = _entry_loc + *cnt;
418:     }
419:     VecGetSize(x,&L);
420:     *cnt = *cnt + L;
421:     return(0);
422:   }

424:   /* Otherwise we have a nest */
425:   bx = (Vec_Nest*)x->data;
426:   nr = bx->nb;

428:   /* now descend recursively */
429:   for (i=0; i<nr; i++) {
430:     VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
431:   }
432:   return(0);
433: }

435: /* supports nested blocks */
436: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
437: {
438:   PetscInt       cnt;

442:   cnt  = 0;
443:   *p   = 0;
444:   *max = PETSC_MIN_REAL;
445:   VecMax_Nest_Recursive(x,&cnt,p,max);
446:   return(0);
447: }

449: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
450: {
451:   Vec_Nest       *bx = (Vec_Nest*)x->data;
452:   PetscInt       i,nr,L,_entry_loc;
453:   PetscBool      isnest;
454:   PetscReal      _entry_val;

458:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
459:   if (!isnest) {
460:     /* Not nest */
461:     VecMin(x,&_entry_loc,&_entry_val);
462:     if (_entry_val < *min) {
463:       *min = _entry_val;
464:       *p   = _entry_loc + *cnt;
465:     }
466:     VecGetSize(x,&L);
467:     *cnt = *cnt + L;
468:     return(0);
469:   }

471:   /* Otherwise we have a nest */
472:   nr = bx->nb;

474:   /* now descend recursively */
475:   for (i=0; i<nr; i++) {
476:     VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
477:   }
478:   return(0);
479: }

481: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
482: {
483:   PetscInt       cnt;

487:   cnt  = 0;
488:   *p   = 0;
489:   *min = PETSC_MAX_REAL;
490:   VecMin_Nest_Recursive(x,&cnt,p,min);
491:   return(0);
492: }

494: /* supports nested blocks */
495: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
496: {
497:   Vec_Nest       *bx = (Vec_Nest*)x->data;
498:   PetscBool      isascii;
499:   PetscInt       i;

503:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
504:   if (isascii) {
505:     PetscViewerASCIIPushTab(viewer);
506:     PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D,  structure: \n",bx->nb);
507:     for (i=0; i<bx->nb; i++) {
508:       VecType  type;
509:       char     name[256] = "",prefix[256] = "";
510:       PetscInt NR;

512:       VecGetSize(bx->v[i],&NR);
513:       VecGetType(bx->v[i],&type);
514:       if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
515:       if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}

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

519:       PetscViewerASCIIPushTab(viewer);             /* push1 */
520:       VecView(bx->v[i],viewer);
521:       PetscViewerASCIIPopTab(viewer);              /* pop1 */
522:     }
523:     PetscViewerASCIIPopTab(viewer);                /* pop0 */
524:   }
525:   return(0);
526: }

528: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
529: {
530:   Vec_Nest       *bx;
531:   PetscInt       size,i,nr;
532:   PetscBool      isnest;

536:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
537:   if (!isnest) {
538:     /* Not nest */
539:     if (globalsize) { VecGetSize(x,&size); }
540:     else {            VecGetLocalSize(x,&size); }
541:     *L = *L + size;
542:     return(0);
543:   }

545:   /* Otherwise we have a nest */
546:   bx = (Vec_Nest*)x->data;
547:   nr = bx->nb;

549:   /* now descend recursively */
550:   for (i=0; i<nr; i++) {
551:     VecSize_Nest_Recursive(bx->v[i],globalsize,L);
552:   }
553:   return(0);
554: }

556: /* Returns the sum of the global size of all the consituent vectors in the nest */
557: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
558: {
560:   *N = x->map->N;
561:   return(0);
562: }

564: /* Returns the sum of the local size of all the consituent vectors in the nest */
565: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
566: {
568:   *n = x->map->n;
569:   return(0);
570: }

572: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
573: {
574:   Vec_Nest       *bx = (Vec_Nest*)x->data;
575:   Vec_Nest       *by = (Vec_Nest*)y->data;
576:   PetscInt       i,nr;
577:   PetscReal      local_max,m;

581:   VecNestCheckCompatible2(x,1,y,2);
582:   nr = bx->nb;
583:   m  = 0.0;
584:   for (i=0; i<nr; i++) {
585:     VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
586:     if (local_max > m) m = local_max;
587:   }
588:   *max = m;
589:   return(0);
590: }

592: static PetscErrorCode  VecGetSubVector_Nest(Vec X,IS is,Vec *x)
593: {
594:   Vec_Nest       *bx = (Vec_Nest*)X->data;
595:   PetscInt       i;

599:   *x = NULL;
600:   for (i=0; i<bx->nb; i++) {
601:     PetscBool issame = PETSC_FALSE;
602:     ISEqual(is,bx->is[i],&issame);
603:     if (issame) {
604:       *x   = bx->v[i];
605:       PetscObjectReference((PetscObject)(*x));
606:       break;
607:     }
608:   }
609:   if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
610:   return(0);
611: }

613: static PetscErrorCode  VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
614: {

618:   VecDestroy(x);
619:   return(0);
620: }

622: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
623: {
624:   Vec_Nest       *bx = (Vec_Nest*)X->data;
625:   PetscInt       i,m,rstart,rend;

629:   VecGetOwnershipRange(X,&rstart,&rend);
630:   VecGetLocalSize(X,&m);
631:   PetscMalloc1(m,x);
632:   for (i=0; i<bx->nb; i++) {
633:     Vec               subvec = bx->v[i];
634:     IS                isy    = bx->is[i];
635:     PetscInt          j,sm;
636:     const PetscInt    *ixy;
637:     const PetscScalar *y;
638:     VecGetLocalSize(subvec,&sm);
639:     VecGetArrayRead(subvec,&y);
640:     ISGetIndices(isy,&ixy);
641:     for (j=0; j<sm; j++) {
642:       PetscInt ix = ixy[j];
643:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
644:       (*x)[ix-rstart] = y[j];
645:     }
646:     ISRestoreIndices(isy,&ixy);
647:     VecRestoreArrayRead(subvec,&y);
648:   }
649:   return(0);
650: }

652: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
653: {
654:   Vec_Nest       *bx = (Vec_Nest*)X->data;
655:   PetscInt       i,m,rstart,rend;

659:   VecGetOwnershipRange(X,&rstart,&rend);
660:   VecGetLocalSize(X,&m);
661:   for (i=0; i<bx->nb; i++) {
662:     Vec            subvec = bx->v[i];
663:     IS             isy    = bx->is[i];
664:     PetscInt       j,sm;
665:     const PetscInt *ixy;
666:     PetscScalar    *y;
667:     VecGetLocalSize(subvec,&sm);
668:     VecGetArray(subvec,&y);
669:     ISGetIndices(isy,&ixy);
670:     for (j=0; j<sm; j++) {
671:       PetscInt ix = ixy[j];
672:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
673:       y[j] = (*x)[ix-rstart];
674:     }
675:     ISRestoreIndices(isy,&ixy);
676:     VecRestoreArray(subvec,&y);
677:   }
678:   PetscFree(*x);
679:   return(0);
680: }


683: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
684: {
686:   /* 0 */
687:   ops->duplicate               = VecDuplicate_Nest;
688:   ops->duplicatevecs           = VecDuplicateVecs_Default;
689:   ops->destroyvecs             = VecDestroyVecs_Default;
690:   ops->dot                     = VecDot_Nest;
691:   ops->mdot                    = VecMDot_Nest;

693:   /* 5 */
694:   ops->norm                    = VecNorm_Nest;
695:   ops->tdot                    = VecTDot_Nest;
696:   ops->mtdot                   = VecMTDot_Nest;
697:   ops->scale                   = VecScale_Nest;
698:   ops->copy                    = VecCopy_Nest;

700:   /* 10 */
701:   ops->set                     = VecSet_Nest;
702:   ops->swap                    = VecSwap_Nest;
703:   ops->axpy                    = VecAXPY_Nest;
704:   ops->axpby                   = VecAXPBY_Nest;
705:   ops->maxpy                   = VecMAXPY_Nest;

707:   /* 15 */
708:   ops->aypx                    = VecAYPX_Nest;
709:   ops->waxpy                   = VecWAXPY_Nest;
710:   ops->axpbypcz                = 0;
711:   ops->pointwisemult           = VecPointwiseMult_Nest;
712:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
713:   /* 20 */
714:   ops->setvalues               = 0;
715:   ops->assemblybegin           = VecAssemblyBegin_Nest;
716:   ops->assemblyend             = VecAssemblyEnd_Nest;
717:   ops->getarray                = VecGetArray_Nest;
718:   ops->getsize                 = VecGetSize_Nest;

720:   /* 25 */
721:   ops->getlocalsize            = VecGetLocalSize_Nest;
722:   ops->restorearray            = VecRestoreArray_Nest;
723:   ops->max                     = VecMax_Nest;
724:   ops->min                     = VecMin_Nest;
725:   ops->setrandom               = 0;

727:   /* 30 */
728:   ops->setoption               = 0;
729:   ops->setvaluesblocked        = 0;
730:   ops->destroy                 = VecDestroy_Nest;
731:   ops->view                    = VecView_Nest;
732:   ops->placearray              = 0;

734:   /* 35 */
735:   ops->replacearray            = 0;
736:   ops->dot_local               = VecDot_Nest;
737:   ops->tdot_local              = VecTDot_Nest;
738:   ops->norm_local              = VecNorm_Nest;
739:   ops->mdot_local              = VecMDot_Nest;

741:   /* 40 */
742:   ops->mtdot_local             = VecMTDot_Nest;
743:   ops->load                    = 0;
744:   ops->reciprocal              = VecReciprocal_Nest;
745:   ops->conjugate               = VecConjugate_Nest;
746:   ops->setlocaltoglobalmapping = 0;

748:   /* 45 */
749:   ops->setvalueslocal          = 0;
750:   ops->resetarray              = 0;
751:   ops->setfromoptions          = 0;
752:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
753:   ops->load                    = 0;

755:   /* 50 */
756:   ops->pointwisemax            = 0;
757:   ops->pointwisemaxabs         = 0;
758:   ops->pointwisemin            = 0;
759:   ops->getvalues               = 0;
760:   ops->sqrt                    = 0;

762:   /* 55 */
763:   ops->abs                     = 0;
764:   ops->exp                     = 0;
765:   ops->shift                   = 0;
766:   ops->create                  = 0;
767:   ops->stridegather            = 0;

769:   /* 60 */
770:   ops->stridescatter           = 0;
771:   ops->dotnorm2                = VecDotNorm2_Nest;
772:   ops->getsubvector            = VecGetSubVector_Nest;
773:   ops->restoresubvector        = VecRestoreSubVector_Nest;
774:   return(0);
775: }

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

784:   if (!m) return(0);
785:   for (i=0; i<m; i++) {
786:     row = idxm[i];
787:     if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
788:     vec[i] = b->v[row];
789:   }
790:   return(0);
791: }

793: PetscErrorCode  VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
794: {

798:   VecNestGetSubVecs_Private(X,1,&idxm,sx);
799:   return(0);
800: }

802: /*@
803:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

805:  Not collective

807:  Input Parameters:
808:  .  X  - nest vector
809:  .  idxm - index of the vector within the nest

811:  Output Parameter:
812:  .  sx - vector at index idxm within the nest

814:  Notes:

816:  Level: developer

818:  .seealso: VecNestGetSize(), VecNestGetSubVecs()
819: @*/
820: PetscErrorCode  VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
821: {

825:   PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
826:   return(0);
827: }

829: PetscErrorCode  VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
830: {
831:   Vec_Nest *b = (Vec_Nest*)X->data;

834:   if (N)  *N  = b->nb;
835:   if (sx) *sx = b->v;
836:   return(0);
837: }

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

842:  Not collective

844:  Input Parameters:
845: .  X  - nest vector

847:  Output Parameter:
848: +  N - number of nested vecs
849: -  sx - array of vectors

851:  Notes:
852:  The user should not free the array sx.

854:  Fortran Notes:
855:  The caller must allocate the array to hold the subvectors.

857:  Level: developer

859:  .seealso: VecNestGetSize(), VecNestGetSubVec()
860: @*/
861: PetscErrorCode  VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
862: {

866:   PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
867:   return(0);
868: }

870: static PetscErrorCode  VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
871: {
872:   Vec_Nest       *bx = (Vec_Nest*)X->data;
873:   PetscInt       i,offset=0,n=0,bs;
874:   IS             is;
876:   PetscBool      issame = PETSC_FALSE;
877:   PetscInt       N=0;

879:   /* check if idxm < bx->nb */
880:   if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);

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

887:   /* check if we need to update the IS for the block */
888:   offset = X->map->rstart;
889:   for (i=0; i<idxm; i++) {
890:     n=0;
891:     VecGetLocalSize(bx->v[i],&n);
892:     offset += n;
893:   }

895:   /* get the local size and block size */
896:   VecGetLocalSize(x,&n);
897:   VecGetBlockSize(x,&bs);

899:   /* create the new IS */
900:   ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
901:   ISSetBlockSize(is,bs);

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

906:   if (!issame) {
907:     /* The IS of given vector has a different layout compared to the existing block vector.
908:      Destroy the existing reference and update the IS. */
909:     ISDestroy(&bx->is[idxm]);
910:     ISDuplicate(is,&bx->is[idxm]);
911:     ISCopy(is,bx->is[idxm]);

913:     offset += n;
914:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
915:     for (i=idxm+1; i<bx->nb; i++) {
916:       /* get the local size and block size */
917:       VecGetLocalSize(bx->v[i],&n);
918:       VecGetBlockSize(bx->v[i],&bs);

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

925:       offset += n;
926:     }

928:     n=0;
929:     VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
930:     VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
931:     PetscLayoutSetSize(X->map,N);
932:     PetscLayoutSetLocalSize(X->map,n);
933:   }

935:   ISDestroy(&is);
936:   return(0);
937: }

939: PetscErrorCode  VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
940: {

944:   VecNestSetSubVec_Private(X,idxm,sx);
945:   return(0);
946: }

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

951:    Not collective

953:    Input Parameters:
954: +  X  - nest vector
955: .  idxm - index of the vector within the nest vector
956: -  sx - vector at index idxm within the nest vector

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

961:    Level: developer

963: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
964: @*/
965: PetscErrorCode  VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
966: {

970:   PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
971:   return(0);
972: }

974: PetscErrorCode  VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
975: {
976:   PetscInt       i;

980:   for (i=0; i<N; i++) {
981:     VecNestSetSubVec_Private(X,idxm[i],sx[i]);
982:   }
983:   return(0);
984: }

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

989:    Not collective

991:    Input Parameters:
992: +  X  - nest vector
993: .  N - number of component vecs in sx
994: .  idxm - indices of component vecs that are to be replaced
995: -  sx - array of vectors

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

1001:    Level: developer

1003: .seealso: VecNestGetSize(), VecNestGetSubVec()
1004: @*/
1005: PetscErrorCode  VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1006: {

1010:   PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1011:   return(0);
1012: }

1014: PetscErrorCode  VecNestGetSize_Nest(Vec X,PetscInt *N)
1015: {
1016:   Vec_Nest *b = (Vec_Nest*)X->data;

1019:   *N = b->nb;
1020:   return(0);
1021: }

1023: /*@
1024:  VecNestGetSize - Returns the size of the nest vector.

1026:  Not collective

1028:  Input Parameters:
1029:  .  X  - nest vector

1031:  Output Parameter:
1032:  .  N - number of nested vecs

1034:  Notes:

1036:  Level: developer

1038:  .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1039: @*/
1040: PetscErrorCode  VecNestGetSize(Vec X,PetscInt *N)
1041: {

1047:   PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1048:   return(0);
1049: }

1051: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1052: {
1053:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1054:   PetscInt       i;

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

1060:   ctx->nb = nb;
1061:   if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");

1063:   /* Create space */
1064:   PetscMalloc1(ctx->nb,&ctx->v);
1065:   for (i=0; i<ctx->nb; i++) {
1066:     ctx->v[i] = x[i];
1067:     PetscObjectReference((PetscObject)x[i]);
1068:     /* Do not allocate memory for internal sub blocks */
1069:   }

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

1073:   ctx->setup_called = PETSC_TRUE;
1074:   return(0);
1075: }

1077: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1078: {
1079:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1080:   PetscInt       i,offset,m,n,M,N;

1084:   if (is) {                     /* Do some consistency checks and reference the is */
1085:     offset = V->map->rstart;
1086:     for (i=0; i<ctx->nb; i++) {
1087:       ISGetSize(is[i],&M);
1088:       VecGetSize(ctx->v[i],&N);
1089:       if (M != N) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"In slot %D, IS of size %D is not compatible with Vec of size %D",i,M,N);
1090:       ISGetLocalSize(is[i],&m);
1091:       VecGetLocalSize(ctx->v[i],&n);
1092:       if (m != n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"In slot %D, IS of local size %D is not compatible with Vec of local size %D",i,m,n);
1093: #if defined(PETSC_USE_DEBUG)
1094:       {                         /* This test can be expensive */
1095:         PetscInt  start;
1096:         PetscBool contiguous;
1097:         ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1098:         if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1099:         if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1100:       }
1101: #endif
1102:       PetscObjectReference((PetscObject)is[i]);
1103:       ctx->is[i] = is[i];
1104:       offset += n;
1105:     }
1106:   } else {                      /* Create a contiguous ISStride for each entry */
1107:     offset = V->map->rstart;
1108:     for (i=0; i<ctx->nb; i++) {
1109:       PetscInt bs;
1110:       VecGetLocalSize(ctx->v[i],&n);
1111:       VecGetBlockSize(ctx->v[i],&bs);
1112:       ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1113:       ISSetBlockSize(ctx->is[i],bs);
1114:       offset += n;
1115:     }
1116:   }
1117:   return(0);
1118: }

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

1123:    Collective on Vec

1125:    Input Parameter:
1126: +  comm - Communicator for the new Vec
1127: .  nb - number of nested blocks
1128: .  is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1129: -  x - array of nb sub-vectors

1131:    Output Parameter:
1132: .  Y - new vector

1134:    Level: advanced

1136: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1137: @*/
1138: PetscErrorCode  VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1139: {
1140:   Vec            V;
1141:   Vec_Nest       *s;
1142:   PetscInt       n,N;

1146:   VecCreate(comm,&V);
1147:   PetscObjectChangeTypeName((PetscObject)V,VECNEST);

1149:   /* allocate and set pointer for implememtation data */
1150:   PetscNew(&s);
1151:   V->data          = (void*)s;
1152:   s->setup_called  = PETSC_FALSE;
1153:   s->nb            = -1;
1154:   s->v             = NULL;

1156:   VecSetUp_Nest_Private(V,nb,x);

1158:   n = N = 0;
1159:   VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1160:   VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1161:   PetscLayoutSetSize(V->map,N);
1162:   PetscLayoutSetLocalSize(V->map,n);
1163:   PetscLayoutSetBlockSize(V->map,1);
1164:   PetscLayoutSetUp(V->map);

1166:   VecSetUp_NestIS_Private(V,nb,is);

1168:   VecNestSetOps_Private(V->ops);
1169:   V->petscnative = PETSC_FALSE;


1172:   /* expose block api's */
1173:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1174:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1175:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1176:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1177:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);

1179:   *Y = V;
1180:   return(0);
1181: }

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

1186:   Level: intermediate

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

1192: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1193: M*/