Actual source code: vecnest.c
petsc-3.10.5 2019-03-28
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: *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: *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: *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: *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: }
699: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
700: {
702: /* 0 */
703: ops->duplicate = VecDuplicate_Nest;
704: ops->duplicatevecs = VecDuplicateVecs_Default;
705: ops->destroyvecs = VecDestroyVecs_Default;
706: ops->dot = VecDot_Nest;
707: ops->mdot = VecMDot_Nest;
709: /* 5 */
710: ops->norm = VecNorm_Nest;
711: ops->tdot = VecTDot_Nest;
712: ops->mtdot = VecMTDot_Nest;
713: ops->scale = VecScale_Nest;
714: ops->copy = VecCopy_Nest;
716: /* 10 */
717: ops->set = VecSet_Nest;
718: ops->swap = VecSwap_Nest;
719: ops->axpy = VecAXPY_Nest;
720: ops->axpby = VecAXPBY_Nest;
721: ops->maxpy = VecMAXPY_Nest;
723: /* 15 */
724: ops->aypx = VecAYPX_Nest;
725: ops->waxpy = VecWAXPY_Nest;
726: ops->axpbypcz = 0;
727: ops->pointwisemult = VecPointwiseMult_Nest;
728: ops->pointwisedivide = VecPointwiseDivide_Nest;
729: /* 20 */
730: ops->setvalues = 0;
731: ops->assemblybegin = VecAssemblyBegin_Nest;
732: ops->assemblyend = VecAssemblyEnd_Nest;
733: ops->getarray = VecGetArray_Nest;
734: ops->getsize = VecGetSize_Nest;
736: /* 25 */
737: ops->getlocalsize = VecGetLocalSize_Nest;
738: ops->restorearray = VecRestoreArray_Nest;
739: ops->max = VecMax_Nest;
740: ops->min = VecMin_Nest;
741: ops->setrandom = 0;
743: /* 30 */
744: ops->setoption = 0;
745: ops->setvaluesblocked = 0;
746: ops->destroy = VecDestroy_Nest;
747: ops->view = VecView_Nest;
748: ops->placearray = 0;
750: /* 35 */
751: ops->replacearray = 0;
752: ops->dot_local = VecDot_Nest;
753: ops->tdot_local = VecTDot_Nest;
754: ops->norm_local = VecNorm_Nest;
755: ops->mdot_local = VecMDot_Nest;
757: /* 40 */
758: ops->mtdot_local = VecMTDot_Nest;
759: ops->load = 0;
760: ops->reciprocal = VecReciprocal_Nest;
761: ops->conjugate = VecConjugate_Nest;
762: ops->setlocaltoglobalmapping = 0;
764: /* 45 */
765: ops->setvalueslocal = 0;
766: ops->resetarray = 0;
767: ops->setfromoptions = 0;
768: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
769: ops->load = 0;
771: /* 50 */
772: ops->pointwisemax = 0;
773: ops->pointwisemaxabs = 0;
774: ops->pointwisemin = 0;
775: ops->getvalues = 0;
776: ops->sqrt = 0;
778: /* 55 */
779: ops->abs = 0;
780: ops->exp = 0;
781: ops->shift = 0;
782: ops->create = 0;
783: ops->stridegather = 0;
785: /* 60 */
786: ops->stridescatter = 0;
787: ops->dotnorm2 = VecDotNorm2_Nest;
788: ops->getsubvector = VecGetSubVector_Nest;
789: ops->restoresubvector = VecRestoreSubVector_Nest;
790: ops->axpbypcz = VecAXPBYPCZ_Nest;
791: return(0);
792: }
794: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
795: {
796: Vec_Nest *b = (Vec_Nest*)x->data;
797: PetscInt i;
798: PetscInt row;
801: if (!m) return(0);
802: for (i=0; i<m; i++) {
803: row = idxm[i];
804: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
805: vec[i] = b->v[row];
806: }
807: return(0);
808: }
810: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
811: {
815: VecNestGetSubVecs_Private(X,1,&idxm,sx);
816: return(0);
817: }
819: /*@
820: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
822: Not collective
824: Input Parameters:
825: . X - nest vector
826: . idxm - index of the vector within the nest
828: Output Parameter:
829: . sx - vector at index idxm within the nest
831: Notes:
833: Level: developer
835: .seealso: VecNestGetSize(), VecNestGetSubVecs()
836: @*/
837: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
838: {
842: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
843: return(0);
844: }
846: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
847: {
848: Vec_Nest *b = (Vec_Nest*)X->data;
851: if (N) *N = b->nb;
852: if (sx) *sx = b->v;
853: return(0);
854: }
856: /*@C
857: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
859: Not collective
861: Input Parameters:
862: . X - nest vector
864: Output Parameter:
865: + N - number of nested vecs
866: - sx - array of vectors
868: Notes:
869: The user should not free the array sx.
871: Fortran Notes:
872: The caller must allocate the array to hold the subvectors.
874: Level: developer
876: .seealso: VecNestGetSize(), VecNestGetSubVec()
877: @*/
878: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
879: {
883: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
884: return(0);
885: }
887: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
888: {
889: Vec_Nest *bx = (Vec_Nest*)X->data;
890: PetscInt i,offset=0,n=0,bs;
891: IS is;
893: PetscBool issame = PETSC_FALSE;
894: PetscInt N=0;
896: /* check if idxm < bx->nb */
897: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
900: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
901: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
902: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
904: /* check if we need to update the IS for the block */
905: offset = X->map->rstart;
906: for (i=0; i<idxm; i++) {
907: n=0;
908: VecGetLocalSize(bx->v[i],&n);
909: offset += n;
910: }
912: /* get the local size and block size */
913: VecGetLocalSize(x,&n);
914: VecGetBlockSize(x,&bs);
916: /* create the new IS */
917: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
918: ISSetBlockSize(is,bs);
920: /* check if they are equal */
921: ISEqual(is,bx->is[idxm],&issame);
923: if (!issame) {
924: /* The IS of given vector has a different layout compared to the existing block vector.
925: Destroy the existing reference and update the IS. */
926: ISDestroy(&bx->is[idxm]);
927: ISDuplicate(is,&bx->is[idxm]);
928: ISCopy(is,bx->is[idxm]);
930: offset += n;
931: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
932: for (i=idxm+1; i<bx->nb; i++) {
933: /* get the local size and block size */
934: VecGetLocalSize(bx->v[i],&n);
935: VecGetBlockSize(bx->v[i],&bs);
937: /* destroy the old and create the new IS */
938: ISDestroy(&bx->is[i]);
939: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
940: ISSetBlockSize(bx->is[i],bs);
942: offset += n;
943: }
945: n=0;
946: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
947: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
948: PetscLayoutSetSize(X->map,N);
949: PetscLayoutSetLocalSize(X->map,n);
950: }
952: ISDestroy(&is);
953: return(0);
954: }
956: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
957: {
961: VecNestSetSubVec_Private(X,idxm,sx);
962: return(0);
963: }
965: /*@
966: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
968: Not collective
970: Input Parameters:
971: + X - nest vector
972: . idxm - index of the vector within the nest vector
973: - sx - vector at index idxm within the nest vector
975: Notes:
976: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
978: Level: developer
980: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
981: @*/
982: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
983: {
987: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
988: return(0);
989: }
991: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
992: {
993: PetscInt i;
997: for (i=0; i<N; i++) {
998: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
999: }
1000: return(0);
1001: }
1003: /*@C
1004: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1006: Not collective
1008: Input Parameters:
1009: + X - nest vector
1010: . N - number of component vecs in sx
1011: . idxm - indices of component vecs that are to be replaced
1012: - sx - array of vectors
1014: Notes:
1015: The components in the vector array sx do not have to be of the same size as corresponding
1016: components in X. The user can also free the array "sx" after the call.
1018: Level: developer
1020: .seealso: VecNestGetSize(), VecNestGetSubVec()
1021: @*/
1022: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1023: {
1027: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1028: return(0);
1029: }
1031: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1032: {
1033: Vec_Nest *b = (Vec_Nest*)X->data;
1036: *N = b->nb;
1037: return(0);
1038: }
1040: /*@
1041: VecNestGetSize - Returns the size of the nest vector.
1043: Not collective
1045: Input Parameters:
1046: . X - nest vector
1048: Output Parameter:
1049: . N - number of nested vecs
1051: Notes:
1053: Level: developer
1055: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1056: @*/
1057: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1058: {
1064: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1065: return(0);
1066: }
1068: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1069: {
1070: Vec_Nest *ctx = (Vec_Nest*)V->data;
1071: PetscInt i;
1075: if (ctx->setup_called) return(0);
1077: ctx->nb = nb;
1078: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1080: /* Create space */
1081: PetscMalloc1(ctx->nb,&ctx->v);
1082: for (i=0; i<ctx->nb; i++) {
1083: ctx->v[i] = x[i];
1084: PetscObjectReference((PetscObject)x[i]);
1085: /* Do not allocate memory for internal sub blocks */
1086: }
1088: PetscMalloc1(ctx->nb,&ctx->is);
1090: ctx->setup_called = PETSC_TRUE;
1091: return(0);
1092: }
1094: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1095: {
1096: Vec_Nest *ctx = (Vec_Nest*)V->data;
1097: PetscInt i,offset,m,n,M,N;
1101: if (is) { /* Do some consistency checks and reference the is */
1102: offset = V->map->rstart;
1103: for (i=0; i<ctx->nb; i++) {
1104: ISGetSize(is[i],&M);
1105: VecGetSize(ctx->v[i],&N);
1106: 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);
1107: ISGetLocalSize(is[i],&m);
1108: VecGetLocalSize(ctx->v[i],&n);
1109: 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);
1110: #if defined(PETSC_USE_DEBUG)
1111: { /* This test can be expensive */
1112: PetscInt start;
1113: PetscBool contiguous;
1114: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1115: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1116: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1117: }
1118: #endif
1119: PetscObjectReference((PetscObject)is[i]);
1120: ctx->is[i] = is[i];
1121: offset += n;
1122: }
1123: } else { /* Create a contiguous ISStride for each entry */
1124: offset = V->map->rstart;
1125: for (i=0; i<ctx->nb; i++) {
1126: PetscInt bs;
1127: VecGetLocalSize(ctx->v[i],&n);
1128: VecGetBlockSize(ctx->v[i],&bs);
1129: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1130: ISSetBlockSize(ctx->is[i],bs);
1131: offset += n;
1132: }
1133: }
1134: return(0);
1135: }
1137: /*@C
1138: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1140: Collective on Vec
1142: Input Parameter:
1143: + comm - Communicator for the new Vec
1144: . nb - number of nested blocks
1145: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1146: - x - array of nb sub-vectors
1148: Output Parameter:
1149: . Y - new vector
1151: Level: advanced
1153: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1154: @*/
1155: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1156: {
1157: Vec V;
1158: Vec_Nest *s;
1159: PetscInt n,N;
1163: VecCreate(comm,&V);
1164: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1166: /* allocate and set pointer for implememtation data */
1167: PetscNew(&s);
1168: V->data = (void*)s;
1169: s->setup_called = PETSC_FALSE;
1170: s->nb = -1;
1171: s->v = NULL;
1173: VecSetUp_Nest_Private(V,nb,x);
1175: n = N = 0;
1176: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1177: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1178: PetscLayoutSetSize(V->map,N);
1179: PetscLayoutSetLocalSize(V->map,n);
1180: PetscLayoutSetBlockSize(V->map,1);
1181: PetscLayoutSetUp(V->map);
1183: VecSetUp_NestIS_Private(V,nb,is);
1185: VecNestSetOps_Private(V->ops);
1186: V->petscnative = PETSC_FALSE;
1189: /* expose block api's */
1190: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1191: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1192: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1193: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1194: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);
1196: *Y = V;
1197: return(0);
1198: }
1200: /*MC
1201: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1203: Level: intermediate
1205: Notes:
1206: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1207: It is usually used with MATNEST and DMComposite via DMSetVecType().
1209: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1210: M*/