Actual source code: vecnest.c
petsc-main 2021-04-15
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 VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
207: {
208: Vec_Nest *bx = (Vec_Nest*)x->data;
209: Vec_Nest *by = (Vec_Nest*)y->data;
210: Vec_Nest *bz = (Vec_Nest*)z->data;
211: PetscInt i,nr;
215: nr = bx->nb;
216: for (i=0; i<nr; i++) {
217: VecAXPBYPCZ(bz->v[i],alpha,beta,gamma,bx->v[i],by->v[i]);
218: }
219: return(0);
220: }
222: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
223: {
224: Vec_Nest *bx = (Vec_Nest*)x->data;
225: PetscInt i,nr;
229: nr = bx->nb;
230: for (i=0; i<nr; i++) {
231: VecScale(bx->v[i],alpha);
232: }
233: return(0);
234: }
236: static PetscErrorCode VecPointwiseMult_Nest(Vec w,Vec x,Vec y)
237: {
238: Vec_Nest *bx = (Vec_Nest*)x->data;
239: Vec_Nest *by = (Vec_Nest*)y->data;
240: Vec_Nest *bw = (Vec_Nest*)w->data;
241: PetscInt i,nr;
245: VecNestCheckCompatible3(w,1,x,3,y,4);
246: nr = bx->nb;
247: for (i=0; i<nr; i++) {
248: VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
249: }
250: return(0);
251: }
253: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
254: {
255: Vec_Nest *bx = (Vec_Nest*)x->data;
256: Vec_Nest *by = (Vec_Nest*)y->data;
257: Vec_Nest *bw = (Vec_Nest*)w->data;
258: PetscInt i,nr;
262: VecNestCheckCompatible3(w,1,x,2,y,3);
264: nr = bx->nb;
265: for (i=0; i<nr; i++) {
266: VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
267: }
268: return(0);
269: }
271: static PetscErrorCode VecReciprocal_Nest(Vec x)
272: {
273: Vec_Nest *bx = (Vec_Nest*)x->data;
274: PetscInt i,nr;
278: nr = bx->nb;
279: for (i=0; i<nr; i++) {
280: VecReciprocal(bx->v[i]);
281: }
282: return(0);
283: }
285: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
286: {
287: Vec_Nest *bx = (Vec_Nest*)xin->data;
288: PetscInt i,nr;
289: PetscReal z_i;
290: PetscReal _z;
294: nr = bx->nb;
295: _z = 0.0;
297: if (type == NORM_2) {
298: PetscScalar dot;
299: VecDot(xin,xin,&dot);
300: _z = PetscAbsScalar(PetscSqrtScalar(dot));
301: } else if (type == NORM_1) {
302: for (i=0; i<nr; i++) {
303: VecNorm(bx->v[i],type,&z_i);
304: _z = _z + z_i;
305: }
306: } else if (type == NORM_INFINITY) {
307: for (i=0; i<nr; i++) {
308: VecNorm(bx->v[i],type,&z_i);
309: if (z_i > _z) _z = z_i;
310: }
311: }
313: *z = _z;
314: return(0);
315: }
317: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
318: {
319: PetscInt v;
323: for (v=0; v<nv; v++) {
324: /* Do axpy on each vector,v */
325: VecAXPY(y,alpha[v],x[v]);
326: }
327: return(0);
328: }
330: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
331: {
332: PetscInt j;
336: for (j=0; j<nv; j++) {
337: VecDot(x,y[j],&val[j]);
338: }
339: return(0);
340: }
342: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
343: {
344: PetscInt j;
348: for (j=0; j<nv; j++) {
349: VecTDot(x,y[j],&val[j]);
350: }
351: return(0);
352: }
354: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
355: {
356: Vec_Nest *bx = (Vec_Nest*)x->data;
357: PetscInt i,nr;
361: nr = bx->nb;
362: for (i=0; i<nr; i++) {
363: VecSet(bx->v[i],alpha);
364: }
365: return(0);
366: }
368: static PetscErrorCode VecConjugate_Nest(Vec x)
369: {
370: Vec_Nest *bx = (Vec_Nest*)x->data;
371: PetscInt j,nr;
375: nr = bx->nb;
376: for (j=0; j<nr; j++) {
377: VecConjugate(bx->v[j]);
378: }
379: return(0);
380: }
382: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
383: {
384: Vec_Nest *bx = (Vec_Nest*)x->data;
385: Vec_Nest *by = (Vec_Nest*)y->data;
386: PetscInt i,nr;
390: VecNestCheckCompatible2(x,1,y,2);
391: nr = bx->nb;
392: for (i=0; i<nr; i++) {
393: VecSwap(bx->v[i],by->v[i]);
394: }
395: return(0);
396: }
398: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
399: {
400: Vec_Nest *bx = (Vec_Nest*)x->data;
401: Vec_Nest *by = (Vec_Nest*)y->data;
402: Vec_Nest *bw = (Vec_Nest*)w->data;
403: PetscInt i,nr;
407: VecNestCheckCompatible3(w,1,x,3,y,4);
409: nr = bx->nb;
410: for (i=0; i<nr; i++) {
411: VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
412: }
413: return(0);
414: }
416: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
417: {
418: Vec_Nest *bx;
419: PetscInt i,nr;
420: PetscBool isnest;
421: PetscInt L;
422: PetscInt _entry_loc;
423: PetscReal _entry_val;
427: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
428: if (!isnest) {
429: /* Not nest */
430: VecMax(x,&_entry_loc,&_entry_val);
431: if (_entry_val > *max) {
432: *max = _entry_val;
433: if (p) *p = _entry_loc + *cnt;
434: }
435: VecGetSize(x,&L);
436: *cnt = *cnt + L;
437: return(0);
438: }
440: /* Otherwise we have a nest */
441: bx = (Vec_Nest*)x->data;
442: nr = bx->nb;
444: /* now descend recursively */
445: for (i=0; i<nr; i++) {
446: VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
447: }
448: return(0);
449: }
451: /* supports nested blocks */
452: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
453: {
454: PetscInt cnt;
458: cnt = 0;
459: if (p) *p = 0;
460: *max = PETSC_MIN_REAL;
461: VecMax_Nest_Recursive(x,&cnt,p,max);
462: return(0);
463: }
465: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
466: {
467: Vec_Nest *bx = (Vec_Nest*)x->data;
468: PetscInt i,nr,L,_entry_loc;
469: PetscBool isnest;
470: PetscReal _entry_val;
474: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
475: if (!isnest) {
476: /* Not nest */
477: VecMin(x,&_entry_loc,&_entry_val);
478: if (_entry_val < *min) {
479: *min = _entry_val;
480: if (p) *p = _entry_loc + *cnt;
481: }
482: VecGetSize(x,&L);
483: *cnt = *cnt + L;
484: return(0);
485: }
487: /* Otherwise we have a nest */
488: nr = bx->nb;
490: /* now descend recursively */
491: for (i=0; i<nr; i++) {
492: VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
493: }
494: return(0);
495: }
497: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
498: {
499: PetscInt cnt;
503: cnt = 0;
504: if (p) *p = 0;
505: *min = PETSC_MAX_REAL;
506: VecMin_Nest_Recursive(x,&cnt,p,min);
507: return(0);
508: }
510: /* supports nested blocks */
511: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
512: {
513: Vec_Nest *bx = (Vec_Nest*)x->data;
514: PetscBool isascii;
515: PetscInt i;
519: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
520: if (isascii) {
521: PetscViewerASCIIPushTab(viewer);
522: PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D, structure: \n",bx->nb);
523: for (i=0; i<bx->nb; i++) {
524: VecType type;
525: char name[256] = "",prefix[256] = "";
526: PetscInt NR;
528: VecGetSize(bx->v[i],&NR);
529: VecGetType(bx->v[i],&type);
530: if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
531: if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}
533: PetscViewerASCIIPrintf(viewer,"(%D) : %s%stype=%s, rows=%D \n",i,name,prefix,type,NR);
535: PetscViewerASCIIPushTab(viewer); /* push1 */
536: VecView(bx->v[i],viewer);
537: PetscViewerASCIIPopTab(viewer); /* pop1 */
538: }
539: PetscViewerASCIIPopTab(viewer); /* pop0 */
540: }
541: return(0);
542: }
544: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
545: {
546: Vec_Nest *bx;
547: PetscInt size,i,nr;
548: PetscBool isnest;
552: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
553: if (!isnest) {
554: /* Not nest */
555: if (globalsize) { VecGetSize(x,&size); }
556: else { VecGetLocalSize(x,&size); }
557: *L = *L + size;
558: return(0);
559: }
561: /* Otherwise we have a nest */
562: bx = (Vec_Nest*)x->data;
563: nr = bx->nb;
565: /* now descend recursively */
566: for (i=0; i<nr; i++) {
567: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
568: }
569: return(0);
570: }
572: /* Returns the sum of the global size of all the consituent vectors in the nest */
573: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
574: {
576: *N = x->map->N;
577: return(0);
578: }
580: /* Returns the sum of the local size of all the consituent vectors in the nest */
581: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
582: {
584: *n = x->map->n;
585: return(0);
586: }
588: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
589: {
590: Vec_Nest *bx = (Vec_Nest*)x->data;
591: Vec_Nest *by = (Vec_Nest*)y->data;
592: PetscInt i,nr;
593: PetscReal local_max,m;
597: VecNestCheckCompatible2(x,1,y,2);
598: nr = bx->nb;
599: m = 0.0;
600: for (i=0; i<nr; i++) {
601: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
602: if (local_max > m) m = local_max;
603: }
604: *max = m;
605: return(0);
606: }
608: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
609: {
610: Vec_Nest *bx = (Vec_Nest*)X->data;
611: PetscInt i;
615: *x = NULL;
616: for (i=0; i<bx->nb; i++) {
617: PetscBool issame = PETSC_FALSE;
618: ISEqual(is,bx->is[i],&issame);
619: if (issame) {
620: *x = bx->v[i];
621: PetscObjectReference((PetscObject)(*x));
622: break;
623: }
624: }
625: if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
626: return(0);
627: }
629: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
630: {
634: VecDestroy(x);
635: return(0);
636: }
638: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
639: {
640: Vec_Nest *bx = (Vec_Nest*)X->data;
641: PetscInt i,m,rstart,rend;
645: VecGetOwnershipRange(X,&rstart,&rend);
646: VecGetLocalSize(X,&m);
647: PetscMalloc1(m,x);
648: for (i=0; i<bx->nb; i++) {
649: Vec subvec = bx->v[i];
650: IS isy = bx->is[i];
651: PetscInt j,sm;
652: const PetscInt *ixy;
653: const PetscScalar *y;
654: VecGetLocalSize(subvec,&sm);
655: VecGetArrayRead(subvec,&y);
656: ISGetIndices(isy,&ixy);
657: for (j=0; j<sm; j++) {
658: PetscInt ix = ixy[j];
659: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
660: (*x)[ix-rstart] = y[j];
661: }
662: ISRestoreIndices(isy,&ixy);
663: VecRestoreArrayRead(subvec,&y);
664: }
665: return(0);
666: }
668: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
669: {
670: Vec_Nest *bx = (Vec_Nest*)X->data;
671: PetscInt i,m,rstart,rend;
675: VecGetOwnershipRange(X,&rstart,&rend);
676: VecGetLocalSize(X,&m);
677: for (i=0; i<bx->nb; i++) {
678: Vec subvec = bx->v[i];
679: IS isy = bx->is[i];
680: PetscInt j,sm;
681: const PetscInt *ixy;
682: PetscScalar *y;
683: VecGetLocalSize(subvec,&sm);
684: VecGetArray(subvec,&y);
685: ISGetIndices(isy,&ixy);
686: for (j=0; j<sm; j++) {
687: PetscInt ix = ixy[j];
688: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
689: y[j] = (*x)[ix-rstart];
690: }
691: ISRestoreIndices(isy,&ixy);
692: VecRestoreArray(subvec,&y);
693: }
694: PetscFree(*x);
695: return(0);
696: }
698: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X,const PetscScalar **x)
699: {
703: PetscFree(*x);
704: return(0);
705: }
707: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
708: {
710: SETERRQ(PetscObjectComm((PetscObject)(*X)), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
711: return(0);
712: }
714: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
715: {
717: ops->duplicate = VecDuplicate_Nest;
718: ops->duplicatevecs = VecDuplicateVecs_Default;
719: ops->destroyvecs = VecDestroyVecs_Default;
720: ops->dot = VecDot_Nest;
721: ops->mdot = VecMDot_Nest;
722: ops->norm = VecNorm_Nest;
723: ops->tdot = VecTDot_Nest;
724: ops->mtdot = VecMTDot_Nest;
725: ops->scale = VecScale_Nest;
726: ops->copy = VecCopy_Nest;
727: ops->set = VecSet_Nest;
728: ops->swap = VecSwap_Nest;
729: ops->axpy = VecAXPY_Nest;
730: ops->axpby = VecAXPBY_Nest;
731: ops->maxpy = VecMAXPY_Nest;
732: ops->aypx = VecAYPX_Nest;
733: ops->waxpy = VecWAXPY_Nest;
734: ops->axpbypcz = NULL;
735: ops->pointwisemult = VecPointwiseMult_Nest;
736: ops->pointwisedivide = VecPointwiseDivide_Nest;
737: ops->setvalues = NULL;
738: ops->assemblybegin = VecAssemblyBegin_Nest;
739: ops->assemblyend = VecAssemblyEnd_Nest;
740: ops->getarray = VecGetArray_Nest;
741: ops->getsize = VecGetSize_Nest;
742: ops->getlocalsize = VecGetLocalSize_Nest;
743: ops->restorearray = VecRestoreArray_Nest;
744: ops->restorearrayread = VecRestoreArrayRead_Nest;
745: ops->max = VecMax_Nest;
746: ops->min = VecMin_Nest;
747: ops->setrandom = NULL;
748: ops->setoption = NULL;
749: ops->setvaluesblocked = NULL;
750: ops->destroy = VecDestroy_Nest;
751: ops->view = VecView_Nest;
752: ops->placearray = NULL;
753: ops->replacearray = NULL;
754: ops->dot_local = VecDot_Nest;
755: ops->tdot_local = VecTDot_Nest;
756: ops->norm_local = VecNorm_Nest;
757: ops->mdot_local = VecMDot_Nest;
758: ops->mtdot_local = VecMTDot_Nest;
759: ops->load = NULL;
760: ops->reciprocal = VecReciprocal_Nest;
761: ops->conjugate = VecConjugate_Nest;
762: ops->setlocaltoglobalmapping = NULL;
763: ops->setvalueslocal = NULL;
764: ops->resetarray = NULL;
765: ops->setfromoptions = NULL;
766: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
767: ops->load = NULL;
768: ops->pointwisemax = NULL;
769: ops->pointwisemaxabs = NULL;
770: ops->pointwisemin = NULL;
771: ops->getvalues = NULL;
772: ops->sqrt = NULL;
773: ops->abs = NULL;
774: ops->exp = NULL;
775: ops->shift = NULL;
776: ops->create = NULL;
777: ops->stridegather = NULL;
778: ops->stridescatter = NULL;
779: ops->dotnorm2 = VecDotNorm2_Nest;
780: ops->getsubvector = VecGetSubVector_Nest;
781: ops->restoresubvector = VecRestoreSubVector_Nest;
782: ops->axpbypcz = VecAXPBYPCZ_Nest;
783: ops->concatenate = VecConcatenate_Nest;
784: return(0);
785: }
787: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
788: {
789: Vec_Nest *b = (Vec_Nest*)x->data;
790: PetscInt i;
791: PetscInt row;
794: if (!m) return(0);
795: for (i=0; i<m; i++) {
796: row = idxm[i];
797: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
798: vec[i] = b->v[row];
799: }
800: return(0);
801: }
803: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
804: {
808: VecNestGetSubVecs_Private(X,1,&idxm,sx);
809: return(0);
810: }
812: /*@
813: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
815: Not collective
817: Input Parameters:
818: + X - nest vector
819: - idxm - index of the vector within the nest
821: Output Parameter:
822: . sx - vector at index idxm within the nest
824: Notes:
826: Level: developer
828: .seealso: VecNestGetSize(), VecNestGetSubVecs()
829: @*/
830: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
831: {
835: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
836: return(0);
837: }
839: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
840: {
841: Vec_Nest *b = (Vec_Nest*)X->data;
844: if (N) *N = b->nb;
845: if (sx) *sx = b->v;
846: return(0);
847: }
849: /*@C
850: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
852: Not collective
854: Input Parameter:
855: . X - nest vector
857: Output Parameters:
858: + N - number of nested vecs
859: - sx - array of vectors
861: Notes:
862: The user should not free the array sx.
864: Fortran Notes:
865: The caller must allocate the array to hold the subvectors.
867: Level: developer
869: .seealso: VecNestGetSize(), VecNestGetSubVec()
870: @*/
871: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
872: {
876: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
877: return(0);
878: }
880: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
881: {
882: Vec_Nest *bx = (Vec_Nest*)X->data;
883: PetscInt i,offset=0,n=0,bs;
884: IS is;
886: PetscBool issame = PETSC_FALSE;
887: PetscInt N=0;
889: /* check if idxm < bx->nb */
890: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
893: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
894: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
895: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
897: /* check if we need to update the IS for the block */
898: offset = X->map->rstart;
899: for (i=0; i<idxm; i++) {
900: n=0;
901: VecGetLocalSize(bx->v[i],&n);
902: offset += n;
903: }
905: /* get the local size and block size */
906: VecGetLocalSize(x,&n);
907: VecGetBlockSize(x,&bs);
909: /* create the new IS */
910: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
911: ISSetBlockSize(is,bs);
913: /* check if they are equal */
914: ISEqual(is,bx->is[idxm],&issame);
916: if (!issame) {
917: /* The IS of given vector has a different layout compared to the existing block vector.
918: Destroy the existing reference and update the IS. */
919: ISDestroy(&bx->is[idxm]);
920: ISDuplicate(is,&bx->is[idxm]);
921: ISCopy(is,bx->is[idxm]);
923: offset += n;
924: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
925: for (i=idxm+1; i<bx->nb; i++) {
926: /* get the local size and block size */
927: VecGetLocalSize(bx->v[i],&n);
928: VecGetBlockSize(bx->v[i],&bs);
930: /* destroy the old and create the new IS */
931: ISDestroy(&bx->is[i]);
932: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
933: ISSetBlockSize(bx->is[i],bs);
935: offset += n;
936: }
938: n=0;
939: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
940: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
941: PetscLayoutSetSize(X->map,N);
942: PetscLayoutSetLocalSize(X->map,n);
943: }
945: ISDestroy(&is);
946: return(0);
947: }
949: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
950: {
954: VecNestSetSubVec_Private(X,idxm,sx);
955: return(0);
956: }
958: /*@
959: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
961: Not collective
963: Input Parameters:
964: + X - nest vector
965: . idxm - index of the vector within the nest vector
966: - sx - vector at index idxm within the nest vector
968: Notes:
969: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
971: Level: developer
973: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
974: @*/
975: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
976: {
980: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
981: return(0);
982: }
984: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
985: {
986: PetscInt i;
990: for (i=0; i<N; i++) {
991: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
992: }
993: return(0);
994: }
996: /*@C
997: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
999: Not collective
1001: Input Parameters:
1002: + X - nest vector
1003: . N - number of component vecs in sx
1004: . idxm - indices of component vecs that are to be replaced
1005: - sx - array of vectors
1007: Notes:
1008: The components in the vector array sx do not have to be of the same size as corresponding
1009: components in X. The user can also free the array "sx" after the call.
1011: Level: developer
1013: .seealso: VecNestGetSize(), VecNestGetSubVec()
1014: @*/
1015: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1016: {
1020: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1021: return(0);
1022: }
1024: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1025: {
1026: Vec_Nest *b = (Vec_Nest*)X->data;
1029: *N = b->nb;
1030: return(0);
1031: }
1033: /*@
1034: VecNestGetSize - Returns the size of the nest vector.
1036: Not collective
1038: Input Parameter:
1039: . X - nest vector
1041: Output Parameter:
1042: . N - number of nested vecs
1044: Notes:
1046: Level: developer
1048: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1049: @*/
1050: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1051: {
1057: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1058: return(0);
1059: }
1061: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1062: {
1063: Vec_Nest *ctx = (Vec_Nest*)V->data;
1064: PetscInt i;
1068: if (ctx->setup_called) return(0);
1070: ctx->nb = nb;
1071: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1073: /* Create space */
1074: PetscMalloc1(ctx->nb,&ctx->v);
1075: for (i=0; i<ctx->nb; i++) {
1076: ctx->v[i] = x[i];
1077: PetscObjectReference((PetscObject)x[i]);
1078: /* Do not allocate memory for internal sub blocks */
1079: }
1081: PetscMalloc1(ctx->nb,&ctx->is);
1083: ctx->setup_called = PETSC_TRUE;
1084: return(0);
1085: }
1087: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1088: {
1089: Vec_Nest *ctx = (Vec_Nest*)V->data;
1090: PetscInt i,offset,m,n,M,N;
1094: if (is) { /* Do some consistency checks and reference the is */
1095: offset = V->map->rstart;
1096: for (i=0; i<ctx->nb; i++) {
1097: ISGetSize(is[i],&M);
1098: VecGetSize(ctx->v[i],&N);
1099: 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);
1100: ISGetLocalSize(is[i],&m);
1101: VecGetLocalSize(ctx->v[i],&n);
1102: 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);
1103: if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
1104: PetscInt start;
1105: PetscBool contiguous;
1106: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1107: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1108: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1109: }
1110: PetscObjectReference((PetscObject)is[i]);
1111: ctx->is[i] = is[i];
1112: offset += n;
1113: }
1114: } else { /* Create a contiguous ISStride for each entry */
1115: offset = V->map->rstart;
1116: for (i=0; i<ctx->nb; i++) {
1117: PetscInt bs;
1118: VecGetLocalSize(ctx->v[i],&n);
1119: VecGetBlockSize(ctx->v[i],&bs);
1120: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1121: ISSetBlockSize(ctx->is[i],bs);
1122: offset += n;
1123: }
1124: }
1125: return(0);
1126: }
1128: /*@C
1129: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1131: Collective on Vec
1133: Input Parameters:
1134: + comm - Communicator for the new Vec
1135: . nb - number of nested blocks
1136: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1137: - x - array of nb sub-vectors
1139: Output Parameter:
1140: . Y - new vector
1142: Level: advanced
1144: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1145: @*/
1146: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1147: {
1148: Vec V;
1149: Vec_Nest *s;
1150: PetscInt n,N;
1154: VecCreate(comm,&V);
1155: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1157: /* allocate and set pointer for implememtation data */
1158: PetscNew(&s);
1159: V->data = (void*)s;
1160: s->setup_called = PETSC_FALSE;
1161: s->nb = -1;
1162: s->v = NULL;
1164: VecSetUp_Nest_Private(V,nb,x);
1166: n = N = 0;
1167: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1168: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1169: PetscLayoutSetSize(V->map,N);
1170: PetscLayoutSetLocalSize(V->map,n);
1171: PetscLayoutSetBlockSize(V->map,1);
1172: PetscLayoutSetUp(V->map);
1174: VecSetUp_NestIS_Private(V,nb,is);
1176: VecNestSetOps_Private(V->ops);
1177: V->petscnative = PETSC_FALSE;
1180: /* expose block api's */
1181: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1182: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1183: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1184: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1185: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);
1187: *Y = V;
1188: return(0);
1189: }
1191: /*MC
1192: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1194: Level: intermediate
1196: Notes:
1197: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1198: It is usually used with MATNEST and DMComposite via DMSetVecType().
1200: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1201: M*/