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;

 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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

207: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
208: {
209:   Vec_Nest       *bx = (Vec_Nest*)x->data;
210:   Vec_Nest       *by = (Vec_Nest*)y->data;
211:   Vec_Nest       *bz = (Vec_Nest*)z->data;
212:   PetscInt       i,nr;

216:   nr = bx->nb;
217:   for (i=0; i<nr; i++) {
218:     VecAXPBYPCZ(bz->v[i],alpha,beta,gamma,bx->v[i],by->v[i]);
219:   }
220:   return(0);
221: }

223: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
224: {
225:   Vec_Nest       *bx = (Vec_Nest*)x->data;
226:   PetscInt       i,nr;

230:   nr = bx->nb;
231:   for (i=0; i<nr; i++) {
232:     VecScale(bx->v[i],alpha);
233:   }
234:   return(0);
235: }

237: static PetscErrorCode VecPointwiseMult_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);
247:   nr = bx->nb;
248:   for (i=0; i<nr; i++) {
249:     VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
250:   }
251:   return(0);
252: }

254: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
255: {
256:   Vec_Nest       *bx = (Vec_Nest*)x->data;
257:   Vec_Nest       *by = (Vec_Nest*)y->data;
258:   Vec_Nest       *bw = (Vec_Nest*)w->data;
259:   PetscInt       i,nr;

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

265:   nr = bx->nb;
266:   for (i=0; i<nr; i++) {
267:     VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
268:   }
269:   return(0);
270: }

272: static PetscErrorCode VecReciprocal_Nest(Vec x)
273: {
274:   Vec_Nest       *bx = (Vec_Nest*)x->data;
275:   PetscInt       i,nr;

279:   nr = bx->nb;
280:   for (i=0; i<nr; i++) {
281:     VecReciprocal(bx->v[i]);
282:   }
283:   return(0);
284: }

286: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
287: {
288:   Vec_Nest       *bx = (Vec_Nest*)xin->data;
289:   PetscInt       i,nr;
290:   PetscReal      z_i;
291:   PetscReal      _z;

295:   nr = bx->nb;
296:   _z = 0.0;

298:   if (type == NORM_2) {
299:     PetscScalar dot;
300:     VecDot(xin,xin,&dot);
301:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
302:   } else if (type == NORM_1) {
303:     for (i=0; i<nr; i++) {
304:       VecNorm(bx->v[i],type,&z_i);
305:       _z = _z + z_i;
306:     }
307:   } else if (type == NORM_INFINITY) {
308:     for (i=0; i<nr; i++) {
309:       VecNorm(bx->v[i],type,&z_i);
310:       if (z_i > _z) _z = z_i;
311:     }
312:   }

314:   *z = _z;
315:   return(0);
316: }

318: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
319: {
320:   PetscInt       v;

324:   for (v=0; v<nv; v++) {
325:     /* Do axpy on each vector,v */
326:     VecAXPY(y,alpha[v],x[v]);
327:   }
328:   return(0);
329: }

331: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
332: {
333:   PetscInt       j;

337:   for (j=0; j<nv; j++) {
338:     VecDot(x,y[j],&val[j]);
339:   }
340:   return(0);
341: }

343: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
344: {
345:   PetscInt       j;

349:   for (j=0; j<nv; j++) {
350:     VecTDot(x,y[j],&val[j]);
351:   }
352:   return(0);
353: }

355: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
356: {
357:   Vec_Nest       *bx = (Vec_Nest*)x->data;
358:   PetscInt       i,nr;

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

369: static PetscErrorCode VecConjugate_Nest(Vec x)
370: {
371:   Vec_Nest       *bx = (Vec_Nest*)x->data;
372:   PetscInt       j,nr;

376:   nr = bx->nb;
377:   for (j=0; j<nr; j++) {
378:     VecConjugate(bx->v[j]);
379:   }
380:   return(0);
381: }

383: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
384: {
385:   Vec_Nest       *bx = (Vec_Nest*)x->data;
386:   Vec_Nest       *by = (Vec_Nest*)y->data;
387:   PetscInt       i,nr;

391:   VecNestCheckCompatible2(x,1,y,2);
392:   nr = bx->nb;
393:   for (i=0; i<nr; i++) {
394:     VecSwap(bx->v[i],by->v[i]);
395:   }
396:   return(0);
397: }

399: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
400: {
401:   Vec_Nest       *bx = (Vec_Nest*)x->data;
402:   Vec_Nest       *by = (Vec_Nest*)y->data;
403:   Vec_Nest       *bw = (Vec_Nest*)w->data;
404:   PetscInt       i,nr;

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

410:   nr = bx->nb;
411:   for (i=0; i<nr; i++) {
412:     VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
413:   }
414:   return(0);
415: }

417: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
418: {
419:   Vec_Nest       *bx;
420:   PetscInt       i,nr;
421:   PetscBool      isnest;
422:   PetscInt       L;
423:   PetscInt       _entry_loc;
424:   PetscReal      _entry_val;

428:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
429:   if (!isnest) {
430:     /* Not nest */
431:     VecMax(x,&_entry_loc,&_entry_val);
432:     if (_entry_val > *max) {
433:       *max = _entry_val;
434:       if (p) *p = _entry_loc + *cnt;
435:     }
436:     VecGetSize(x,&L);
437:     *cnt = *cnt + L;
438:     return(0);
439:   }

441:   /* Otherwise we have a nest */
442:   bx = (Vec_Nest*)x->data;
443:   nr = bx->nb;

445:   /* now descend recursively */
446:   for (i=0; i<nr; i++) {
447:     VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
448:   }
449:   return(0);
450: }

452: /* supports nested blocks */
453: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
454: {
455:   PetscInt       cnt;

459:   cnt  = 0;
460:   if (p) *p = 0;
461:   *max = PETSC_MIN_REAL;
462:   VecMax_Nest_Recursive(x,&cnt,p,max);
463:   return(0);
464: }

466: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
467: {
468:   Vec_Nest       *bx = (Vec_Nest*)x->data;
469:   PetscInt       i,nr,L,_entry_loc;
470:   PetscBool      isnest;
471:   PetscReal      _entry_val;

475:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
476:   if (!isnest) {
477:     /* Not nest */
478:     VecMin(x,&_entry_loc,&_entry_val);
479:     if (_entry_val < *min) {
480:       *min = _entry_val;
481:       if (p) *p = _entry_loc + *cnt;
482:     }
483:     VecGetSize(x,&L);
484:     *cnt = *cnt + L;
485:     return(0);
486:   }

488:   /* Otherwise we have a nest */
489:   nr = bx->nb;

491:   /* now descend recursively */
492:   for (i=0; i<nr; i++) {
493:     VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
494:   }
495:   return(0);
496: }

498: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
499: {
500:   PetscInt       cnt;

504:   cnt  = 0;
505:   if (p) *p = 0;
506:   *min = PETSC_MAX_REAL;
507:   VecMin_Nest_Recursive(x,&cnt,p,min);
508:   return(0);
509: }

511: /* supports nested blocks */
512: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
513: {
514:   Vec_Nest       *bx = (Vec_Nest*)x->data;
515:   PetscBool      isascii;
516:   PetscInt       i;

520:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
521:   if (isascii) {
522:     PetscViewerASCIIPushTab(viewer);
523:     PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D,  structure: \n",bx->nb);
524:     for (i=0; i<bx->nb; i++) {
525:       VecType  type;
526:       char     name[256] = "",prefix[256] = "";
527:       PetscInt NR;

529:       VecGetSize(bx->v[i],&NR);
530:       VecGetType(bx->v[i],&type);
531:       if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
532:       if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}

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

536:       PetscViewerASCIIPushTab(viewer);             /* push1 */
537:       VecView(bx->v[i],viewer);
538:       PetscViewerASCIIPopTab(viewer);              /* pop1 */
539:     }
540:     PetscViewerASCIIPopTab(viewer);                /* pop0 */
541:   }
542:   return(0);
543: }

545: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
546: {
547:   Vec_Nest       *bx;
548:   PetscInt       size,i,nr;
549:   PetscBool      isnest;

553:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
554:   if (!isnest) {
555:     /* Not nest */
556:     if (globalsize) { VecGetSize(x,&size); }
557:     else {            VecGetLocalSize(x,&size); }
558:     *L = *L + size;
559:     return(0);
560:   }

562:   /* Otherwise we have a nest */
563:   bx = (Vec_Nest*)x->data;
564:   nr = bx->nb;

566:   /* now descend recursively */
567:   for (i=0; i<nr; i++) {
568:     VecSize_Nest_Recursive(bx->v[i],globalsize,L);
569:   }
570:   return(0);
571: }

573: /* Returns the sum of the global size of all the consituent vectors in the nest */
574: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
575: {
577:   *N = x->map->N;
578:   return(0);
579: }

581: /* Returns the sum of the local size of all the consituent vectors in the nest */
582: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
583: {
585:   *n = x->map->n;
586:   return(0);
587: }

589: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
590: {
591:   Vec_Nest       *bx = (Vec_Nest*)x->data;
592:   Vec_Nest       *by = (Vec_Nest*)y->data;
593:   PetscInt       i,nr;
594:   PetscReal      local_max,m;

598:   VecNestCheckCompatible2(x,1,y,2);
599:   nr = bx->nb;
600:   m  = 0.0;
601:   for (i=0; i<nr; i++) {
602:     VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
603:     if (local_max > m) m = local_max;
604:   }
605:   *max = m;
606:   return(0);
607: }

609: static PetscErrorCode  VecGetSubVector_Nest(Vec X,IS is,Vec *x)
610: {
611:   Vec_Nest       *bx = (Vec_Nest*)X->data;
612:   PetscInt       i;

616:   *x = NULL;
617:   for (i=0; i<bx->nb; i++) {
618:     PetscBool issame = PETSC_FALSE;
619:     ISEqual(is,bx->is[i],&issame);
620:     if (issame) {
621:       *x   = bx->v[i];
622:       PetscObjectReference((PetscObject)(*x));
623:       break;
624:     }
625:   }
626:   if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
627:   return(0);
628: }

630: static PetscErrorCode  VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
631: {

635:   VecDestroy(x);
636:   return(0);
637: }

639: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
640: {
641:   Vec_Nest       *bx = (Vec_Nest*)X->data;
642:   PetscInt       i,m,rstart,rend;

646:   VecGetOwnershipRange(X,&rstart,&rend);
647:   VecGetLocalSize(X,&m);
648:   PetscMalloc1(m,x);
649:   for (i=0; i<bx->nb; i++) {
650:     Vec               subvec = bx->v[i];
651:     IS                isy    = bx->is[i];
652:     PetscInt          j,sm;
653:     const PetscInt    *ixy;
654:     const PetscScalar *y;
655:     VecGetLocalSize(subvec,&sm);
656:     VecGetArrayRead(subvec,&y);
657:     ISGetIndices(isy,&ixy);
658:     for (j=0; j<sm; j++) {
659:       PetscInt ix = ixy[j];
660:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
661:       (*x)[ix-rstart] = y[j];
662:     }
663:     ISRestoreIndices(isy,&ixy);
664:     VecRestoreArrayRead(subvec,&y);
665:   }
666:   return(0);
667: }

669: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
670: {
671:   Vec_Nest       *bx = (Vec_Nest*)X->data;
672:   PetscInt       i,m,rstart,rend;

676:   VecGetOwnershipRange(X,&rstart,&rend);
677:   VecGetLocalSize(X,&m);
678:   for (i=0; i<bx->nb; i++) {
679:     Vec            subvec = bx->v[i];
680:     IS             isy    = bx->is[i];
681:     PetscInt       j,sm;
682:     const PetscInt *ixy;
683:     PetscScalar    *y;
684:     VecGetLocalSize(subvec,&sm);
685:     VecGetArray(subvec,&y);
686:     ISGetIndices(isy,&ixy);
687:     for (j=0; j<sm; j++) {
688:       PetscInt ix = ixy[j];
689:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
690:       y[j] = (*x)[ix-rstart];
691:     }
692:     ISRestoreIndices(isy,&ixy);
693:     VecRestoreArray(subvec,&y);
694:   }
695:   PetscFree(*x);
696:   return(0);
697: }

699: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X,const PetscScalar **x)
700: {

704:   PetscFree(*x);
705:   return(0);
706: }

708: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
709: {
711:   if (nx > 0) SETERRQ(PetscObjectComm((PetscObject)(*X)), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
712:   return(0);
713: }

715: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
716: {
718:   ops->duplicate               = VecDuplicate_Nest;
719:   ops->duplicatevecs           = VecDuplicateVecs_Default;
720:   ops->destroyvecs             = VecDestroyVecs_Default;
721:   ops->dot                     = VecDot_Nest;
722:   ops->mdot                    = VecMDot_Nest;
723:   ops->norm                    = VecNorm_Nest;
724:   ops->tdot                    = VecTDot_Nest;
725:   ops->mtdot                   = VecMTDot_Nest;
726:   ops->scale                   = VecScale_Nest;
727:   ops->copy                    = VecCopy_Nest;
728:   ops->set                     = VecSet_Nest;
729:   ops->swap                    = VecSwap_Nest;
730:   ops->axpy                    = VecAXPY_Nest;
731:   ops->axpby                   = VecAXPBY_Nest;
732:   ops->maxpy                   = VecMAXPY_Nest;
733:   ops->aypx                    = VecAYPX_Nest;
734:   ops->waxpy                   = VecWAXPY_Nest;
735:   ops->axpbypcz                = NULL;
736:   ops->pointwisemult           = VecPointwiseMult_Nest;
737:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
738:   ops->setvalues               = NULL;
739:   ops->assemblybegin           = VecAssemblyBegin_Nest;
740:   ops->assemblyend             = VecAssemblyEnd_Nest;
741:   ops->getarray                = VecGetArray_Nest;
742:   ops->getsize                 = VecGetSize_Nest;
743:   ops->getlocalsize            = VecGetLocalSize_Nest;
744:   ops->restorearray            = VecRestoreArray_Nest;
745:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
746:   ops->max                     = VecMax_Nest;
747:   ops->min                     = VecMin_Nest;
748:   ops->setrandom               = NULL;
749:   ops->setoption               = NULL;
750:   ops->setvaluesblocked        = NULL;
751:   ops->destroy                 = VecDestroy_Nest;
752:   ops->view                    = VecView_Nest;
753:   ops->placearray              = NULL;
754:   ops->replacearray            = NULL;
755:   ops->dot_local               = VecDot_Nest;
756:   ops->tdot_local              = VecTDot_Nest;
757:   ops->norm_local              = VecNorm_Nest;
758:   ops->mdot_local              = VecMDot_Nest;
759:   ops->mtdot_local             = VecMTDot_Nest;
760:   ops->load                    = NULL;
761:   ops->reciprocal              = VecReciprocal_Nest;
762:   ops->conjugate               = VecConjugate_Nest;
763:   ops->setlocaltoglobalmapping = NULL;
764:   ops->setvalueslocal          = NULL;
765:   ops->resetarray              = NULL;
766:   ops->setfromoptions          = NULL;
767:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
768:   ops->load                    = NULL;
769:   ops->pointwisemax            = NULL;
770:   ops->pointwisemaxabs         = NULL;
771:   ops->pointwisemin            = NULL;
772:   ops->getvalues               = NULL;
773:   ops->sqrt                    = NULL;
774:   ops->abs                     = NULL;
775:   ops->exp                     = NULL;
776:   ops->shift                   = NULL;
777:   ops->create                  = NULL;
778:   ops->stridegather            = NULL;
779:   ops->stridescatter           = NULL;
780:   ops->dotnorm2                = VecDotNorm2_Nest;
781:   ops->getsubvector            = VecGetSubVector_Nest;
782:   ops->restoresubvector        = VecRestoreSubVector_Nest;
783:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
784:   ops->concatenate             = VecConcatenate_Nest;
785:   return(0);
786: }

788: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
789: {
790:   Vec_Nest *b = (Vec_Nest*)x->data;
791:   PetscInt i;
792:   PetscInt row;

795:   if (!m) return(0);
796:   for (i=0; i<m; i++) {
797:     row = idxm[i];
798:     if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
799:     vec[i] = b->v[row];
800:   }
801:   return(0);
802: }

804: PetscErrorCode  VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
805: {

809:   VecNestGetSubVecs_Private(X,1,&idxm,sx);
810:   return(0);
811: }

813: /*@
814:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

816:  Not collective

818:  Input Parameters:
819: +  X  - nest vector
820: -  idxm - index of the vector within the nest

822:  Output Parameter:
823: .  sx - vector at index idxm within the nest

825:  Notes:

827:  Level: developer

829: .seealso: VecNestGetSize(), VecNestGetSubVecs()
830: @*/
831: PetscErrorCode  VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
832: {

836:   PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
837:   return(0);
838: }

840: PetscErrorCode  VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
841: {
842:   Vec_Nest *b = (Vec_Nest*)X->data;

845:   if (N)  *N  = b->nb;
846:   if (sx) *sx = b->v;
847:   return(0);
848: }

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

853:  Not collective

855:  Input Parameter:
856: .  X  - nest vector

858:  Output Parameters:
859: +  N - number of nested vecs
860: -  sx - array of vectors

862:  Notes:
863:  The user should not free the array sx.

865:  Fortran Notes:
866:  The caller must allocate the array to hold the subvectors.

868:  Level: developer

870: .seealso: VecNestGetSize(), VecNestGetSubVec()
871: @*/
872: PetscErrorCode  VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
873: {

877:   PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
878:   return(0);
879: }

881: static PetscErrorCode  VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
882: {
883:   Vec_Nest       *bx = (Vec_Nest*)X->data;
884:   PetscInt       i,offset=0,n=0,bs;
885:   IS             is;
887:   PetscBool      issame = PETSC_FALSE;
888:   PetscInt       N=0;

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

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

898:   /* check if we need to update the IS for the block */
899:   offset = X->map->rstart;
900:   for (i=0; i<idxm; i++) {
901:     n=0;
902:     VecGetLocalSize(bx->v[i],&n);
903:     offset += n;
904:   }

906:   /* get the local size and block size */
907:   VecGetLocalSize(x,&n);
908:   VecGetBlockSize(x,&bs);

910:   /* create the new IS */
911:   ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
912:   ISSetBlockSize(is,bs);

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

917:   if (!issame) {
918:     /* The IS of given vector has a different layout compared to the existing block vector.
919:      Destroy the existing reference and update the IS. */
920:     ISDestroy(&bx->is[idxm]);
921:     ISDuplicate(is,&bx->is[idxm]);
922:     ISCopy(is,bx->is[idxm]);

924:     offset += n;
925:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
926:     for (i=idxm+1; i<bx->nb; i++) {
927:       /* get the local size and block size */
928:       VecGetLocalSize(bx->v[i],&n);
929:       VecGetBlockSize(bx->v[i],&bs);

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

936:       offset += n;
937:     }

939:     n=0;
940:     VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
941:     VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
942:     PetscLayoutSetSize(X->map,N);
943:     PetscLayoutSetLocalSize(X->map,n);
944:   }

946:   ISDestroy(&is);
947:   return(0);
948: }

950: PetscErrorCode  VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
951: {

955:   VecNestSetSubVec_Private(X,idxm,sx);
956:   return(0);
957: }

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

962:    Not collective

964:    Input Parameters:
965: +  X  - nest vector
966: .  idxm - index of the vector within the nest vector
967: -  sx - vector at index idxm within the nest vector

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

972:    Level: developer

974: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
975: @*/
976: PetscErrorCode  VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
977: {

981:   PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
982:   return(0);
983: }

985: PetscErrorCode  VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
986: {
987:   PetscInt       i;

991:   for (i=0; i<N; i++) {
992:     VecNestSetSubVec_Private(X,idxm[i],sx[i]);
993:   }
994:   return(0);
995: }

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

1000:    Not collective

1002:    Input Parameters:
1003: +  X  - nest vector
1004: .  N - number of component vecs in sx
1005: .  idxm - indices of component vecs that are to be replaced
1006: -  sx - array of vectors

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

1012:    Level: developer

1014: .seealso: VecNestGetSize(), VecNestGetSubVec()
1015: @*/
1016: PetscErrorCode  VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1017: {

1021:   PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1022:   return(0);
1023: }

1025: PetscErrorCode  VecNestGetSize_Nest(Vec X,PetscInt *N)
1026: {
1027:   Vec_Nest *b = (Vec_Nest*)X->data;

1030:   *N = b->nb;
1031:   return(0);
1032: }

1034: /*@
1035:  VecNestGetSize - Returns the size of the nest vector.

1037:  Not collective

1039:  Input Parameter:
1040: .  X  - nest vector

1042:  Output Parameter:
1043: .  N - number of nested vecs

1045:  Notes:

1047:  Level: developer

1049: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1050: @*/
1051: PetscErrorCode  VecNestGetSize(Vec X,PetscInt *N)
1052: {

1058:   PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1059:   return(0);
1060: }

1062: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1063: {
1064:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1065:   PetscInt       i;

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

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

1074:   /* Create space */
1075:   PetscMalloc1(ctx->nb,&ctx->v);
1076:   for (i=0; i<ctx->nb; i++) {
1077:     ctx->v[i] = x[i];
1078:     PetscObjectReference((PetscObject)x[i]);
1079:     /* Do not allocate memory for internal sub blocks */
1080:   }

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

1084:   ctx->setup_called = PETSC_TRUE;
1085:   return(0);
1086: }

1088: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1089: {
1090:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1091:   PetscInt       i,offset,m,n,M,N;

1095:   if (is) {                     /* Do some consistency checks and reference the is */
1096:     offset = V->map->rstart;
1097:     for (i=0; i<ctx->nb; i++) {
1098:       ISGetSize(is[i],&M);
1099:       VecGetSize(ctx->v[i],&N);
1100:       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);
1101:       ISGetLocalSize(is[i],&m);
1102:       VecGetLocalSize(ctx->v[i],&n);
1103:       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);
1104:       if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
1105:         PetscInt  start;
1106:         PetscBool contiguous;
1107:         ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1108:         if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1109:         if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1110:       }
1111:       PetscObjectReference((PetscObject)is[i]);
1112:       ctx->is[i] = is[i];
1113:       offset += n;
1114:     }
1115:   } else {                      /* Create a contiguous ISStride for each entry */
1116:     offset = V->map->rstart;
1117:     for (i=0; i<ctx->nb; i++) {
1118:       PetscInt bs;
1119:       VecGetLocalSize(ctx->v[i],&n);
1120:       VecGetBlockSize(ctx->v[i],&bs);
1121:       ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1122:       ISSetBlockSize(ctx->is[i],bs);
1123:       offset += n;
1124:     }
1125:   }
1126:   return(0);
1127: }

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

1132:    Collective on Vec

1134:    Input Parameters:
1135: +  comm - Communicator for the new Vec
1136: .  nb - number of nested blocks
1137: .  is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1138: -  x - array of nb sub-vectors

1140:    Output Parameter:
1141: .  Y - new vector

1143:    Level: advanced

1145: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1146: @*/
1147: PetscErrorCode  VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1148: {
1149:   Vec            V;
1150:   Vec_Nest       *s;
1151:   PetscInt       n,N;

1155:   VecCreate(comm,&V);
1156:   PetscObjectChangeTypeName((PetscObject)V,VECNEST);

1158:   /* allocate and set pointer for implememtation data */
1159:   PetscNew(&s);
1160:   V->data          = (void*)s;
1161:   s->setup_called  = PETSC_FALSE;
1162:   s->nb            = -1;
1163:   s->v             = NULL;

1165:   VecSetUp_Nest_Private(V,nb,x);

1167:   n = N = 0;
1168:   VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1169:   VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1170:   PetscLayoutSetSize(V->map,N);
1171:   PetscLayoutSetLocalSize(V->map,n);
1172:   PetscLayoutSetBlockSize(V->map,1);
1173:   PetscLayoutSetUp(V->map);

1175:   VecSetUp_NestIS_Private(V,nb,is);

1177:   VecNestSetOps_Private(V->ops);
1178:   V->petscnative = PETSC_FALSE;

1180:   /* expose block api's */
1181:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1182:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1183:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1184:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1185:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);

1187:   *Y = V;
1188:   return(0);
1189: }

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

1194:   Level: intermediate

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

1200: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1201: M*/