Actual source code: vecnest.c
petsc-master 2020-08-09
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 VecNestSetOps_Private(struct _VecOps *ops)
708: {
710: ops->duplicate = VecDuplicate_Nest;
711: ops->duplicatevecs = VecDuplicateVecs_Default;
712: ops->destroyvecs = VecDestroyVecs_Default;
713: ops->dot = VecDot_Nest;
714: ops->mdot = VecMDot_Nest;
715: ops->norm = VecNorm_Nest;
716: ops->tdot = VecTDot_Nest;
717: ops->mtdot = VecMTDot_Nest;
718: ops->scale = VecScale_Nest;
719: ops->copy = VecCopy_Nest;
720: ops->set = VecSet_Nest;
721: ops->swap = VecSwap_Nest;
722: ops->axpy = VecAXPY_Nest;
723: ops->axpby = VecAXPBY_Nest;
724: ops->maxpy = VecMAXPY_Nest;
725: ops->aypx = VecAYPX_Nest;
726: ops->waxpy = VecWAXPY_Nest;
727: ops->axpbypcz = 0;
728: ops->pointwisemult = VecPointwiseMult_Nest;
729: ops->pointwisedivide = VecPointwiseDivide_Nest;
730: ops->setvalues = 0;
731: ops->assemblybegin = VecAssemblyBegin_Nest;
732: ops->assemblyend = VecAssemblyEnd_Nest;
733: ops->getarray = VecGetArray_Nest;
734: ops->getsize = VecGetSize_Nest;
735: ops->getlocalsize = VecGetLocalSize_Nest;
736: ops->restorearray = VecRestoreArray_Nest;
737: ops->restorearrayread = VecRestoreArrayRead_Nest;
738: ops->max = VecMax_Nest;
739: ops->min = VecMin_Nest;
740: ops->setrandom = 0;
741: ops->setoption = 0;
742: ops->setvaluesblocked = 0;
743: ops->destroy = VecDestroy_Nest;
744: ops->view = VecView_Nest;
745: ops->placearray = 0;
746: ops->replacearray = 0;
747: ops->dot_local = VecDot_Nest;
748: ops->tdot_local = VecTDot_Nest;
749: ops->norm_local = VecNorm_Nest;
750: ops->mdot_local = VecMDot_Nest;
751: ops->mtdot_local = VecMTDot_Nest;
752: ops->load = 0;
753: ops->reciprocal = VecReciprocal_Nest;
754: ops->conjugate = VecConjugate_Nest;
755: ops->setlocaltoglobalmapping = 0;
756: ops->setvalueslocal = 0;
757: ops->resetarray = 0;
758: ops->setfromoptions = 0;
759: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
760: ops->load = 0;
761: ops->pointwisemax = 0;
762: ops->pointwisemaxabs = 0;
763: ops->pointwisemin = 0;
764: ops->getvalues = 0;
765: ops->sqrt = 0;
766: ops->abs = 0;
767: ops->exp = 0;
768: ops->shift = 0;
769: ops->create = 0;
770: ops->stridegather = 0;
771: ops->stridescatter = 0;
772: ops->dotnorm2 = VecDotNorm2_Nest;
773: ops->getsubvector = VecGetSubVector_Nest;
774: ops->restoresubvector = VecRestoreSubVector_Nest;
775: ops->axpbypcz = VecAXPBYPCZ_Nest;
776: return(0);
777: }
779: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
780: {
781: Vec_Nest *b = (Vec_Nest*)x->data;
782: PetscInt i;
783: PetscInt row;
786: if (!m) return(0);
787: for (i=0; i<m; i++) {
788: row = idxm[i];
789: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
790: vec[i] = b->v[row];
791: }
792: return(0);
793: }
795: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
796: {
800: VecNestGetSubVecs_Private(X,1,&idxm,sx);
801: return(0);
802: }
804: /*@
805: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
807: Not collective
809: Input Parameters:
810: . X - nest vector
811: . idxm - index of the vector within the nest
813: Output Parameter:
814: . sx - vector at index idxm within the nest
816: Notes:
818: Level: developer
820: .seealso: VecNestGetSize(), VecNestGetSubVecs()
821: @*/
822: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
823: {
827: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
828: return(0);
829: }
831: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
832: {
833: Vec_Nest *b = (Vec_Nest*)X->data;
836: if (N) *N = b->nb;
837: if (sx) *sx = b->v;
838: return(0);
839: }
841: /*@C
842: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
844: Not collective
846: Input Parameters:
847: . X - nest vector
849: Output Parameter:
850: + N - number of nested vecs
851: - sx - array of vectors
853: Notes:
854: The user should not free the array sx.
856: Fortran Notes:
857: The caller must allocate the array to hold the subvectors.
859: Level: developer
861: .seealso: VecNestGetSize(), VecNestGetSubVec()
862: @*/
863: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
864: {
868: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
869: return(0);
870: }
872: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
873: {
874: Vec_Nest *bx = (Vec_Nest*)X->data;
875: PetscInt i,offset=0,n=0,bs;
876: IS is;
878: PetscBool issame = PETSC_FALSE;
879: PetscInt N=0;
881: /* check if idxm < bx->nb */
882: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
885: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
886: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
887: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
889: /* check if we need to update the IS for the block */
890: offset = X->map->rstart;
891: for (i=0; i<idxm; i++) {
892: n=0;
893: VecGetLocalSize(bx->v[i],&n);
894: offset += n;
895: }
897: /* get the local size and block size */
898: VecGetLocalSize(x,&n);
899: VecGetBlockSize(x,&bs);
901: /* create the new IS */
902: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
903: ISSetBlockSize(is,bs);
905: /* check if they are equal */
906: ISEqual(is,bx->is[idxm],&issame);
908: if (!issame) {
909: /* The IS of given vector has a different layout compared to the existing block vector.
910: Destroy the existing reference and update the IS. */
911: ISDestroy(&bx->is[idxm]);
912: ISDuplicate(is,&bx->is[idxm]);
913: ISCopy(is,bx->is[idxm]);
915: offset += n;
916: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
917: for (i=idxm+1; i<bx->nb; i++) {
918: /* get the local size and block size */
919: VecGetLocalSize(bx->v[i],&n);
920: VecGetBlockSize(bx->v[i],&bs);
922: /* destroy the old and create the new IS */
923: ISDestroy(&bx->is[i]);
924: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
925: ISSetBlockSize(bx->is[i],bs);
927: offset += n;
928: }
930: n=0;
931: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
932: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
933: PetscLayoutSetSize(X->map,N);
934: PetscLayoutSetLocalSize(X->map,n);
935: }
937: ISDestroy(&is);
938: return(0);
939: }
941: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
942: {
946: VecNestSetSubVec_Private(X,idxm,sx);
947: return(0);
948: }
950: /*@
951: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
953: Not collective
955: Input Parameters:
956: + X - nest vector
957: . idxm - index of the vector within the nest vector
958: - sx - vector at index idxm within the nest vector
960: Notes:
961: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
963: Level: developer
965: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
966: @*/
967: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
968: {
972: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
973: return(0);
974: }
976: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
977: {
978: PetscInt i;
982: for (i=0; i<N; i++) {
983: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
984: }
985: return(0);
986: }
988: /*@C
989: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
991: Not collective
993: Input Parameters:
994: + X - nest vector
995: . N - number of component vecs in sx
996: . idxm - indices of component vecs that are to be replaced
997: - sx - array of vectors
999: Notes:
1000: The components in the vector array sx do not have to be of the same size as corresponding
1001: components in X. The user can also free the array "sx" after the call.
1003: Level: developer
1005: .seealso: VecNestGetSize(), VecNestGetSubVec()
1006: @*/
1007: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1008: {
1012: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1013: return(0);
1014: }
1016: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1017: {
1018: Vec_Nest *b = (Vec_Nest*)X->data;
1021: *N = b->nb;
1022: return(0);
1023: }
1025: /*@
1026: VecNestGetSize - Returns the size of the nest vector.
1028: Not collective
1030: Input Parameters:
1031: . X - nest vector
1033: Output Parameter:
1034: . N - number of nested vecs
1036: Notes:
1038: Level: developer
1040: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1041: @*/
1042: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1043: {
1049: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1050: return(0);
1051: }
1053: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1054: {
1055: Vec_Nest *ctx = (Vec_Nest*)V->data;
1056: PetscInt i;
1060: if (ctx->setup_called) return(0);
1062: ctx->nb = nb;
1063: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1065: /* Create space */
1066: PetscMalloc1(ctx->nb,&ctx->v);
1067: for (i=0; i<ctx->nb; i++) {
1068: ctx->v[i] = x[i];
1069: PetscObjectReference((PetscObject)x[i]);
1070: /* Do not allocate memory for internal sub blocks */
1071: }
1073: PetscMalloc1(ctx->nb,&ctx->is);
1075: ctx->setup_called = PETSC_TRUE;
1076: return(0);
1077: }
1079: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1080: {
1081: Vec_Nest *ctx = (Vec_Nest*)V->data;
1082: PetscInt i,offset,m,n,M,N;
1086: if (is) { /* Do some consistency checks and reference the is */
1087: offset = V->map->rstart;
1088: for (i=0; i<ctx->nb; i++) {
1089: ISGetSize(is[i],&M);
1090: VecGetSize(ctx->v[i],&N);
1091: 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);
1092: ISGetLocalSize(is[i],&m);
1093: VecGetLocalSize(ctx->v[i],&n);
1094: 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);
1095: if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
1096: PetscInt start;
1097: PetscBool contiguous;
1098: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1099: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1100: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1101: }
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*/