Actual source code: vecnest.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: #include <../src/vec/vec/impls/nest/vecnestimpl.h>   /*I  "petscvec.h"   I*/

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

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

 23: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
 24: {
 25:   Vec_Nest       *vs = (Vec_Nest*)v->data;
 26:   PetscInt       i;

 30:   for (i=0; i<vs->nb; i++) {
 31:     VecAssemblyEnd(vs->v[i]);
 32:   }
 33:   return(0);
 34: }

 38: static PetscErrorCode VecDestroy_Nest(Vec v)
 39: {
 40:   Vec_Nest       *vs = (Vec_Nest*)v->data;
 41:   PetscInt       i;

 45:   if (vs->v) {
 46:     for (i=0; i<vs->nb; i++) {
 47:       VecDestroy(&vs->v[i]);
 48:     }
 49:     PetscFree(vs->v);
 50:   }
 51:   for (i=0; i<vs->nb; i++) {
 52:     ISDestroy(&vs->is[i]);
 53:   }
 54:   PetscFree(vs->is);
 55:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 56:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 57:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 58:   PetscObjectComposeFunction((PetscObject)v,"",NULL);
 59:   PetscObjectComposeFunction((PetscObject)v,"",NULL);

 61:   PetscFree(vs);
 62:   return(0);
 63: }

 65: /* supports nested blocks */
 68: static PetscErrorCode VecCopy_Nest(Vec x,Vec y)
 69: {
 70:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 71:   Vec_Nest       *by = (Vec_Nest*)y->data;
 72:   PetscInt       i;

 76:   VecNestCheckCompatible2(x,1,y,2);
 77:   for (i=0; i<bx->nb; i++) {
 78:     VecCopy(bx->v[i],by->v[i]);
 79:   }
 80:   return(0);
 81: }

 83: /* supports nested blocks */
 86: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
 87: {
 88:   Vec_Nest       *bx = (Vec_Nest*)x->data;
 89:   Vec            Y;
 90:   Vec            *sub;
 91:   PetscInt       i;

 95:   PetscMalloc(sizeof(Vec)*bx->nb,&sub);
 96:   for (i=0; i<bx->nb; i++) {
 97:     VecDuplicate(bx->v[i],&sub[i]);
 98:   }
 99:   VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
100:   for (i=0; i<bx->nb; i++) {
101:     VecDestroy(&sub[i]);
102:   }
103:   PetscFree(sub);
104:   *y   = Y;
105:   return(0);
106: }

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

120:   nr   = bx->nb;
121:   _val = 0.0;
122:   for (i=0; i<nr; i++) {
123:     VecDot(bx->v[i],by->v[i],&x_dot_y);
124:     _val = _val + x_dot_y;
125:   }
126:   *val = _val;
127:   return(0);
128: }

130: /* supports nested blocks */
133: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
134: {
135:   Vec_Nest       *bx = (Vec_Nest*)x->data;
136:   Vec_Nest       *by = (Vec_Nest*)y->data;
137:   PetscInt       i,nr;
138:   PetscScalar    x_dot_y,_val;

142:   nr   = bx->nb;
143:   _val = 0.0;
144:   for (i=0; i<nr; i++) {
145:     VecTDot(bx->v[i],by->v[i],&x_dot_y);
146:     _val = _val + x_dot_y;
147:   }
148:   *val = _val;
149:   return(0);
150: }

154: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
155: {
156:   Vec_Nest       *bx = (Vec_Nest*)x->data;
157:   Vec_Nest       *by = (Vec_Nest*)y->data;
158:   PetscInt       i,nr;
159:   PetscScalar    x_dot_y,_dp,_nm;
160:   PetscReal      norm2_y;

164:   nr  = bx->nb;
165:   _dp = 0.0;
166:   _nm = 0.0;
167:   for (i=0; i<nr; i++) {
168:     VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
169:     _dp += x_dot_y;
170:     _nm += norm2_y;
171:   }
172:   *dp = _dp;
173:   *nm = _nm;
174:   return(0);
175: }

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

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

196: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
197: {
198:   Vec_Nest       *bx = (Vec_Nest*)x->data;
199:   Vec_Nest       *by = (Vec_Nest*)y->data;
200:   PetscInt       i,nr;

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

213: static PetscErrorCode VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
214: {
215:   Vec_Nest       *bx = (Vec_Nest*)x->data;
216:   Vec_Nest       *by = (Vec_Nest*)y->data;
217:   PetscInt       i,nr;

221:   nr = bx->nb;
222:   for (i=0; i<nr; i++) {
223:     VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
224:   }
225:   return(0);
226: }

230: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
231: {
232:   Vec_Nest       *bx = (Vec_Nest*)x->data;
233:   PetscInt       i,nr;

237:   nr = bx->nb;
238:   for (i=0; i<nr; i++) {
239:     VecScale(bx->v[i],alpha);
240:   }
241:   return(0);
242: }

246: static PetscErrorCode VecPointwiseMult_Nest(Vec w,Vec x,Vec y)
247: {
248:   Vec_Nest       *bx = (Vec_Nest*)x->data;
249:   Vec_Nest       *by = (Vec_Nest*)y->data;
250:   Vec_Nest       *bw = (Vec_Nest*)w->data;
251:   PetscInt       i,nr;

255:   VecNestCheckCompatible3(w,1,x,3,y,4);
256:   nr = bx->nb;
257:   for (i=0; i<nr; i++) {
258:     VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
259:   }
260:   return(0);
261: }

265: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
266: {
267:   Vec_Nest       *bx = (Vec_Nest*)x->data;
268:   Vec_Nest       *by = (Vec_Nest*)y->data;
269:   Vec_Nest       *bw = (Vec_Nest*)w->data;
270:   PetscInt       i,nr;

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

276:   nr = bx->nb;
277:   for (i=0; i<nr; i++) {
278:     VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
279:   }
280:   return(0);
281: }

285: static PetscErrorCode VecReciprocal_Nest(Vec x)
286: {
287:   Vec_Nest       *bx = (Vec_Nest*)x->data;
288:   PetscInt       i,nr;

292:   nr = bx->nb;
293:   for (i=0; i<nr; i++) {
294:     VecReciprocal(bx->v[i]);
295:   }
296:   return(0);
297: }

301: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
302: {
303:   Vec_Nest       *bx = (Vec_Nest*)xin->data;
304:   PetscInt       i,nr;
305:   PetscReal      z_i;
306:   PetscReal      _z;

310:   nr = bx->nb;
311:   _z = 0.0;

313:   if (type == NORM_2) {
314:     PetscScalar dot;
315:     VecDot(xin,xin,&dot);
316:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
317:   } else if (type == NORM_1) {
318:     for (i=0; i<nr; i++) {
319:       VecNorm(bx->v[i],type,&z_i);
320:       _z = _z + z_i;
321:     }
322:   } else if (type == NORM_INFINITY) {
323:     for (i=0; i<nr; i++) {
324:       VecNorm(bx->v[i],type,&z_i);
325:       if (z_i > _z) _z = z_i;
326:     }
327:   }

329:   *z = _z;
330:   return(0);
331: }

335: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
336: {
337:   PetscInt       v;

341:   for (v=0; v<nv; v++) {
342:     /* Do axpy on each vector,v */
343:     VecAXPY(y,alpha[v],x[v]);
344:   }
345:   return(0);
346: }

350: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
351: {
352:   PetscInt       j;

356:   for (j=0; j<nv; j++) {
357:     VecDot(x,y[j],&val[j]);
358:   }
359:   return(0);
360: }

364: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
365: {
366:   PetscInt       j;

370:   for (j=0; j<nv; j++) {
371:     VecTDot(x,y[j],&val[j]);
372:   }
373:   return(0);
374: }

378: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
379: {
380:   Vec_Nest       *bx = (Vec_Nest*)x->data;
381:   PetscInt       i,nr;

385:   nr = bx->nb;
386:   for (i=0; i<nr; i++) {
387:     VecSet(bx->v[i],alpha);
388:   }
389:   return(0);
390: }

394: static PetscErrorCode VecConjugate_Nest(Vec x)
395: {
396:   Vec_Nest       *bx = (Vec_Nest*)x->data;
397:   PetscInt       j,nr;

401:   nr = bx->nb;
402:   for (j=0; j<nr; j++) {
403:     VecConjugate(bx->v[j]);
404:   }
405:   return(0);
406: }

410: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
411: {
412:   Vec_Nest       *bx = (Vec_Nest*)x->data;
413:   Vec_Nest       *by = (Vec_Nest*)y->data;
414:   PetscInt       i,nr;

418:   VecNestCheckCompatible2(x,1,y,2);
419:   nr = bx->nb;
420:   for (i=0; i<nr; i++) {
421:     VecSwap(bx->v[i],by->v[i]);
422:   }
423:   return(0);
424: }

428: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
429: {
430:   Vec_Nest       *bx = (Vec_Nest*)x->data;
431:   Vec_Nest       *by = (Vec_Nest*)y->data;
432:   Vec_Nest       *bw = (Vec_Nest*)w->data;
433:   PetscInt       i,nr;

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

439:   nr = bx->nb;
440:   for (i=0; i<nr; i++) {
441:     VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
442:   }
443:   return(0);
444: }

448: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
449: {
450:   Vec_Nest       *bx;
451:   PetscInt       i,nr;
452:   PetscBool      isnest;
453:   PetscInt       L;
454:   PetscInt       _entry_loc;
455:   PetscReal      _entry_val;

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

472:   /* Otherwise we have a nest */
473:   bx = (Vec_Nest*)x->data;
474:   nr = bx->nb;

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

483: /* supports nested blocks */
486: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
487: {
488:   PetscInt       cnt;

492:   cnt  = 0;
493:   *p   = 0;
494:   *max = 0.0;
495:   VecMax_Nest_Recursive(x,&cnt,p,max);
496:   return(0);
497: }

501: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
502: {
503:   Vec_Nest       *bx = (Vec_Nest*)x->data;
504:   PetscInt       i,nr,L,_entry_loc;
505:   PetscBool      isnest;
506:   PetscReal      _entry_val;

510:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
511:   if (!isnest) {
512:     /* Not nest */
513:     VecMin(x,&_entry_loc,&_entry_val);
514:     if (_entry_val < *min) {
515:       *min = _entry_val;
516:       *p   = _entry_loc + *cnt;
517:     }
518:     VecGetSize(x,&L);
519:     *cnt = *cnt + L;
520:     return(0);
521:   }

523:   /* Otherwise we have a nest */
524:   bx = (Vec_Nest*)x->data;
525:   nr = bx->nb;

527:   /* now descend recursively */
528:   for (i=0; i<nr; i++) {
529:     VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
530:   }
531:   return(0);
532: }

536: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
537: {
538:   PetscInt       cnt;

542:   cnt  = 0;
543:   *p   = 0;
544:   *min = 1.0e308;
545:   VecMin_Nest_Recursive(x,&cnt,p,min);
546:   return(0);
547: }

549: /* supports nested blocks */
552: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
553: {
554:   Vec_Nest       *bx = (Vec_Nest*)x->data;
555:   PetscBool      isascii;
556:   PetscInt       i;

560:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
561:   if (isascii) {
562:     PetscViewerASCIIPushTab(viewer);
563:     PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D,  structure: \n",bx->nb);
564:     for (i=0; i<bx->nb; i++) {
565:       VecType  type;
566:       char     name[256] = "",prefix[256] = "";
567:       PetscInt NR;

569:       VecGetSize(bx->v[i],&NR);
570:       VecGetType(bx->v[i],&type);
571:       if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
572:       if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}

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

576:       PetscViewerASCIIPushTab(viewer);             /* push1 */
577:       VecView(bx->v[i],viewer);
578:       PetscViewerASCIIPopTab(viewer);              /* pop1 */
579:     }
580:     PetscViewerASCIIPopTab(viewer);                /* pop0 */
581:   }
582:   return(0);
583: }

587: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
588: {
589:   Vec_Nest       *bx;
590:   PetscInt       size,i,nr;
591:   PetscBool      isnest;

595:   PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
596:   if (!isnest) {
597:     /* Not nest */
598:     if (globalsize) { VecGetSize(x,&size); }
599:     else {            VecGetLocalSize(x,&size); }
600:     *L = *L + size;
601:     return(0);
602:   }

604:   /* Otherwise we have a nest */
605:   bx = (Vec_Nest*)x->data;
606:   nr = bx->nb;

608:   /* now descend recursively */
609:   for (i=0; i<nr; i++) {
610:     VecSize_Nest_Recursive(bx->v[i],globalsize,L);
611:   }
612:   return(0);
613: }

615: /* Returns the sum of the global size of all the consituent vectors in the nest */
618: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
619: {
621:   *N = x->map->N;
622:   return(0);
623: }

625: /* Returns the sum of the local size of all the consituent vectors in the nest */
628: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
629: {
631:   *n = x->map->n;
632:   return(0);
633: }

637: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
638: {
639:   Vec_Nest       *bx = (Vec_Nest*)x->data;
640:   Vec_Nest       *by = (Vec_Nest*)y->data;
641:   PetscInt       i,nr;
642:   PetscReal      local_max,m;

646:   VecNestCheckCompatible2(x,1,y,2);
647:   nr = bx->nb;
648:   m  = 0.0;
649:   for (i=0; i<nr; i++) {
650:     VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
651:     if (local_max > m) m = local_max;
652:   }
653:   *max = m;
654:   return(0);
655: }

659: static PetscErrorCode  VecGetSubVector_Nest(Vec X,IS is,Vec *x)
660: {
661:   Vec_Nest       *bx = (Vec_Nest*)X->data;
662:   PetscInt       i;

666:   *x = NULL;
667:   for (i=0; i<bx->nb; i++) {
668:     PetscBool issame = PETSC_FALSE;
669:     ISEqual(is,bx->is[i],&issame);
670:     if (issame) {
671:       *x   = bx->v[i];
672:       PetscObjectReference((PetscObject)(*x));
673:       break;
674:     }
675:   }
676:   if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
677:   return(0);
678: }

682: static PetscErrorCode  VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
683: {

687:   VecDestroy(x);
688:   return(0);
689: }

693: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
694: {
695:   Vec_Nest       *bx = (Vec_Nest*)X->data;
696:   PetscInt       i,m,rstart,rend;

700:   VecGetOwnershipRange(X,&rstart,&rend);
701:   VecGetLocalSize(X,&m);
702:   PetscMalloc1(m,x);
703:   for (i=0; i<bx->nb; i++) {
704:     Vec               subvec = bx->v[i];
705:     IS                isy    = bx->is[i];
706:     PetscInt          j,sm;
707:     const PetscInt    *ixy;
708:     const PetscScalar *y;
709:     VecGetLocalSize(subvec,&sm);
710:     VecGetArrayRead(subvec,&y);
711:     ISGetIndices(isy,&ixy);
712:     for (j=0; j<sm; j++) {
713:       PetscInt ix = ixy[j];
714:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
715:       (*x)[ix-rstart] = y[j];
716:     }
717:     ISRestoreIndices(isy,&ixy);
718:     VecRestoreArrayRead(subvec,&y);
719:   }
720:   return(0);
721: }

725: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
726: {
727:   Vec_Nest       *bx = (Vec_Nest*)X->data;
728:   PetscInt       i,m,rstart,rend;

732:   VecGetOwnershipRange(X,&rstart,&rend);
733:   VecGetLocalSize(X,&m);
734:   for (i=0; i<bx->nb; i++) {
735:     Vec            subvec = bx->v[i];
736:     IS             isy    = bx->is[i];
737:     PetscInt       j,sm;
738:     const PetscInt *ixy;
739:     PetscScalar    *y;
740:     VecGetLocalSize(subvec,&sm);
741:     VecGetArray(subvec,&y);
742:     ISGetIndices(isy,&ixy);
743:     for (j=0; j<sm; j++) {
744:       PetscInt ix = ixy[j];
745:       if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
746:       y[j] = (*x)[ix-rstart];
747:     }
748:     ISRestoreIndices(isy,&ixy);
749:     VecRestoreArray(subvec,&y);
750:   }
751:   PetscFree(*x);
752:   return(0);
753: }


758: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
759: {
761:   /* 0 */
762:   ops->duplicate               = VecDuplicate_Nest;
763:   ops->duplicatevecs           = VecDuplicateVecs_Default;
764:   ops->destroyvecs             = VecDestroyVecs_Default;
765:   ops->dot                     = VecDot_Nest;
766:   ops->mdot                    = VecMDot_Nest;

768:   /* 5 */
769:   ops->norm                    = VecNorm_Nest;
770:   ops->tdot                    = VecTDot_Nest;
771:   ops->mtdot                   = VecMTDot_Nest;
772:   ops->scale                   = VecScale_Nest;
773:   ops->copy                    = VecCopy_Nest;

775:   /* 10 */
776:   ops->set                     = VecSet_Nest;
777:   ops->swap                    = VecSwap_Nest;
778:   ops->axpy                    = VecAXPY_Nest;
779:   ops->axpby                   = VecAXPBY_Nest;
780:   ops->maxpy                   = VecMAXPY_Nest;

782:   /* 15 */
783:   ops->aypx                    = VecAYPX_Nest;
784:   ops->waxpy                   = VecWAXPY_Nest;
785:   ops->axpbypcz                = 0;
786:   ops->pointwisemult           = VecPointwiseMult_Nest;
787:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
788:   /* 20 */
789:   ops->setvalues               = 0;
790:   ops->assemblybegin           = VecAssemblyBegin_Nest;
791:   ops->assemblyend             = VecAssemblyEnd_Nest;
792:   ops->getarray                = VecGetArray_Nest;
793:   ops->getsize                 = VecGetSize_Nest;

795:   /* 25 */
796:   ops->getlocalsize            = VecGetLocalSize_Nest;
797:   ops->restorearray            = VecRestoreArray_Nest;
798:   ops->max                     = VecMax_Nest;
799:   ops->min                     = VecMin_Nest;
800:   ops->setrandom               = 0;

802:   /* 30 */
803:   ops->setoption               = 0;
804:   ops->setvaluesblocked        = 0;
805:   ops->destroy                 = VecDestroy_Nest;
806:   ops->view                    = VecView_Nest;
807:   ops->placearray              = 0;

809:   /* 35 */
810:   ops->replacearray            = 0;
811:   ops->dot_local               = VecDot_Nest;
812:   ops->tdot_local              = VecTDot_Nest;
813:   ops->norm_local              = VecNorm_Nest;
814:   ops->mdot_local              = VecMDot_Nest;

816:   /* 40 */
817:   ops->mtdot_local             = VecMTDot_Nest;
818:   ops->load                    = 0;
819:   ops->reciprocal              = VecReciprocal_Nest;
820:   ops->conjugate               = VecConjugate_Nest;
821:   ops->setlocaltoglobalmapping = 0;

823:   /* 45 */
824:   ops->setvalueslocal          = 0;
825:   ops->resetarray              = 0;
826:   ops->setfromoptions          = 0;
827:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
828:   ops->load                    = 0;

830:   /* 50 */
831:   ops->pointwisemax            = 0;
832:   ops->pointwisemaxabs         = 0;
833:   ops->pointwisemin            = 0;
834:   ops->getvalues               = 0;
835:   ops->sqrt                    = 0;

837:   /* 55 */
838:   ops->abs                     = 0;
839:   ops->exp                     = 0;
840:   ops->shift                   = 0;
841:   ops->create                  = 0;
842:   ops->stridegather            = 0;

844:   /* 60 */
845:   ops->stridescatter           = 0;
846:   ops->dotnorm2                = VecDotNorm2_Nest;
847:   ops->getsubvector            = VecGetSubVector_Nest;
848:   ops->restoresubvector        = VecRestoreSubVector_Nest;
849:   return(0);
850: }

854: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
855: {
856:   Vec_Nest *b = (Vec_Nest*)x->data;
857:   PetscInt i;
858:   PetscInt row;

861:   if (!m) return(0);
862:   for (i=0; i<m; i++) {
863:     row = idxm[i];
864:     if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
865:     vec[i] = b->v[row];
866:   }
867:   return(0);
868: }

872: PetscErrorCode  VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
873: {

877:   VecNestGetSubVecs_Private(X,1,&idxm,sx);
878:   return(0);
879: }

883: /*@
884:  VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

886:  Not collective

888:  Input Parameters:
889:  .  X  - nest vector
890:  .  idxm - index of the vector within the nest

892:  Output Parameter:
893:  .  sx - vector at index idxm within the nest

895:  Notes:

897:  Level: developer

899:  .seealso: VecNestGetSize(), VecNestGetSubVecs()
900: @*/
901: PetscErrorCode  VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
902: {

906:   PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
907:   return(0);
908: }

912: PetscErrorCode  VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
913: {
914:   Vec_Nest *b = (Vec_Nest*)X->data;

917:   if (N)  *N  = b->nb;
918:   if (sx) *sx = b->v;
919:   return(0);
920: }

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

927:  Not collective

929:  Input Parameters:
930: .  X  - nest vector

932:  Output Parameter:
933: +  N - number of nested vecs
934: -  sx - array of vectors

936:  Notes:
937:  The user should not free the array sx.

939:  Fortran Notes:
940:  The caller must allocate the array to hold the subvectors.

942:  Level: developer

944:  .seealso: VecNestGetSize(), VecNestGetSubVec()
945: @*/
946: PetscErrorCode  VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
947: {

951:   PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
952:   return(0);
953: }

957: static PetscErrorCode  VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
958: {
959:   Vec_Nest       *bx = (Vec_Nest*)X->data;
960:   PetscInt       i,offset=0,n=0,bs;
961:   IS             is;
963:   PetscBool      issame = PETSC_FALSE;
964:   PetscInt       N=0;

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

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

974:   /* check if we need to update the IS for the block */
975:   offset = X->map->rstart;
976:   for (i=0; i<idxm; i++) {
977:     n=0;
978:     VecGetLocalSize(bx->v[i],&n);
979:     offset += n;
980:   }

982:   /* get the local size and block size */
983:   VecGetLocalSize(x,&n);
984:   VecGetBlockSize(x,&bs);

986:   /* create the new IS */
987:   ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
988:   ISSetBlockSize(is,bs);

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

993:   if (!issame) {
994:     /* The IS of given vector has a different layout compared to the existing block vector.
995:      Destroy the existing reference and update the IS. */
996:     ISDestroy(&bx->is[idxm]);
997:     ISDuplicate(is,&bx->is[idxm]);
998:     ISCopy(is,bx->is[idxm]);

1000:     offset += n;
1001:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
1002:     for (i=idxm+1; i<bx->nb; i++) {
1003:       /* get the local size and block size */
1004:       VecGetLocalSize(bx->v[i],&n);
1005:       VecGetBlockSize(bx->v[i],&bs);

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

1012:       offset += n;
1013:     }

1015:     n=0;
1016:     VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
1017:     VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
1018:     PetscLayoutSetSize(X->map,N);
1019:     PetscLayoutSetLocalSize(X->map,n);
1020:   }

1022:   ISDestroy(&is);
1023:   return(0);
1024: }

1028: PetscErrorCode  VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
1029: {

1033:   VecNestSetSubVec_Private(X,idxm,sx);
1034:   return(0);
1035: }

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

1042:    Not collective

1044:    Input Parameters:
1045: +  X  - nest vector
1046: .  idxm - index of the vector within the nest vector
1047: -  sx - vector at index idxm within the nest vector

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

1052:    Level: developer

1054: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
1055: @*/
1056: PetscErrorCode  VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
1057: {

1061:   PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
1062:   return(0);
1063: }

1067: PetscErrorCode  VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1068: {
1069:   PetscInt       i;

1073:   for (i=0; i<N; i++) {
1074:     VecNestSetSubVec_Private(X,idxm[i],sx[i]);
1075:   }
1076:   return(0);
1077: }

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

1084:    Not collective

1086:    Input Parameters:
1087: +  X  - nest vector
1088: .  N - number of component vecs in sx
1089: .  idxm - indices of component vecs that are to be replaced
1090: -  sx - array of vectors

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

1096:    Level: developer

1098: .seealso: VecNestGetSize(), VecNestGetSubVec()
1099: @*/
1100: PetscErrorCode  VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1101: {

1105:   PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1106:   return(0);
1107: }

1111: PetscErrorCode  VecNestGetSize_Nest(Vec X,PetscInt *N)
1112: {
1113:   Vec_Nest *b = (Vec_Nest*)X->data;

1116:   *N = b->nb;
1117:   return(0);
1118: }

1122: /*@
1123:  VecNestGetSize - Returns the size of the nest vector.

1125:  Not collective

1127:  Input Parameters:
1128:  .  X  - nest vector

1130:  Output Parameter:
1131:  .  N - number of nested vecs

1133:  Notes:

1135:  Level: developer

1137:  .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1138: @*/
1139: PetscErrorCode  VecNestGetSize(Vec X,PetscInt *N)
1140: {

1146:   PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1147:   return(0);
1148: }

1152: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1153: {
1154:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1155:   PetscInt       i;

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

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

1164:   /* Create space */
1165:   PetscMalloc1(ctx->nb,&ctx->v);
1166:   for (i=0; i<ctx->nb; i++) {
1167:     ctx->v[i] = x[i];
1168:     PetscObjectReference((PetscObject)x[i]);
1169:     /* Do not allocate memory for internal sub blocks */
1170:   }

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

1174:   ctx->setup_called = PETSC_TRUE;
1175:   return(0);
1176: }

1180: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1181: {
1182:   Vec_Nest       *ctx = (Vec_Nest*)V->data;
1183:   PetscInt       i,offset,m,n,M,N;

1187:   if (is) {                     /* Do some consistency checks and reference the is */
1188:     offset = V->map->rstart;
1189:     for (i=0; i<ctx->nb; i++) {
1190:       ISGetSize(is[i],&M);
1191:       VecGetSize(ctx->v[i],&N);
1192:       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);
1193:       ISGetLocalSize(is[i],&m);
1194:       VecGetLocalSize(ctx->v[i],&n);
1195:       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);
1196: #if defined(PETSC_USE_DEBUG)
1197:       {                         /* This test can be expensive */
1198:         PetscInt  start;
1199:         PetscBool contiguous;
1200:         ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1201:         if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1202:         if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1203:       }
1204: #endif
1205:       PetscObjectReference((PetscObject)is[i]);
1206:       ctx->is[i] = is[i];
1207:       offset += n;
1208:     }
1209:   } else {                      /* Create a contiguous ISStride for each entry */
1210:     offset = V->map->rstart;
1211:     for (i=0; i<ctx->nb; i++) {
1212:       PetscInt bs;
1213:       VecGetLocalSize(ctx->v[i],&n);
1214:       VecGetBlockSize(ctx->v[i],&bs);
1215:       ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1216:       ISSetBlockSize(ctx->is[i],bs);
1217:       offset += n;
1218:     }
1219:   }
1220:   return(0);
1221: }

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

1228:    Collective on Vec

1230:    Input Parameter:
1231: +  comm - Communicator for the new Vec
1232: .  nb - number of nested blocks
1233: .  is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1234: -  x - array of nb sub-vectors

1236:    Output Parameter:
1237: .  Y - new vector

1239:    Level: advanced

1241: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1242: @*/
1243: PetscErrorCode  VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1244: {
1245:   Vec            V;
1246:   Vec_Nest       *s;
1247:   PetscInt       n,N;

1251:   VecCreate(comm,&V);
1252:   PetscObjectChangeTypeName((PetscObject)V,VECNEST);

1254:   /* allocate and set pointer for implememtation data */
1255:   PetscMalloc(sizeof(Vec_Nest),&s);
1256:   PetscMemzero(s,sizeof(Vec_Nest));
1257:   V->data          = (void*)s;
1258:   s->setup_called  = PETSC_FALSE;
1259:   s->nb            = -1;
1260:   s->v             = NULL;

1262:   VecSetUp_Nest_Private(V,nb,x);

1264:   n = N = 0;
1265:   VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1266:   VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1267:   PetscLayoutSetSize(V->map,N);
1268:   PetscLayoutSetLocalSize(V->map,n);
1269:   PetscLayoutSetBlockSize(V->map,1);
1270:   PetscLayoutSetUp(V->map);

1272:   VecSetUp_NestIS_Private(V,nb,is);

1274:   VecNestSetOps_Private(V->ops);
1275:   V->petscnative = PETSC_FALSE;


1278:   /* expose block api's */
1279:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1280:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1281:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1282:   PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1283:   PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);

1285:   *Y = V;
1286:   return(0);
1287: }

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

1292:   Level: intermediate

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

1298: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1299: M*/