Actual source code: vecnest.c
petsc-3.4.5 2014-06-29
2: #include ../src/vec/vec/impls/nest/vecnestimpl.h
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: PetscViewerASCIIPrintf(viewer,"Vector Object:\n");
563: PetscViewerASCIIPushTab(viewer); /* push0 */
564: PetscViewerASCIIPrintf(viewer,"type=nest, rows=%d \n",bx->nb);
566: PetscViewerASCIIPrintf(viewer,"VecNest structure: \n");
567: for (i=0; i<bx->nb; i++) {
568: VecType type;
569: char name[256] = "",prefix[256] = "";
570: PetscInt NR;
572: VecGetSize(bx->v[i],&NR);
573: VecGetType(bx->v[i],&type);
574: if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
575: if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}
577: PetscViewerASCIIPrintf(viewer,"(%D) : %s%stype=%s, rows=%D \n",i,name,prefix,type,NR);
579: PetscViewerASCIIPushTab(viewer); /* push1 */
580: VecView(bx->v[i],viewer);
581: PetscViewerASCIIPopTab(viewer); /* pop1 */
582: }
583: PetscViewerASCIIPopTab(viewer); /* pop0 */
584: }
585: return(0);
586: }
590: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
591: {
592: Vec_Nest *bx;
593: PetscInt size,i,nr;
594: PetscBool isnest;
598: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
599: if (!isnest) {
600: /* Not nest */
601: if (globalsize) { VecGetSize(x,&size); }
602: else { VecGetLocalSize(x,&size); }
603: *L = *L + size;
604: return(0);
605: }
607: /* Otherwise we have a nest */
608: bx = (Vec_Nest*)x->data;
609: nr = bx->nb;
611: /* now descend recursively */
612: for (i=0; i<nr; i++) {
613: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
614: }
615: return(0);
616: }
618: /* Returns the sum of the global size of all the consituent vectors in the nest */
621: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
622: {
624: *N = x->map->N;
625: return(0);
626: }
628: /* Returns the sum of the local size of all the consituent vectors in the nest */
631: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
632: {
634: *n = x->map->n;
635: return(0);
636: }
640: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
641: {
642: Vec_Nest *bx = (Vec_Nest*)x->data;
643: Vec_Nest *by = (Vec_Nest*)y->data;
644: PetscInt i,nr;
645: PetscReal local_max,m;
649: VecNestCheckCompatible2(x,1,y,2);
650: nr = bx->nb;
651: m = 0.0;
652: for (i=0; i<nr; i++) {
653: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
654: if (local_max > m) m = local_max;
655: }
656: *max = m;
657: return(0);
658: }
662: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
663: {
664: Vec_Nest *bx = (Vec_Nest*)X->data;
665: PetscInt i;
669: *x = NULL;
670: for (i=0; i<bx->nb; i++) {
671: PetscBool issame = PETSC_FALSE;
672: ISEqual(is,bx->is[i],&issame);
673: if (issame) {
674: *x = bx->v[i];
675: PetscObjectReference((PetscObject)(*x));
676: break;
677: }
678: }
679: if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
680: return(0);
681: }
685: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
686: {
690: VecDestroy(x);
691: return(0);
692: }
696: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
697: {
698: Vec_Nest *bx = (Vec_Nest*)X->data;
699: PetscInt i,m,rstart,rend;
703: VecGetOwnershipRange(X,&rstart,&rend);
704: VecGetLocalSize(X,&m);
705: PetscMalloc(m*sizeof(PetscScalar),x);
706: for (i=0; i<bx->nb; i++) {
707: Vec subvec = bx->v[i];
708: IS isy = bx->is[i];
709: PetscInt j,sm;
710: const PetscInt *ixy;
711: const PetscScalar *y;
712: VecGetLocalSize(subvec,&sm);
713: VecGetArrayRead(subvec,&y);
714: ISGetIndices(isy,&ixy);
715: for (j=0; j<sm; j++) {
716: PetscInt ix = ixy[j];
717: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
718: (*x)[ix-rstart] = y[j];
719: }
720: ISRestoreIndices(isy,&ixy);
721: VecRestoreArrayRead(subvec,&y);
722: }
723: return(0);
724: }
728: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
729: {
730: Vec_Nest *bx = (Vec_Nest*)X->data;
731: PetscInt i,m,rstart,rend;
735: VecGetOwnershipRange(X,&rstart,&rend);
736: VecGetLocalSize(X,&m);
737: for (i=0; i<bx->nb; i++) {
738: Vec subvec = bx->v[i];
739: IS isy = bx->is[i];
740: PetscInt j,sm;
741: const PetscInt *ixy;
742: PetscScalar *y;
743: VecGetLocalSize(subvec,&sm);
744: VecGetArray(subvec,&y);
745: ISGetIndices(isy,&ixy);
746: for (j=0; j<sm; j++) {
747: PetscInt ix = ixy[j];
748: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
749: y[j] = (*x)[ix-rstart];
750: }
751: ISRestoreIndices(isy,&ixy);
752: VecRestoreArray(subvec,&y);
753: }
754: PetscFree(*x);
755: return(0);
756: }
761: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
762: {
764: /* 0 */
765: ops->duplicate = VecDuplicate_Nest;
766: ops->duplicatevecs = VecDuplicateVecs_Default;
767: ops->destroyvecs = VecDestroyVecs_Default;
768: ops->dot = VecDot_Nest;
769: ops->mdot = VecMDot_Nest;
771: /* 5 */
772: ops->norm = VecNorm_Nest;
773: ops->tdot = VecTDot_Nest;
774: ops->mtdot = VecMTDot_Nest;
775: ops->scale = VecScale_Nest;
776: ops->copy = VecCopy_Nest;
778: /* 10 */
779: ops->set = VecSet_Nest;
780: ops->swap = VecSwap_Nest;
781: ops->axpy = VecAXPY_Nest;
782: ops->axpby = VecAXPBY_Nest;
783: ops->maxpy = VecMAXPY_Nest;
785: /* 15 */
786: ops->aypx = VecAYPX_Nest;
787: ops->waxpy = VecWAXPY_Nest;
788: ops->axpbypcz = 0;
789: ops->pointwisemult = VecPointwiseMult_Nest;
790: ops->pointwisedivide = VecPointwiseDivide_Nest;
791: /* 20 */
792: ops->setvalues = 0;
793: ops->assemblybegin = VecAssemblyBegin_Nest;
794: ops->assemblyend = VecAssemblyEnd_Nest;
795: ops->getarray = VecGetArray_Nest;
796: ops->getsize = VecGetSize_Nest;
798: /* 25 */
799: ops->getlocalsize = VecGetLocalSize_Nest;
800: ops->restorearray = VecRestoreArray_Nest;
801: ops->max = VecMax_Nest;
802: ops->min = VecMin_Nest;
803: ops->setrandom = 0;
805: /* 30 */
806: ops->setoption = 0;
807: ops->setvaluesblocked = 0;
808: ops->destroy = VecDestroy_Nest;
809: ops->view = VecView_Nest;
810: ops->placearray = 0;
812: /* 35 */
813: ops->replacearray = 0;
814: ops->dot_local = VecDot_Nest;
815: ops->tdot_local = VecTDot_Nest;
816: ops->norm_local = VecNorm_Nest;
817: ops->mdot_local = VecMDot_Nest;
819: /* 40 */
820: ops->mtdot_local = VecMTDot_Nest;
821: ops->load = 0;
822: ops->reciprocal = VecReciprocal_Nest;
823: ops->conjugate = VecConjugate_Nest;
824: ops->setlocaltoglobalmapping = 0;
826: /* 45 */
827: ops->setvalueslocal = 0;
828: ops->resetarray = 0;
829: ops->setfromoptions = 0;
830: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
831: ops->load = 0;
833: /* 50 */
834: ops->pointwisemax = 0;
835: ops->pointwisemaxabs = 0;
836: ops->pointwisemin = 0;
837: ops->getvalues = 0;
838: ops->sqrt = 0;
840: /* 55 */
841: ops->abs = 0;
842: ops->exp = 0;
843: ops->shift = 0;
844: ops->create = 0;
845: ops->stridegather = 0;
847: /* 60 */
848: ops->stridescatter = 0;
849: ops->dotnorm2 = VecDotNorm2_Nest;
850: ops->getsubvector = VecGetSubVector_Nest;
851: ops->restoresubvector = VecRestoreSubVector_Nest;
852: return(0);
853: }
857: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
858: {
859: Vec_Nest *b = (Vec_Nest*)x->data;
860: PetscInt i;
861: PetscInt row;
864: if (!m) return(0);
865: for (i=0; i<m; i++) {
866: row = idxm[i];
867: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
868: vec[i] = b->v[row];
869: }
870: return(0);
871: }
875: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
876: {
880: VecNestGetSubVecs_Private(X,1,&idxm,sx);
881: return(0);
882: }
886: /*@
887: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
889: Not collective
891: Input Parameters:
892: . X - nest vector
893: . idxm - index of the vector within the nest
895: Output Parameter:
896: . sx - vector at index idxm within the nest
898: Notes:
900: Level: developer
902: .seealso: VecNestGetSize(), VecNestGetSubVecs()
903: @*/
904: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
905: {
909: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
910: return(0);
911: }
915: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
916: {
917: Vec_Nest *b = (Vec_Nest*)X->data;
920: if (N) *N = b->nb;
921: if (sx) *sx = b->v;
922: return(0);
923: }
927: /*@C
928: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
930: Not collective
932: Input Parameters:
933: . X - nest vector
935: Output Parameter:
936: + N - number of nested vecs
937: - sx - array of vectors
939: Notes:
940: The user should not free the array sx.
942: Fortran Notes:
943: The caller must allocate the array to hold the subvectors.
945: Level: developer
947: .seealso: VecNestGetSize(), VecNestGetSubVec()
948: @*/
949: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
950: {
954: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
955: return(0);
956: }
960: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
961: {
962: Vec_Nest *bx = (Vec_Nest*)X->data;
963: PetscInt i,offset=0,n=0,bs;
964: IS is;
966: PetscBool issame = PETSC_FALSE;
967: PetscInt N=0;
969: /* check if idxm < bx->nb */
970: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
973: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
974: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
975: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
977: /* check if we need to update the IS for the block */
978: offset = X->map->rstart;
979: for (i=0; i<idxm; i++) {
980: n=0;
981: VecGetLocalSize(bx->v[i],&n);
982: offset += n;
983: }
985: /* get the local size and block size */
986: VecGetLocalSize(x,&n);
987: VecGetBlockSize(x,&bs);
989: /* create the new IS */
990: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
991: ISSetBlockSize(is,bs);
993: /* check if they are equal */
994: ISEqual(is,bx->is[idxm],&issame);
996: if (!issame) {
997: /* The IS of given vector has a different layout compared to the existing block vector.
998: Destroy the existing reference and update the IS. */
999: ISDestroy(&bx->is[idxm]);
1000: ISDuplicate(is,&bx->is[idxm]);
1001: ISCopy(is,bx->is[idxm]);
1003: offset += n;
1004: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
1005: for (i=idxm+1; i<bx->nb; i++) {
1006: /* get the local size and block size */
1007: VecGetLocalSize(bx->v[i],&n);
1008: VecGetBlockSize(bx->v[i],&bs);
1010: /* destroy the old and create the new IS */
1011: ISDestroy(&bx->is[i]);
1012: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
1013: ISSetBlockSize(bx->is[i],bs);
1015: offset += n;
1016: }
1018: n=0;
1019: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
1020: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
1021: PetscLayoutSetSize(X->map,N);
1022: PetscLayoutSetLocalSize(X->map,n);
1023: }
1025: ISDestroy(&is);
1026: return(0);
1027: }
1031: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
1032: {
1036: VecNestSetSubVec_Private(X,idxm,sx);
1037: return(0);
1038: }
1042: /*@
1043: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
1045: Not collective
1047: Input Parameters:
1048: + X - nest vector
1049: . idxm - index of the vector within the nest vector
1050: - sx - vector at index idxm within the nest vector
1052: Notes:
1053: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1055: Level: developer
1057: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
1058: @*/
1059: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
1060: {
1064: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
1065: return(0);
1066: }
1070: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1071: {
1072: PetscInt i;
1076: for (i=0; i<N; i++) {
1077: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
1078: }
1079: return(0);
1080: }
1084: /*@C
1085: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1087: Not collective
1089: Input Parameters:
1090: + X - nest vector
1091: . N - number of component vecs in sx
1092: . idxm - indices of component vecs that are to be replaced
1093: - sx - array of vectors
1095: Notes:
1096: The components in the vector array sx do not have to be of the same size as corresponding
1097: components in X. The user can also free the array "sx" after the call.
1099: Level: developer
1101: .seealso: VecNestGetSize(), VecNestGetSubVec()
1102: @*/
1103: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1104: {
1108: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1109: return(0);
1110: }
1114: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1115: {
1116: Vec_Nest *b = (Vec_Nest*)X->data;
1119: *N = b->nb;
1120: return(0);
1121: }
1125: /*@
1126: VecNestGetSize - Returns the size of the nest vector.
1128: Not collective
1130: Input Parameters:
1131: . X - nest vector
1133: Output Parameter:
1134: . N - number of nested vecs
1136: Notes:
1138: Level: developer
1140: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1141: @*/
1142: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1143: {
1149: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1150: return(0);
1151: }
1155: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1156: {
1157: Vec_Nest *ctx = (Vec_Nest*)V->data;
1158: PetscInt i;
1162: if (ctx->setup_called) return(0);
1164: ctx->nb = nb;
1165: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1167: /* Create space */
1168: PetscMalloc(ctx->nb*sizeof(Vec),&ctx->v);
1169: for (i=0; i<ctx->nb; i++) {
1170: ctx->v[i] = x[i];
1171: PetscObjectReference((PetscObject)x[i]);
1172: /* Do not allocate memory for internal sub blocks */
1173: }
1175: PetscMalloc(ctx->nb*sizeof(IS),&ctx->is);
1177: ctx->setup_called = PETSC_TRUE;
1178: return(0);
1179: }
1183: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1184: {
1185: Vec_Nest *ctx = (Vec_Nest*)V->data;
1186: PetscInt i,offset,m,n,M,N;
1190: if (is) { /* Do some consistency checks and reference the is */
1191: offset = V->map->rstart;
1192: for (i=0; i<ctx->nb; i++) {
1193: ISGetSize(is[i],&M);
1194: VecGetSize(ctx->v[i],&N);
1195: 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);
1196: ISGetLocalSize(is[i],&m);
1197: VecGetLocalSize(ctx->v[i],&n);
1198: 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);
1199: #if defined(PETSC_USE_DEBUG)
1200: { /* This test can be expensive */
1201: PetscInt start;
1202: PetscBool contiguous;
1203: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1204: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1205: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1206: }
1207: #endif
1208: PetscObjectReference((PetscObject)is[i]);
1209: ctx->is[i] = is[i];
1210: offset += n;
1211: }
1212: } else { /* Create a contiguous ISStride for each entry */
1213: offset = V->map->rstart;
1214: for (i=0; i<ctx->nb; i++) {
1215: PetscInt bs;
1216: VecGetLocalSize(ctx->v[i],&n);
1217: VecGetBlockSize(ctx->v[i],&bs);
1218: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1219: ISSetBlockSize(ctx->is[i],bs);
1220: offset += n;
1221: }
1222: }
1223: return(0);
1224: }
1228: /*@C
1229: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1231: Collective on Vec
1233: Input Parameter:
1234: + comm - Communicator for the new Vec
1235: . nb - number of nested blocks
1236: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1237: - x - array of nb sub-vectors
1239: Output Parameter:
1240: . Y - new vector
1242: Level: advanced
1244: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1245: @*/
1246: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1247: {
1248: Vec V;
1249: Vec_Nest *s;
1250: PetscInt n,N;
1254: VecCreate(comm,&V);
1255: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1257: /* allocate and set pointer for implememtation data */
1258: PetscMalloc(sizeof(Vec_Nest),&s);
1259: PetscMemzero(s,sizeof(Vec_Nest));
1260: V->data = (void*)s;
1261: s->setup_called = PETSC_FALSE;
1262: s->nb = -1;
1263: s->v = NULL;
1265: VecSetUp_Nest_Private(V,nb,x);
1267: n = N = 0;
1268: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1269: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1270: PetscLayoutSetSize(V->map,N);
1271: PetscLayoutSetLocalSize(V->map,n);
1272: PetscLayoutSetBlockSize(V->map,1);
1273: PetscLayoutSetUp(V->map);
1275: VecSetUp_NestIS_Private(V,nb,is);
1277: VecNestSetOps_Private(V->ops);
1278: V->petscnative = PETSC_FALSE;
1281: /* expose block api's */
1282: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1283: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1284: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1285: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1286: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);
1288: *Y = V;
1289: return(0);
1290: }
1292: /*MC
1293: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1295: Level: intermediate
1297: Notes:
1298: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1299: It is usually used with MATNEST and DMComposite via DMSetVecType().
1301: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1302: M*/