Actual source code: vecnest.c
petsc-3.9.4 2018-09-11
2: #include <../src/vec/vec/impls/nest/vecnestimpl.h>
4: /* check all blocks are filled */
5: static PetscErrorCode VecAssemblyBegin_Nest(Vec v)
6: {
7: Vec_Nest *vs = (Vec_Nest*)v->data;
8: PetscInt i;
12: for (i=0; i<vs->nb; i++) {
13: if (!vs->v[i]) SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_SUP,"Nest vector cannot contain NULL blocks");
14: VecAssemblyBegin(vs->v[i]);
15: }
16: return(0);
17: }
19: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
20: {
21: Vec_Nest *vs = (Vec_Nest*)v->data;
22: PetscInt i;
26: for (i=0; i<vs->nb; i++) {
27: VecAssemblyEnd(vs->v[i]);
28: }
29: return(0);
30: }
32: static PetscErrorCode VecDestroy_Nest(Vec v)
33: {
34: Vec_Nest *vs = (Vec_Nest*)v->data;
35: PetscInt i;
39: if (vs->v) {
40: for (i=0; i<vs->nb; i++) {
41: VecDestroy(&vs->v[i]);
42: }
43: PetscFree(vs->v);
44: }
45: for (i=0; i<vs->nb; i++) {
46: ISDestroy(&vs->is[i]);
47: }
48: PetscFree(vs->is);
49: PetscObjectComposeFunction((PetscObject)v,"",NULL);
50: PetscObjectComposeFunction((PetscObject)v,"",NULL);
51: PetscObjectComposeFunction((PetscObject)v,"",NULL);
52: PetscObjectComposeFunction((PetscObject)v,"",NULL);
53: PetscObjectComposeFunction((PetscObject)v,"",NULL);
55: PetscFree(vs);
56: return(0);
57: }
59: /* supports nested blocks */
60: static PetscErrorCode VecCopy_Nest(Vec x,Vec y)
61: {
62: Vec_Nest *bx = (Vec_Nest*)x->data;
63: Vec_Nest *by = (Vec_Nest*)y->data;
64: PetscInt i;
68: VecNestCheckCompatible2(x,1,y,2);
69: for (i=0; i<bx->nb; i++) {
70: VecCopy(bx->v[i],by->v[i]);
71: }
72: return(0);
73: }
75: /* supports nested blocks */
76: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
77: {
78: Vec_Nest *bx = (Vec_Nest*)x->data;
79: Vec Y;
80: Vec *sub;
81: PetscInt i;
85: PetscMalloc1(bx->nb,&sub);
86: for (i=0; i<bx->nb; i++) {
87: VecDuplicate(bx->v[i],&sub[i]);
88: }
89: VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
90: for (i=0; i<bx->nb; i++) {
91: VecDestroy(&sub[i]);
92: }
93: PetscFree(sub);
94: *y = Y;
95: return(0);
96: }
98: /* supports nested blocks */
99: static PetscErrorCode VecDot_Nest(Vec x,Vec y,PetscScalar *val)
100: {
101: Vec_Nest *bx = (Vec_Nest*)x->data;
102: Vec_Nest *by = (Vec_Nest*)y->data;
103: PetscInt i,nr;
104: PetscScalar x_dot_y,_val;
108: nr = bx->nb;
109: _val = 0.0;
110: for (i=0; i<nr; i++) {
111: VecDot(bx->v[i],by->v[i],&x_dot_y);
112: _val = _val + x_dot_y;
113: }
114: *val = _val;
115: return(0);
116: }
118: /* supports nested blocks */
119: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
120: {
121: Vec_Nest *bx = (Vec_Nest*)x->data;
122: Vec_Nest *by = (Vec_Nest*)y->data;
123: PetscInt i,nr;
124: PetscScalar x_dot_y,_val;
128: nr = bx->nb;
129: _val = 0.0;
130: for (i=0; i<nr; i++) {
131: VecTDot(bx->v[i],by->v[i],&x_dot_y);
132: _val = _val + x_dot_y;
133: }
134: *val = _val;
135: return(0);
136: }
138: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
139: {
140: Vec_Nest *bx = (Vec_Nest*)x->data;
141: Vec_Nest *by = (Vec_Nest*)y->data;
142: PetscInt i,nr;
143: PetscScalar x_dot_y,_dp,_nm;
144: PetscReal norm2_y;
148: nr = bx->nb;
149: _dp = 0.0;
150: _nm = 0.0;
151: for (i=0; i<nr; i++) {
152: VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
153: _dp += x_dot_y;
154: _nm += norm2_y;
155: }
156: *dp = _dp;
157: *nm = _nm;
158: return(0);
159: }
161: static PetscErrorCode VecAXPY_Nest(Vec y,PetscScalar alpha,Vec x)
162: {
163: Vec_Nest *bx = (Vec_Nest*)x->data;
164: Vec_Nest *by = (Vec_Nest*)y->data;
165: PetscInt i,nr;
169: nr = bx->nb;
170: for (i=0; i<nr; i++) {
171: VecAXPY(by->v[i],alpha,bx->v[i]);
172: }
173: return(0);
174: }
176: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
177: {
178: Vec_Nest *bx = (Vec_Nest*)x->data;
179: Vec_Nest *by = (Vec_Nest*)y->data;
180: PetscInt i,nr;
184: nr = bx->nb;
185: for (i=0; i<nr; i++) {
186: VecAYPX(by->v[i],alpha,bx->v[i]);
187: }
188: return(0);
189: }
191: static PetscErrorCode VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
192: {
193: Vec_Nest *bx = (Vec_Nest*)x->data;
194: Vec_Nest *by = (Vec_Nest*)y->data;
195: PetscInt i,nr;
199: nr = bx->nb;
200: for (i=0; i<nr; i++) {
201: VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
202: }
203: return(0);
204: }
206: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
207: {
208: Vec_Nest *bx = (Vec_Nest*)x->data;
209: PetscInt i,nr;
213: nr = bx->nb;
214: for (i=0; i<nr; i++) {
215: VecScale(bx->v[i],alpha);
216: }
217: return(0);
218: }
220: static PetscErrorCode VecPointwiseMult_Nest(Vec w,Vec x,Vec y)
221: {
222: Vec_Nest *bx = (Vec_Nest*)x->data;
223: Vec_Nest *by = (Vec_Nest*)y->data;
224: Vec_Nest *bw = (Vec_Nest*)w->data;
225: PetscInt i,nr;
229: VecNestCheckCompatible3(w,1,x,3,y,4);
230: nr = bx->nb;
231: for (i=0; i<nr; i++) {
232: VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
233: }
234: return(0);
235: }
237: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
238: {
239: Vec_Nest *bx = (Vec_Nest*)x->data;
240: Vec_Nest *by = (Vec_Nest*)y->data;
241: Vec_Nest *bw = (Vec_Nest*)w->data;
242: PetscInt i,nr;
246: VecNestCheckCompatible3(w,1,x,2,y,3);
248: nr = bx->nb;
249: for (i=0; i<nr; i++) {
250: VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
251: }
252: return(0);
253: }
255: static PetscErrorCode VecReciprocal_Nest(Vec x)
256: {
257: Vec_Nest *bx = (Vec_Nest*)x->data;
258: PetscInt i,nr;
262: nr = bx->nb;
263: for (i=0; i<nr; i++) {
264: VecReciprocal(bx->v[i]);
265: }
266: return(0);
267: }
269: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
270: {
271: Vec_Nest *bx = (Vec_Nest*)xin->data;
272: PetscInt i,nr;
273: PetscReal z_i;
274: PetscReal _z;
278: nr = bx->nb;
279: _z = 0.0;
281: if (type == NORM_2) {
282: PetscScalar dot;
283: VecDot(xin,xin,&dot);
284: _z = PetscAbsScalar(PetscSqrtScalar(dot));
285: } else if (type == NORM_1) {
286: for (i=0; i<nr; i++) {
287: VecNorm(bx->v[i],type,&z_i);
288: _z = _z + z_i;
289: }
290: } else if (type == NORM_INFINITY) {
291: for (i=0; i<nr; i++) {
292: VecNorm(bx->v[i],type,&z_i);
293: if (z_i > _z) _z = z_i;
294: }
295: }
297: *z = _z;
298: return(0);
299: }
301: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
302: {
303: PetscInt v;
307: for (v=0; v<nv; v++) {
308: /* Do axpy on each vector,v */
309: VecAXPY(y,alpha[v],x[v]);
310: }
311: return(0);
312: }
314: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
315: {
316: PetscInt j;
320: for (j=0; j<nv; j++) {
321: VecDot(x,y[j],&val[j]);
322: }
323: return(0);
324: }
326: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
327: {
328: PetscInt j;
332: for (j=0; j<nv; j++) {
333: VecTDot(x,y[j],&val[j]);
334: }
335: return(0);
336: }
338: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
339: {
340: Vec_Nest *bx = (Vec_Nest*)x->data;
341: PetscInt i,nr;
345: nr = bx->nb;
346: for (i=0; i<nr; i++) {
347: VecSet(bx->v[i],alpha);
348: }
349: return(0);
350: }
352: static PetscErrorCode VecConjugate_Nest(Vec x)
353: {
354: Vec_Nest *bx = (Vec_Nest*)x->data;
355: PetscInt j,nr;
359: nr = bx->nb;
360: for (j=0; j<nr; j++) {
361: VecConjugate(bx->v[j]);
362: }
363: return(0);
364: }
366: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
367: {
368: Vec_Nest *bx = (Vec_Nest*)x->data;
369: Vec_Nest *by = (Vec_Nest*)y->data;
370: PetscInt i,nr;
374: VecNestCheckCompatible2(x,1,y,2);
375: nr = bx->nb;
376: for (i=0; i<nr; i++) {
377: VecSwap(bx->v[i],by->v[i]);
378: }
379: return(0);
380: }
382: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
383: {
384: Vec_Nest *bx = (Vec_Nest*)x->data;
385: Vec_Nest *by = (Vec_Nest*)y->data;
386: Vec_Nest *bw = (Vec_Nest*)w->data;
387: PetscInt i,nr;
391: VecNestCheckCompatible3(w,1,x,3,y,4);
393: nr = bx->nb;
394: for (i=0; i<nr; i++) {
395: VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
396: }
397: return(0);
398: }
400: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
401: {
402: Vec_Nest *bx;
403: PetscInt i,nr;
404: PetscBool isnest;
405: PetscInt L;
406: PetscInt _entry_loc;
407: PetscReal _entry_val;
411: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
412: if (!isnest) {
413: /* Not nest */
414: VecMax(x,&_entry_loc,&_entry_val);
415: if (_entry_val > *max) {
416: *max = _entry_val;
417: *p = _entry_loc + *cnt;
418: }
419: VecGetSize(x,&L);
420: *cnt = *cnt + L;
421: return(0);
422: }
424: /* Otherwise we have a nest */
425: bx = (Vec_Nest*)x->data;
426: nr = bx->nb;
428: /* now descend recursively */
429: for (i=0; i<nr; i++) {
430: VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
431: }
432: return(0);
433: }
435: /* supports nested blocks */
436: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
437: {
438: PetscInt cnt;
442: cnt = 0;
443: *p = 0;
444: *max = PETSC_MIN_REAL;
445: VecMax_Nest_Recursive(x,&cnt,p,max);
446: return(0);
447: }
449: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
450: {
451: Vec_Nest *bx = (Vec_Nest*)x->data;
452: PetscInt i,nr,L,_entry_loc;
453: PetscBool isnest;
454: PetscReal _entry_val;
458: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
459: if (!isnest) {
460: /* Not nest */
461: VecMin(x,&_entry_loc,&_entry_val);
462: if (_entry_val < *min) {
463: *min = _entry_val;
464: *p = _entry_loc + *cnt;
465: }
466: VecGetSize(x,&L);
467: *cnt = *cnt + L;
468: return(0);
469: }
471: /* Otherwise we have a nest */
472: nr = bx->nb;
474: /* now descend recursively */
475: for (i=0; i<nr; i++) {
476: VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
477: }
478: return(0);
479: }
481: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
482: {
483: PetscInt cnt;
487: cnt = 0;
488: *p = 0;
489: *min = PETSC_MAX_REAL;
490: VecMin_Nest_Recursive(x,&cnt,p,min);
491: return(0);
492: }
494: /* supports nested blocks */
495: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
496: {
497: Vec_Nest *bx = (Vec_Nest*)x->data;
498: PetscBool isascii;
499: PetscInt i;
503: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
504: if (isascii) {
505: PetscViewerASCIIPushTab(viewer);
506: PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D, structure: \n",bx->nb);
507: for (i=0; i<bx->nb; i++) {
508: VecType type;
509: char name[256] = "",prefix[256] = "";
510: PetscInt NR;
512: VecGetSize(bx->v[i],&NR);
513: VecGetType(bx->v[i],&type);
514: if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
515: if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}
517: PetscViewerASCIIPrintf(viewer,"(%D) : %s%stype=%s, rows=%D \n",i,name,prefix,type,NR);
519: PetscViewerASCIIPushTab(viewer); /* push1 */
520: VecView(bx->v[i],viewer);
521: PetscViewerASCIIPopTab(viewer); /* pop1 */
522: }
523: PetscViewerASCIIPopTab(viewer); /* pop0 */
524: }
525: return(0);
526: }
528: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
529: {
530: Vec_Nest *bx;
531: PetscInt size,i,nr;
532: PetscBool isnest;
536: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
537: if (!isnest) {
538: /* Not nest */
539: if (globalsize) { VecGetSize(x,&size); }
540: else { VecGetLocalSize(x,&size); }
541: *L = *L + size;
542: return(0);
543: }
545: /* Otherwise we have a nest */
546: bx = (Vec_Nest*)x->data;
547: nr = bx->nb;
549: /* now descend recursively */
550: for (i=0; i<nr; i++) {
551: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
552: }
553: return(0);
554: }
556: /* Returns the sum of the global size of all the consituent vectors in the nest */
557: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
558: {
560: *N = x->map->N;
561: return(0);
562: }
564: /* Returns the sum of the local size of all the consituent vectors in the nest */
565: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
566: {
568: *n = x->map->n;
569: return(0);
570: }
572: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
573: {
574: Vec_Nest *bx = (Vec_Nest*)x->data;
575: Vec_Nest *by = (Vec_Nest*)y->data;
576: PetscInt i,nr;
577: PetscReal local_max,m;
581: VecNestCheckCompatible2(x,1,y,2);
582: nr = bx->nb;
583: m = 0.0;
584: for (i=0; i<nr; i++) {
585: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
586: if (local_max > m) m = local_max;
587: }
588: *max = m;
589: return(0);
590: }
592: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
593: {
594: Vec_Nest *bx = (Vec_Nest*)X->data;
595: PetscInt i;
599: *x = NULL;
600: for (i=0; i<bx->nb; i++) {
601: PetscBool issame = PETSC_FALSE;
602: ISEqual(is,bx->is[i],&issame);
603: if (issame) {
604: *x = bx->v[i];
605: PetscObjectReference((PetscObject)(*x));
606: break;
607: }
608: }
609: if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
610: return(0);
611: }
613: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
614: {
618: VecDestroy(x);
619: return(0);
620: }
622: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
623: {
624: Vec_Nest *bx = (Vec_Nest*)X->data;
625: PetscInt i,m,rstart,rend;
629: VecGetOwnershipRange(X,&rstart,&rend);
630: VecGetLocalSize(X,&m);
631: PetscMalloc1(m,x);
632: for (i=0; i<bx->nb; i++) {
633: Vec subvec = bx->v[i];
634: IS isy = bx->is[i];
635: PetscInt j,sm;
636: const PetscInt *ixy;
637: const PetscScalar *y;
638: VecGetLocalSize(subvec,&sm);
639: VecGetArrayRead(subvec,&y);
640: ISGetIndices(isy,&ixy);
641: for (j=0; j<sm; j++) {
642: PetscInt ix = ixy[j];
643: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
644: (*x)[ix-rstart] = y[j];
645: }
646: ISRestoreIndices(isy,&ixy);
647: VecRestoreArrayRead(subvec,&y);
648: }
649: return(0);
650: }
652: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
653: {
654: Vec_Nest *bx = (Vec_Nest*)X->data;
655: PetscInt i,m,rstart,rend;
659: VecGetOwnershipRange(X,&rstart,&rend);
660: VecGetLocalSize(X,&m);
661: for (i=0; i<bx->nb; i++) {
662: Vec subvec = bx->v[i];
663: IS isy = bx->is[i];
664: PetscInt j,sm;
665: const PetscInt *ixy;
666: PetscScalar *y;
667: VecGetLocalSize(subvec,&sm);
668: VecGetArray(subvec,&y);
669: ISGetIndices(isy,&ixy);
670: for (j=0; j<sm; j++) {
671: PetscInt ix = ixy[j];
672: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
673: y[j] = (*x)[ix-rstart];
674: }
675: ISRestoreIndices(isy,&ixy);
676: VecRestoreArray(subvec,&y);
677: }
678: PetscFree(*x);
679: return(0);
680: }
683: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
684: {
686: /* 0 */
687: ops->duplicate = VecDuplicate_Nest;
688: ops->duplicatevecs = VecDuplicateVecs_Default;
689: ops->destroyvecs = VecDestroyVecs_Default;
690: ops->dot = VecDot_Nest;
691: ops->mdot = VecMDot_Nest;
693: /* 5 */
694: ops->norm = VecNorm_Nest;
695: ops->tdot = VecTDot_Nest;
696: ops->mtdot = VecMTDot_Nest;
697: ops->scale = VecScale_Nest;
698: ops->copy = VecCopy_Nest;
700: /* 10 */
701: ops->set = VecSet_Nest;
702: ops->swap = VecSwap_Nest;
703: ops->axpy = VecAXPY_Nest;
704: ops->axpby = VecAXPBY_Nest;
705: ops->maxpy = VecMAXPY_Nest;
707: /* 15 */
708: ops->aypx = VecAYPX_Nest;
709: ops->waxpy = VecWAXPY_Nest;
710: ops->axpbypcz = 0;
711: ops->pointwisemult = VecPointwiseMult_Nest;
712: ops->pointwisedivide = VecPointwiseDivide_Nest;
713: /* 20 */
714: ops->setvalues = 0;
715: ops->assemblybegin = VecAssemblyBegin_Nest;
716: ops->assemblyend = VecAssemblyEnd_Nest;
717: ops->getarray = VecGetArray_Nest;
718: ops->getsize = VecGetSize_Nest;
720: /* 25 */
721: ops->getlocalsize = VecGetLocalSize_Nest;
722: ops->restorearray = VecRestoreArray_Nest;
723: ops->max = VecMax_Nest;
724: ops->min = VecMin_Nest;
725: ops->setrandom = 0;
727: /* 30 */
728: ops->setoption = 0;
729: ops->setvaluesblocked = 0;
730: ops->destroy = VecDestroy_Nest;
731: ops->view = VecView_Nest;
732: ops->placearray = 0;
734: /* 35 */
735: ops->replacearray = 0;
736: ops->dot_local = VecDot_Nest;
737: ops->tdot_local = VecTDot_Nest;
738: ops->norm_local = VecNorm_Nest;
739: ops->mdot_local = VecMDot_Nest;
741: /* 40 */
742: ops->mtdot_local = VecMTDot_Nest;
743: ops->load = 0;
744: ops->reciprocal = VecReciprocal_Nest;
745: ops->conjugate = VecConjugate_Nest;
746: ops->setlocaltoglobalmapping = 0;
748: /* 45 */
749: ops->setvalueslocal = 0;
750: ops->resetarray = 0;
751: ops->setfromoptions = 0;
752: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
753: ops->load = 0;
755: /* 50 */
756: ops->pointwisemax = 0;
757: ops->pointwisemaxabs = 0;
758: ops->pointwisemin = 0;
759: ops->getvalues = 0;
760: ops->sqrt = 0;
762: /* 55 */
763: ops->abs = 0;
764: ops->exp = 0;
765: ops->shift = 0;
766: ops->create = 0;
767: ops->stridegather = 0;
769: /* 60 */
770: ops->stridescatter = 0;
771: ops->dotnorm2 = VecDotNorm2_Nest;
772: ops->getsubvector = VecGetSubVector_Nest;
773: ops->restoresubvector = VecRestoreSubVector_Nest;
774: return(0);
775: }
777: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
778: {
779: Vec_Nest *b = (Vec_Nest*)x->data;
780: PetscInt i;
781: PetscInt row;
784: if (!m) return(0);
785: for (i=0; i<m; i++) {
786: row = idxm[i];
787: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
788: vec[i] = b->v[row];
789: }
790: return(0);
791: }
793: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
794: {
798: VecNestGetSubVecs_Private(X,1,&idxm,sx);
799: return(0);
800: }
802: /*@
803: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
805: Not collective
807: Input Parameters:
808: . X - nest vector
809: . idxm - index of the vector within the nest
811: Output Parameter:
812: . sx - vector at index idxm within the nest
814: Notes:
816: Level: developer
818: .seealso: VecNestGetSize(), VecNestGetSubVecs()
819: @*/
820: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
821: {
825: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
826: return(0);
827: }
829: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
830: {
831: Vec_Nest *b = (Vec_Nest*)X->data;
834: if (N) *N = b->nb;
835: if (sx) *sx = b->v;
836: return(0);
837: }
839: /*@C
840: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
842: Not collective
844: Input Parameters:
845: . X - nest vector
847: Output Parameter:
848: + N - number of nested vecs
849: - sx - array of vectors
851: Notes:
852: The user should not free the array sx.
854: Fortran Notes:
855: The caller must allocate the array to hold the subvectors.
857: Level: developer
859: .seealso: VecNestGetSize(), VecNestGetSubVec()
860: @*/
861: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
862: {
866: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
867: return(0);
868: }
870: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
871: {
872: Vec_Nest *bx = (Vec_Nest*)X->data;
873: PetscInt i,offset=0,n=0,bs;
874: IS is;
876: PetscBool issame = PETSC_FALSE;
877: PetscInt N=0;
879: /* check if idxm < bx->nb */
880: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
883: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
884: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
885: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
887: /* check if we need to update the IS for the block */
888: offset = X->map->rstart;
889: for (i=0; i<idxm; i++) {
890: n=0;
891: VecGetLocalSize(bx->v[i],&n);
892: offset += n;
893: }
895: /* get the local size and block size */
896: VecGetLocalSize(x,&n);
897: VecGetBlockSize(x,&bs);
899: /* create the new IS */
900: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
901: ISSetBlockSize(is,bs);
903: /* check if they are equal */
904: ISEqual(is,bx->is[idxm],&issame);
906: if (!issame) {
907: /* The IS of given vector has a different layout compared to the existing block vector.
908: Destroy the existing reference and update the IS. */
909: ISDestroy(&bx->is[idxm]);
910: ISDuplicate(is,&bx->is[idxm]);
911: ISCopy(is,bx->is[idxm]);
913: offset += n;
914: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
915: for (i=idxm+1; i<bx->nb; i++) {
916: /* get the local size and block size */
917: VecGetLocalSize(bx->v[i],&n);
918: VecGetBlockSize(bx->v[i],&bs);
920: /* destroy the old and create the new IS */
921: ISDestroy(&bx->is[i]);
922: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
923: ISSetBlockSize(bx->is[i],bs);
925: offset += n;
926: }
928: n=0;
929: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
930: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
931: PetscLayoutSetSize(X->map,N);
932: PetscLayoutSetLocalSize(X->map,n);
933: }
935: ISDestroy(&is);
936: return(0);
937: }
939: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
940: {
944: VecNestSetSubVec_Private(X,idxm,sx);
945: return(0);
946: }
948: /*@
949: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
951: Not collective
953: Input Parameters:
954: + X - nest vector
955: . idxm - index of the vector within the nest vector
956: - sx - vector at index idxm within the nest vector
958: Notes:
959: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
961: Level: developer
963: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
964: @*/
965: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
966: {
970: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
971: return(0);
972: }
974: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
975: {
976: PetscInt i;
980: for (i=0; i<N; i++) {
981: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
982: }
983: return(0);
984: }
986: /*@C
987: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
989: Not collective
991: Input Parameters:
992: + X - nest vector
993: . N - number of component vecs in sx
994: . idxm - indices of component vecs that are to be replaced
995: - sx - array of vectors
997: Notes:
998: The components in the vector array sx do not have to be of the same size as corresponding
999: components in X. The user can also free the array "sx" after the call.
1001: Level: developer
1003: .seealso: VecNestGetSize(), VecNestGetSubVec()
1004: @*/
1005: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1006: {
1010: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1011: return(0);
1012: }
1014: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1015: {
1016: Vec_Nest *b = (Vec_Nest*)X->data;
1019: *N = b->nb;
1020: return(0);
1021: }
1023: /*@
1024: VecNestGetSize - Returns the size of the nest vector.
1026: Not collective
1028: Input Parameters:
1029: . X - nest vector
1031: Output Parameter:
1032: . N - number of nested vecs
1034: Notes:
1036: Level: developer
1038: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1039: @*/
1040: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1041: {
1047: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1048: return(0);
1049: }
1051: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1052: {
1053: Vec_Nest *ctx = (Vec_Nest*)V->data;
1054: PetscInt i;
1058: if (ctx->setup_called) return(0);
1060: ctx->nb = nb;
1061: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1063: /* Create space */
1064: PetscMalloc1(ctx->nb,&ctx->v);
1065: for (i=0; i<ctx->nb; i++) {
1066: ctx->v[i] = x[i];
1067: PetscObjectReference((PetscObject)x[i]);
1068: /* Do not allocate memory for internal sub blocks */
1069: }
1071: PetscMalloc1(ctx->nb,&ctx->is);
1073: ctx->setup_called = PETSC_TRUE;
1074: return(0);
1075: }
1077: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1078: {
1079: Vec_Nest *ctx = (Vec_Nest*)V->data;
1080: PetscInt i,offset,m,n,M,N;
1084: if (is) { /* Do some consistency checks and reference the is */
1085: offset = V->map->rstart;
1086: for (i=0; i<ctx->nb; i++) {
1087: ISGetSize(is[i],&M);
1088: VecGetSize(ctx->v[i],&N);
1089: if (M != N) SETERRQ3(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"In slot %D, IS of size %D is not compatible with Vec of size %D",i,M,N);
1090: ISGetLocalSize(is[i],&m);
1091: VecGetLocalSize(ctx->v[i],&n);
1092: if (m != n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"In slot %D, IS of local size %D is not compatible with Vec of local size %D",i,m,n);
1093: #if defined(PETSC_USE_DEBUG)
1094: { /* This test can be expensive */
1095: PetscInt start;
1096: PetscBool contiguous;
1097: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1098: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1099: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1100: }
1101: #endif
1102: PetscObjectReference((PetscObject)is[i]);
1103: ctx->is[i] = is[i];
1104: offset += n;
1105: }
1106: } else { /* Create a contiguous ISStride for each entry */
1107: offset = V->map->rstart;
1108: for (i=0; i<ctx->nb; i++) {
1109: PetscInt bs;
1110: VecGetLocalSize(ctx->v[i],&n);
1111: VecGetBlockSize(ctx->v[i],&bs);
1112: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1113: ISSetBlockSize(ctx->is[i],bs);
1114: offset += n;
1115: }
1116: }
1117: return(0);
1118: }
1120: /*@C
1121: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1123: Collective on Vec
1125: Input Parameter:
1126: + comm - Communicator for the new Vec
1127: . nb - number of nested blocks
1128: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1129: - x - array of nb sub-vectors
1131: Output Parameter:
1132: . Y - new vector
1134: Level: advanced
1136: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1137: @*/
1138: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1139: {
1140: Vec V;
1141: Vec_Nest *s;
1142: PetscInt n,N;
1146: VecCreate(comm,&V);
1147: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1149: /* allocate and set pointer for implememtation data */
1150: PetscNew(&s);
1151: V->data = (void*)s;
1152: s->setup_called = PETSC_FALSE;
1153: s->nb = -1;
1154: s->v = NULL;
1156: VecSetUp_Nest_Private(V,nb,x);
1158: n = N = 0;
1159: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1160: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1161: PetscLayoutSetSize(V->map,N);
1162: PetscLayoutSetLocalSize(V->map,n);
1163: PetscLayoutSetBlockSize(V->map,1);
1164: PetscLayoutSetUp(V->map);
1166: VecSetUp_NestIS_Private(V,nb,is);
1168: VecNestSetOps_Private(V->ops);
1169: V->petscnative = PETSC_FALSE;
1172: /* expose block api's */
1173: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1174: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1175: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1176: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1177: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);
1179: *Y = V;
1180: return(0);
1181: }
1183: /*MC
1184: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1186: Level: intermediate
1188: Notes:
1189: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1190: It is usually used with MATNEST and DMComposite via DMSetVecType().
1192: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1193: M*/