Actual source code: vecnest.c
petsc-3.3-p7 2013-05-11
2: #include 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(((PetscObject)v)->comm,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: PetscObjectComposeFunctionDynamic((PetscObject)v,"","",PETSC_NULL);
56: PetscObjectComposeFunctionDynamic((PetscObject)v,"","",PETSC_NULL);
57: PetscObjectComposeFunctionDynamic((PetscObject)v,"","",PETSC_NULL);
58: PetscObjectComposeFunctionDynamic((PetscObject)v,"","",PETSC_NULL);
59: PetscObjectComposeFunctionDynamic((PetscObject)v,"","",PETSC_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;
96: PetscMalloc(sizeof(Vec)*bx->nb,&sub);
97: for (i=0; i<bx->nb; i++ ) {
98: VecDuplicate(bx->v[i],&sub[i]);
99: }
100: VecCreateNest(((PetscObject)x)->comm,bx->nb,bx->is,sub,&Y);
101: for (i=0; i<bx->nb; i++ ) {
102: VecDestroy(&sub[i]);
103: }
104: PetscFree(sub);
105: *y = Y;
106: return(0);
107: }
109: /* supports nested blocks */
112: static PetscErrorCode VecDot_Nest(Vec x,Vec y,PetscScalar *val)
113: {
114: Vec_Nest *bx = (Vec_Nest*)x->data;
115: Vec_Nest *by = (Vec_Nest*)y->data;
116: PetscInt i,nr;
117: PetscScalar x_dot_y,_val;
121: nr = bx->nb;
122: _val = 0.0;
123: for (i=0; i<nr; i++) {
124: VecDot(bx->v[i],by->v[i],&x_dot_y);
125: _val = _val + x_dot_y;
126: }
127: *val = _val;
128: return(0);
129: }
131: /* supports nested blocks */
134: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
135: {
136: Vec_Nest *bx = (Vec_Nest*)x->data;
137: Vec_Nest *by = (Vec_Nest*)y->data;
138: PetscInt i,nr;
139: PetscScalar x_dot_y,_val;
143: nr = bx->nb;
144: _val = 0.0;
145: for (i=0; i<nr; i++) {
146: VecTDot(bx->v[i],by->v[i],&x_dot_y);
147: _val = _val + x_dot_y;
148: }
149: *val = _val;
150: return(0);
151: }
155: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
156: {
157: Vec_Nest *bx = (Vec_Nest*)x->data;
158: Vec_Nest *by = (Vec_Nest*)y->data;
159: PetscInt i,nr;
160: PetscScalar x_dot_y,norm2_y,_dp,_nm;
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: }
297: return(0);
298: }
302: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal* z)
303: {
304: Vec_Nest *bx = (Vec_Nest*)xin->data;
305: PetscInt i,nr;
306: PetscReal z_i;
307: PetscReal _z;
311: nr = bx->nb;
312: _z = 0.0;
314: if (type == NORM_2) {
315: PetscScalar dot;
316: VecDot(xin,xin,&dot);
317: _z = PetscAbsScalar(PetscSqrtScalar(dot));
318: } else if (type == NORM_1) {
319: for (i=0; i<nr; i++) {
320: VecNorm(bx->v[i],type,&z_i);
321: _z = _z + z_i;
322: }
323: } else if (type == NORM_INFINITY) {
324: for (i=0; i<nr; i++) {
325: VecNorm(bx->v[i],type,&z_i);
326: if (z_i > _z) {
327: _z = z_i;
328: }
329: }
330: }
332: *z = _z;
333: return(0);
334: }
338: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
339: {
340: PetscInt v;
344: for (v=0; v<nv; v++) {
345: /* Do axpy on each vector,v */
346: VecAXPY(y,alpha[v],x[v]);
347: }
348: return(0);
349: }
353: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
354: {
355: PetscInt j;
359: for (j=0; j<nv; j++) {
360: VecDot(x,y[j],&val[j]);
361: }
362: return(0);
363: }
367: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
368: {
369: PetscInt j;
373: for (j=0; j<nv; j++) {
374: VecTDot(x,y[j],&val[j]);
375: }
376: return(0);
377: }
381: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
382: {
383: Vec_Nest *bx = (Vec_Nest*)x->data;
384: PetscInt i,nr;
388: nr = bx->nb;
389: for (i=0; i<nr; i++) {
390: VecSet(bx->v[i],alpha);
391: }
392: return(0);
393: }
397: static PetscErrorCode VecConjugate_Nest(Vec x)
398: {
399: Vec_Nest *bx = (Vec_Nest*)x->data;
400: PetscInt j,nr;
404: nr = bx->nb;
405: for (j=0; j<nr; j++) {
406: VecConjugate(bx->v[j]);
407: }
408: return(0);
409: }
413: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
414: {
415: Vec_Nest *bx = (Vec_Nest*)x->data;
416: Vec_Nest *by = (Vec_Nest*)y->data;
417: PetscInt i,nr;
421: VecNestCheckCompatible2(x,1,y,2);
422: nr = bx->nb;
423: for (i=0; i<nr; i++) {
424: VecSwap(bx->v[i],by->v[i]);
425: }
426: return(0);
427: }
431: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
432: {
433: Vec_Nest *bx = (Vec_Nest*)x->data;
434: Vec_Nest *by = (Vec_Nest*)y->data;
435: Vec_Nest *bw = (Vec_Nest*)w->data;
436: PetscInt i,nr;
440: VecNestCheckCompatible3(w,1,x,3,y,4);
442: nr = bx->nb;
443: for (i=0; i<nr; i++) {
444: VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
445: }
446: return(0);
447: }
451: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
452: {
453: Vec_Nest *bx = (Vec_Nest*)x->data;
455: PetscInt i,nr;
456: PetscBool isnest;
457: PetscInt L;
458: PetscInt _entry_loc;
459: PetscReal _entry_val;
463: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
464: if (!isnest) {
465: /* Not nest */
466: VecMax(x,&_entry_loc,&_entry_val);
467: if (_entry_val > *max) {
468: *max = _entry_val;
469: *p = _entry_loc + *cnt;
470: }
471: VecGetSize(x,&L);
472: *cnt = *cnt + L;
473: return(0);
474: }
476: /* Otherwise we have a nest */
477: bx = (Vec_Nest*)x->data;
478: nr = bx->nb;
480: /* now descend recursively */
481: for (i=0; i<nr; i++) {
482: VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
483: }
484: return(0);
485: }
487: /* supports nested blocks */
490: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
491: {
492: PetscInt cnt;
496: cnt = 0;
497: *p = 0;
498: *max = 0.0;
499: VecMax_Nest_Recursive(x,&cnt,p,max);
500: return(0);
501: }
505: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
506: {
507: Vec_Nest *bx = (Vec_Nest*)x->data;
508: PetscInt i,nr,L,_entry_loc;
509: PetscBool isnest;
510: PetscReal _entry_val;
514: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
515: if (!isnest) {
516: /* Not nest */
517: VecMin(x,&_entry_loc,&_entry_val);
518: if (_entry_val < *min) {
519: *min = _entry_val;
520: *p = _entry_loc + *cnt;
521: }
522: VecGetSize(x,&L);
523: *cnt = *cnt + L;
524: return(0);
525: }
527: /* Otherwise we have a nest */
528: bx = (Vec_Nest*)x->data;
529: nr = bx->nb;
531: /* now descend recursively */
532: for (i=0; i<nr; i++) {
533: VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
534: }
535: return(0);
536: }
540: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
541: {
542: PetscInt cnt;
546: cnt = 0;
547: *p = 0;
548: *min = 1.0e308;
549: VecMin_Nest_Recursive(x,&cnt,p,min);
550: return(0);
551: }
553: /* supports nested blocks */
556: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
557: {
558: Vec_Nest *bx = (Vec_Nest*)x->data;
559: PetscBool isascii;
560: PetscInt i;
564: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
565: if (isascii) {
566: PetscViewerASCIIPrintf(viewer,"Vector Object:\n");
567: PetscViewerASCIIPushTab(viewer); /* push0 */
568: PetscViewerASCIIPrintf(viewer,"type=nest, rows=%d \n",bx->nb);
570: PetscViewerASCIIPrintf(viewer,"VecNest structure: \n");
571: for (i=0; i<bx->nb; i++) {
572: const VecType type;
573: char name[256] = "",prefix[256] = "";
574: PetscInt NR;
576: VecGetSize(bx->v[i],&NR);
577: VecGetType(bx->v[i],&type);
578: if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof name,"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
579: if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof prefix,"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}
581: PetscViewerASCIIPrintf(viewer,"(%D) : %s%stype=%s, rows=%D \n",i,name,prefix,type,NR);
583: PetscViewerASCIIPushTab(viewer); /* push1 */
584: VecView(bx->v[i],viewer);
585: PetscViewerASCIIPopTab(viewer); /* pop1 */
586: }
587: PetscViewerASCIIPopTab(viewer); /* pop0 */
588: }
589: return(0);
590: }
594: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
595: {
596: Vec_Nest *bx = (Vec_Nest*)x->data;
597: PetscInt size,i,nr;
598: PetscBool isnest;
602: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
603: if (!isnest) {
604: /* Not nest */
605: if (globalsize) { VecGetSize(x,&size); }
606: else { VecGetLocalSize(x,&size); }
607: *L = *L + size;
608: return(0);
609: }
611: /* Otherwise we have a nest */
612: bx = (Vec_Nest*)x->data;
613: nr = bx->nb;
615: /* now descend recursively */
616: for (i=0; i<nr; i++) {
617: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
618: }
619: return(0);
620: }
622: /* Returns the sum of the global size of all the consituent vectors in the nest */
625: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
626: {
628: *N = x->map->N;
629: return(0);
630: }
632: /* Returns the sum of the local size of all the consituent vectors in the nest */
635: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
636: {
638: *n = x->map->n;
639: return(0);
640: }
644: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
645: {
646: Vec_Nest *bx = (Vec_Nest*)x->data;
647: Vec_Nest *by = (Vec_Nest*)y->data;
648: PetscInt i,nr;
649: PetscReal local_max,m;
653: VecNestCheckCompatible2(x,1,y,2);
654: nr = bx->nb;
655: m = 0.0;
656: for (i=0; i<nr; i++) {
657: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
658: if (local_max > m) {
659: m = local_max;
660: }
661: }
662: *max = m;
663: return(0);
664: }
668: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
669: {
670: Vec_Nest *bx = (Vec_Nest*)X->data;
671: PetscInt i;
675: *x = PETSC_NULL;
676: for (i=0;i<bx->nb; i++) {
677: PetscBool issame = PETSC_FALSE;
678: ISEqual(is,bx->is[i],&issame);
679: if (issame) {
680: *x = bx->v[i];
681: PetscObjectReference((PetscObject)(*x));
682: break;
683: }
684: }
685: if (!*x) SETERRQ(((PetscObject)is)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
686: return(0);
687: }
691: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
692: {
696: VecDestroy(x);
697: return(0);
698: }
702: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
703: {
704: Vec_Nest *bx = (Vec_Nest*)X->data;
705: PetscInt i,m,rstart,rend;
709: VecGetOwnershipRange(X,&rstart,&rend);
710: VecGetLocalSize(X,&m);
711: PetscMalloc(m*sizeof(PetscScalar),x);
712: for (i=0; i<bx->nb; i++) {
713: Vec subvec = bx->v[i];
714: IS isy = bx->is[i];
715: PetscInt j,sm;
716: const PetscInt *ixy;
717: const PetscScalar *y;
718: VecGetLocalSize(subvec,&sm);
719: VecGetArrayRead(subvec,&y);
720: ISGetIndices(isy,&ixy);
721: for (j=0; j<sm; j++) {
722: PetscInt ix = ixy[j];
723: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
724: (*x)[ix-rstart] = y[j];
725: }
726: ISRestoreIndices(isy,&ixy);
727: VecRestoreArrayRead(subvec,&y);
728: }
729: return(0);
730: }
734: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
735: {
736: Vec_Nest *bx = (Vec_Nest*)X->data;
737: PetscInt i,m,rstart,rend;
741: VecGetOwnershipRange(X,&rstart,&rend);
742: VecGetLocalSize(X,&m);
743: for (i=0; i<bx->nb; i++) {
744: Vec subvec = bx->v[i];
745: IS isy = bx->is[i];
746: PetscInt j,sm;
747: const PetscInt *ixy;
748: PetscScalar *y;
749: VecGetLocalSize(subvec,&sm);
750: VecGetArray(subvec,&y);
751: ISGetIndices(isy,&ixy);
752: for (j=0; j<sm; j++) {
753: PetscInt ix = ixy[j];
754: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
755: y[j] = (*x)[ix-rstart];
756: }
757: ISRestoreIndices(isy,&ixy);
758: VecRestoreArray(subvec,&y);
759: }
760: PetscFree(*x);
761: return(0);
762: }
767: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
768: {
771: /* 0 */
772: ops->duplicate = VecDuplicate_Nest;
773: ops->duplicatevecs = VecDuplicateVecs_Default;
774: ops->destroyvecs = VecDestroyVecs_Default;
775: ops->dot = VecDot_Nest;
776: ops->mdot = VecMDot_Nest;
778: /* 5 */
779: ops->norm = VecNorm_Nest;
780: ops->tdot = VecTDot_Nest;
781: ops->mtdot = VecMTDot_Nest;
782: ops->scale = VecScale_Nest;
783: ops->copy = VecCopy_Nest;
785: /* 10 */
786: ops->set = VecSet_Nest;
787: ops->swap = VecSwap_Nest;
788: ops->axpy = VecAXPY_Nest;
789: ops->axpby = VecAXPBY_Nest;
790: ops->maxpy = VecMAXPY_Nest;
792: /* 15 */
793: ops->aypx = VecAYPX_Nest;
794: ops->waxpy = VecWAXPY_Nest;
795: ops->axpbypcz = 0;
796: ops->pointwisemult = VecPointwiseMult_Nest;
797: ops->pointwisedivide = VecPointwiseDivide_Nest;
798: /* 20 */
799: ops->setvalues = 0;
800: ops->assemblybegin = VecAssemblyBegin_Nest;
801: ops->assemblyend = VecAssemblyEnd_Nest;
802: ops->getarray = VecGetArray_Nest;
803: ops->getsize = VecGetSize_Nest;
805: /* 25 */
806: ops->getlocalsize = VecGetLocalSize_Nest;
807: ops->restorearray = VecRestoreArray_Nest;
808: ops->max = VecMax_Nest;
809: ops->min = VecMin_Nest;
810: ops->setrandom = 0;
812: /* 30 */
813: ops->setoption = 0;
814: ops->setvaluesblocked = 0;
815: ops->destroy = VecDestroy_Nest;
816: ops->view = VecView_Nest;
817: ops->placearray = 0;
819: /* 35 */
820: ops->replacearray = 0;
821: ops->dot_local = VecDot_Nest;
822: ops->tdot_local = VecTDot_Nest;
823: ops->norm_local = VecNorm_Nest;
824: ops->mdot_local = VecMDot_Nest;
826: /* 40 */
827: ops->mtdot_local = VecMTDot_Nest;
828: ops->load = 0;
829: ops->reciprocal = VecReciprocal_Nest;
830: ops->conjugate = VecConjugate_Nest;
831: ops->setlocaltoglobalmapping = 0;
833: /* 45 */
834: ops->setvalueslocal = 0;
835: ops->resetarray = 0;
836: ops->setfromoptions = 0;
837: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
838: ops->load = 0;
840: /* 50 */
841: ops->pointwisemax = 0;
842: ops->pointwisemaxabs = 0;
843: ops->pointwisemin = 0;
844: ops->getvalues = 0;
845: ops->sqrt = 0;
847: /* 55 */
848: ops->abs = 0;
849: ops->exp = 0;
850: ops->shift = 0;
851: ops->create = 0;
852: ops->stridegather = 0;
854: /* 60 */
855: ops->stridescatter = 0;
856: ops->dotnorm2 = VecDotNorm2_Nest;
857: ops->getsubvector = VecGetSubVector_Nest;
858: ops->restoresubvector = VecRestoreSubVector_Nest;
860: return(0);
861: }
865: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
866: {
867: Vec_Nest *b = (Vec_Nest*)x->data;
868: PetscInt i;
869: PetscInt row;
872: if (!m) return(0);
873: for (i=0; i<m; i++) {
874: row = idxm[i];
875: if (row >= b->nb) SETERRQ2(((PetscObject)x)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
876: vec[i] = b->v[row];
877: }
878: return(0);
879: }
881: EXTERN_C_BEGIN
884: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
885: {
889: VecNestGetSubVecs_Private(X,1,&idxm,sx);
890: return(0);
891: }
892: EXTERN_C_END
896: /*@C
897: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
899: Not collective
901: Input Parameters:
902: . X - nest vector
903: . idxm - index of the vector within the nest
905: Output Parameter:
906: . sx - vector at index idxm within the nest
908: Notes:
910: Level: developer
912: .seealso: VecNestGetSize(), VecNestGetSubVecs()
913: @*/
914: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
915: {
919: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
920: return(0);
921: }
923: EXTERN_C_BEGIN
926: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
927: {
928: Vec_Nest *b = (Vec_Nest*)X->data;
930: if (N) { *N = b->nb; }
931: if (sx) { *sx = b->v; }
932: return(0);
933: }
934: EXTERN_C_END
938: /*@C
939: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
941: Not collective
943: Input Parameters:
944: . X - nest vector
946: Output Parameter:
947: . N - number of nested vecs
948: . sx - array of vectors
950: Notes:
952: The user should not free the array sx.
954: Level: developer
956: .seealso: VecNestGetSize(), VecNestGetSubVec()
957: @*/
958: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
959: {
963: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
964: return(0);
965: }
969: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
970: {
971: Vec_Nest *bx = (Vec_Nest*)X->data;
972: PetscInt i,offset=0,n=0,bs;
973: IS is;
975: PetscBool issame = PETSC_FALSE;
976: PetscInt N=0;
977:
978: /* check if idxm < bx->nb */
979: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
980:
982:
983: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
984: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
985: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
986:
987: /* check if we need to update the IS for the block */
988: offset = X->map->rstart;
989: for (i=0; i<idxm; i++) {
990: n=0;
991: VecGetLocalSize(bx->v[i],&n);
992: offset += n;
993: }
994:
995: /* get the local size and block size */
996: VecGetLocalSize(x,&n);
997: VecGetBlockSize(x,&bs);
998:
999: /* create the new IS */
1000: ISCreateStride(((PetscObject)x)->comm,n,offset,1,&is);
1001: ISSetBlockSize(is,bs);
1002:
1003: /* check if they are equal */
1004: ISEqual(is,bx->is[idxm],&issame);
1005:
1006: if (!issame) {
1007: /* The IS of given vector has a different layout compared to the existing block vector.
1008: Destroy the existing reference and update the IS. */
1009: ISDestroy(&bx->is[idxm]);
1010: ISDuplicate(is,&bx->is[idxm]);
1011: ISCopy(is,bx->is[idxm]);
1012:
1013: offset += n;
1014: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
1015: for (i=idxm+1; i<bx->nb; i++) {
1016: /* get the local size and block size */
1017: VecGetLocalSize(bx->v[i],&n);
1018: VecGetBlockSize(bx->v[i],&bs);
1019:
1020: /* destroy the old and create the new IS */
1021: ISDestroy(&bx->is[i]);
1022: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
1023: ISSetBlockSize(bx->is[i],bs);
1024:
1025: offset += n;
1026: }
1027:
1028: n=0;
1029: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
1030: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
1031: PetscLayoutSetSize(X->map,N);
1032: PetscLayoutSetLocalSize(X->map,n);
1033: }
1034:
1035: ISDestroy(&is);
1036:
1037: return(0);
1038: }
1040: EXTERN_C_BEGIN
1043: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
1044: {
1046:
1048: VecNestSetSubVec_Private(X,idxm,sx);
1049: return(0);
1050: }
1051: EXTERN_C_END
1055: /*@C
1056: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
1057:
1058: Not collective
1059:
1060: Input Parameters:
1061: . X - nest vector
1062: . idxm - index of the vector within the nest vector
1063: . sx - vector at index idxm within the nest vector
1065: Notes:
1066:
1067: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1068:
1069: Level: developer
1070:
1071: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
1072: @*/
1073: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
1074: {
1076:
1078: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
1079: return(0);
1080: }
1082: EXTERN_C_BEGIN
1085: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1086: {
1087: PetscInt i;
1089:
1091: for (i=0; i<N; i++) {
1092: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
1093: }
1095: return(0);
1096: }
1097: EXTERN_C_END
1101: /*@C
1102: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1103:
1104: Not collective
1105:
1106: Input Parameters:
1107: . X - nest vector
1108: . N - number of component vecs in sx
1109: . idxm - indices of component vecs that are to be replaced
1110: . sx - array of vectors
1111:
1112: Notes:
1113:
1114: The components in the vector array sx do not have to be of the same size as corresponding
1115: components in X. The user can also free the array "sx" after the call.
1116:
1117: Level: developer
1118:
1119: .seealso: VecNestGetSize(), VecNestGetSubVec()
1120: @*/
1121: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1122: {
1124:
1126: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1127: return(0);
1128: }
1130: EXTERN_C_BEGIN
1133: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1134: {
1135: Vec_Nest *b = (Vec_Nest*)X->data;
1138: *N = b->nb;
1139: return(0);
1140: }
1141: EXTERN_C_END
1145: /*@C
1146: VecNestGetSize - Returns the size of the nest vector.
1148: Not collective
1150: Input Parameters:
1151: . X - nest vector
1153: Output Parameter:
1154: . N - number of nested vecs
1156: Notes:
1158: Level: developer
1160: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1161: @*/
1162: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1163: {
1169: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1170: return(0);
1171: }
1175: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1176: {
1177: Vec_Nest *ctx = (Vec_Nest*)V->data;
1178: PetscInt i;
1182: if (ctx->setup_called) return(0);
1184: ctx->nb = nb;
1185: if (ctx->nb < 0) SETERRQ(((PetscObject)V)->comm,PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1187: /* Create space */
1188: PetscMalloc(ctx->nb*sizeof(Vec),&ctx->v);
1189: for (i=0; i<ctx->nb; i++) {
1190: ctx->v[i] = x[i];
1191: PetscObjectReference((PetscObject)x[i]);
1192: /* Do not allocate memory for internal sub blocks */
1193: }
1195: PetscMalloc(ctx->nb*sizeof(IS),&ctx->is);
1197: ctx->setup_called = PETSC_TRUE;
1198: return(0);
1199: }
1203: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1204: {
1205: Vec_Nest *ctx = (Vec_Nest*)V->data;
1206: PetscInt i,offset,m,n,M,N;
1210: if (is) { /* Do some consistency checks and reference the is */
1211: offset = V->map->rstart;
1212: for (i=0; i<ctx->nb; i++) {
1213: ISGetSize(is[i],&M);
1214: VecGetSize(ctx->v[i],&N);
1215: if (M != N) SETERRQ3(((PetscObject)V)->comm,PETSC_ERR_ARG_INCOMP,"In slot %D, IS of size %D is not compatible with Vec of size %D",i,M,N);
1216: ISGetLocalSize(is[i],&m);
1217: VecGetLocalSize(ctx->v[i],&n);
1218: 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);
1219: #if defined(PETSC_USE_DEBUG)
1220: { /* This test can be expensive */
1221: PetscInt start;
1222: PetscBool contiguous;
1223: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1224: if (!contiguous) SETERRQ1(((PetscObject)V)->comm,PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1225: if (start != 0) SETERRQ1(((PetscObject)V)->comm,PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1226: }
1227: #endif
1228: PetscObjectReference((PetscObject)is[i]);
1229: ctx->is[i] = is[i];
1230: offset += n;
1231: }
1232: } else { /* Create a contiguous ISStride for each entry */
1233: offset = V->map->rstart;
1234: for (i=0; i<ctx->nb; i++) {
1235: PetscInt bs;
1236: VecGetLocalSize(ctx->v[i],&n);
1237: VecGetBlockSize(ctx->v[i],&bs);
1238: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1239: ISSetBlockSize(ctx->is[i],bs);
1240: offset += n;
1241: }
1242: }
1243: return(0);
1244: }
1248: /*@
1249: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1251: Collective on Vec
1253: Input Parameter:
1254: + comm - Communicator for the new Vec
1255: . nb - number of nested blocks
1256: . is - array of nb index sets describing each nested block, or PETSC_NULL to pack subvectors contiguously
1257: - x - array of nb sub-vectors
1259: Output Parameter:
1260: . Y - new vector
1262: Level: advanced
1264: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1265: @*/
1266: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1267: {
1268: Vec V;
1269: Vec_Nest *s;
1270: PetscInt n,N;
1274: VecCreate(comm,&V);
1275: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1277: /* allocate and set pointer for implememtation data */
1278: PetscMalloc(sizeof(Vec_Nest),&s);
1279: PetscMemzero(s,sizeof(Vec_Nest));
1280: V->data = (void*)s;
1281: s->setup_called = PETSC_FALSE;
1282: s->nb = -1;
1283: s->v = PETSC_NULL;
1285: VecSetUp_Nest_Private(V,nb,x);
1287: n = N = 0;
1288: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1289: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1290: PetscLayoutSetSize(V->map,N);
1291: PetscLayoutSetLocalSize(V->map,n);
1292: PetscLayoutSetBlockSize(V->map,1);
1293: PetscLayoutSetUp(V->map);
1295: VecSetUp_NestIS_Private(V,nb,is);
1297: VecNestSetOps_Private(V->ops);
1298: V->petscnative = PETSC_FALSE;
1301: /* expose block api's */
1302: PetscObjectComposeFunctionDynamic((PetscObject)V,"VecNestGetSubVec_C","VecNestGetSubVec_Nest",VecNestGetSubVec_Nest);
1303: PetscObjectComposeFunctionDynamic((PetscObject)V,"VecNestGetSubVecs_C","VecNestGetSubVecs_Nest",VecNestGetSubVecs_Nest);
1304: PetscObjectComposeFunctionDynamic((PetscObject)V,"VecNestSetSubVec_C","VecNestSetSubVec_Nest",VecNestSetSubVec_Nest);
1305: PetscObjectComposeFunctionDynamic((PetscObject)V,"VecNestSetSubVecs_C","VecNestSetSubVecs_Nest",VecNestSetSubVecs_Nest);
1306: PetscObjectComposeFunctionDynamic((PetscObject)V,"VecNestGetSize_C","VecNestGetSize_Nest",VecNestGetSize_Nest);
1308: *Y = V;
1309: return(0);
1310: }
1312: /*MC
1313: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1315: Level: intermediate
1317: Notes:
1318: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1319: It is usually used with MATNEST and DMComposite via DMSetVecType().
1321: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1322: M*/