Actual source code: vecnest.c
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;
69: VecNestCheckCompatible2(x,1,y,2);
70: for (i=0; i<bx->nb; i++) {
71: VecCopy(bx->v[i],by->v[i]);
72: }
73: return(0);
74: }
76: /* supports nested blocks */
77: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
78: {
79: Vec_Nest *bx = (Vec_Nest*)x->data;
80: Vec Y;
81: Vec *sub;
82: PetscInt i;
86: PetscMalloc1(bx->nb,&sub);
87: for (i=0; i<bx->nb; i++) {
88: VecDuplicate(bx->v[i],&sub[i]);
89: }
90: VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
91: for (i=0; i<bx->nb; i++) {
92: VecDestroy(&sub[i]);
93: }
94: PetscFree(sub);
95: *y = Y;
96: return(0);
97: }
99: /* supports nested blocks */
100: static PetscErrorCode VecDot_Nest(Vec x,Vec y,PetscScalar *val)
101: {
102: Vec_Nest *bx = (Vec_Nest*)x->data;
103: Vec_Nest *by = (Vec_Nest*)y->data;
104: PetscInt i,nr;
105: PetscScalar x_dot_y,_val;
109: nr = bx->nb;
110: _val = 0.0;
111: for (i=0; i<nr; i++) {
112: VecDot(bx->v[i],by->v[i],&x_dot_y);
113: _val = _val + x_dot_y;
114: }
115: *val = _val;
116: return(0);
117: }
119: /* supports nested blocks */
120: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
121: {
122: Vec_Nest *bx = (Vec_Nest*)x->data;
123: Vec_Nest *by = (Vec_Nest*)y->data;
124: PetscInt i,nr;
125: PetscScalar x_dot_y,_val;
129: nr = bx->nb;
130: _val = 0.0;
131: for (i=0; i<nr; i++) {
132: VecTDot(bx->v[i],by->v[i],&x_dot_y);
133: _val = _val + x_dot_y;
134: }
135: *val = _val;
136: return(0);
137: }
139: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
140: {
141: Vec_Nest *bx = (Vec_Nest*)x->data;
142: Vec_Nest *by = (Vec_Nest*)y->data;
143: PetscInt i,nr;
144: PetscScalar x_dot_y,_dp,_nm;
145: PetscReal norm2_y;
149: nr = bx->nb;
150: _dp = 0.0;
151: _nm = 0.0;
152: for (i=0; i<nr; i++) {
153: VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
154: _dp += x_dot_y;
155: _nm += norm2_y;
156: }
157: *dp = _dp;
158: *nm = _nm;
159: return(0);
160: }
162: static PetscErrorCode VecAXPY_Nest(Vec y,PetscScalar alpha,Vec x)
163: {
164: Vec_Nest *bx = (Vec_Nest*)x->data;
165: Vec_Nest *by = (Vec_Nest*)y->data;
166: PetscInt i,nr;
170: nr = bx->nb;
171: for (i=0; i<nr; i++) {
172: VecAXPY(by->v[i],alpha,bx->v[i]);
173: }
174: return(0);
175: }
177: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
178: {
179: Vec_Nest *bx = (Vec_Nest*)x->data;
180: Vec_Nest *by = (Vec_Nest*)y->data;
181: PetscInt i,nr;
185: nr = bx->nb;
186: for (i=0; i<nr; i++) {
187: VecAYPX(by->v[i],alpha,bx->v[i]);
188: }
189: return(0);
190: }
192: static PetscErrorCode VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
193: {
194: Vec_Nest *bx = (Vec_Nest*)x->data;
195: Vec_Nest *by = (Vec_Nest*)y->data;
196: PetscInt i,nr;
200: nr = bx->nb;
201: for (i=0; i<nr; i++) {
202: VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
203: }
204: return(0);
205: }
207: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
208: {
209: Vec_Nest *bx = (Vec_Nest*)x->data;
210: Vec_Nest *by = (Vec_Nest*)y->data;
211: Vec_Nest *bz = (Vec_Nest*)z->data;
212: PetscInt i,nr;
216: nr = bx->nb;
217: for (i=0; i<nr; i++) {
218: VecAXPBYPCZ(bz->v[i],alpha,beta,gamma,bx->v[i],by->v[i]);
219: }
220: return(0);
221: }
223: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
224: {
225: Vec_Nest *bx = (Vec_Nest*)x->data;
226: PetscInt i,nr;
230: nr = bx->nb;
231: for (i=0; i<nr; i++) {
232: VecScale(bx->v[i],alpha);
233: }
234: return(0);
235: }
237: static PetscErrorCode VecPointwiseMult_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);
247: nr = bx->nb;
248: for (i=0; i<nr; i++) {
249: VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
250: }
251: return(0);
252: }
254: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
255: {
256: Vec_Nest *bx = (Vec_Nest*)x->data;
257: Vec_Nest *by = (Vec_Nest*)y->data;
258: Vec_Nest *bw = (Vec_Nest*)w->data;
259: PetscInt i,nr;
263: VecNestCheckCompatible3(w,1,x,2,y,3);
265: nr = bx->nb;
266: for (i=0; i<nr; i++) {
267: VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
268: }
269: return(0);
270: }
272: static PetscErrorCode VecReciprocal_Nest(Vec x)
273: {
274: Vec_Nest *bx = (Vec_Nest*)x->data;
275: PetscInt i,nr;
279: nr = bx->nb;
280: for (i=0; i<nr; i++) {
281: VecReciprocal(bx->v[i]);
282: }
283: return(0);
284: }
286: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
287: {
288: Vec_Nest *bx = (Vec_Nest*)xin->data;
289: PetscInt i,nr;
290: PetscReal z_i;
291: PetscReal _z;
295: nr = bx->nb;
296: _z = 0.0;
298: if (type == NORM_2) {
299: PetscScalar dot;
300: VecDot(xin,xin,&dot);
301: _z = PetscAbsScalar(PetscSqrtScalar(dot));
302: } else if (type == NORM_1) {
303: for (i=0; i<nr; i++) {
304: VecNorm(bx->v[i],type,&z_i);
305: _z = _z + z_i;
306: }
307: } else if (type == NORM_INFINITY) {
308: for (i=0; i<nr; i++) {
309: VecNorm(bx->v[i],type,&z_i);
310: if (z_i > _z) _z = z_i;
311: }
312: }
314: *z = _z;
315: return(0);
316: }
318: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
319: {
320: PetscInt v;
324: for (v=0; v<nv; v++) {
325: /* Do axpy on each vector,v */
326: VecAXPY(y,alpha[v],x[v]);
327: }
328: return(0);
329: }
331: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
332: {
333: PetscInt j;
337: for (j=0; j<nv; j++) {
338: VecDot(x,y[j],&val[j]);
339: }
340: return(0);
341: }
343: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
344: {
345: PetscInt j;
349: for (j=0; j<nv; j++) {
350: VecTDot(x,y[j],&val[j]);
351: }
352: return(0);
353: }
355: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
356: {
357: Vec_Nest *bx = (Vec_Nest*)x->data;
358: PetscInt i,nr;
362: nr = bx->nb;
363: for (i=0; i<nr; i++) {
364: VecSet(bx->v[i],alpha);
365: }
366: return(0);
367: }
369: static PetscErrorCode VecConjugate_Nest(Vec x)
370: {
371: Vec_Nest *bx = (Vec_Nest*)x->data;
372: PetscInt j,nr;
376: nr = bx->nb;
377: for (j=0; j<nr; j++) {
378: VecConjugate(bx->v[j]);
379: }
380: return(0);
381: }
383: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
384: {
385: Vec_Nest *bx = (Vec_Nest*)x->data;
386: Vec_Nest *by = (Vec_Nest*)y->data;
387: PetscInt i,nr;
391: VecNestCheckCompatible2(x,1,y,2);
392: nr = bx->nb;
393: for (i=0; i<nr; i++) {
394: VecSwap(bx->v[i],by->v[i]);
395: }
396: return(0);
397: }
399: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
400: {
401: Vec_Nest *bx = (Vec_Nest*)x->data;
402: Vec_Nest *by = (Vec_Nest*)y->data;
403: Vec_Nest *bw = (Vec_Nest*)w->data;
404: PetscInt i,nr;
408: VecNestCheckCompatible3(w,1,x,3,y,4);
410: nr = bx->nb;
411: for (i=0; i<nr; i++) {
412: VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
413: }
414: return(0);
415: }
417: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
418: {
419: Vec_Nest *bx;
420: PetscInt i,nr;
421: PetscBool isnest;
422: PetscInt L;
423: PetscInt _entry_loc;
424: PetscReal _entry_val;
428: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
429: if (!isnest) {
430: /* Not nest */
431: VecMax(x,&_entry_loc,&_entry_val);
432: if (_entry_val > *max) {
433: *max = _entry_val;
434: if (p) *p = _entry_loc + *cnt;
435: }
436: VecGetSize(x,&L);
437: *cnt = *cnt + L;
438: return(0);
439: }
441: /* Otherwise we have a nest */
442: bx = (Vec_Nest*)x->data;
443: nr = bx->nb;
445: /* now descend recursively */
446: for (i=0; i<nr; i++) {
447: VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
448: }
449: return(0);
450: }
452: /* supports nested blocks */
453: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
454: {
455: PetscInt cnt;
459: cnt = 0;
460: if (p) *p = 0;
461: *max = PETSC_MIN_REAL;
462: VecMax_Nest_Recursive(x,&cnt,p,max);
463: return(0);
464: }
466: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
467: {
468: Vec_Nest *bx = (Vec_Nest*)x->data;
469: PetscInt i,nr,L,_entry_loc;
470: PetscBool isnest;
471: PetscReal _entry_val;
475: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
476: if (!isnest) {
477: /* Not nest */
478: VecMin(x,&_entry_loc,&_entry_val);
479: if (_entry_val < *min) {
480: *min = _entry_val;
481: if (p) *p = _entry_loc + *cnt;
482: }
483: VecGetSize(x,&L);
484: *cnt = *cnt + L;
485: return(0);
486: }
488: /* Otherwise we have a nest */
489: nr = bx->nb;
491: /* now descend recursively */
492: for (i=0; i<nr; i++) {
493: VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
494: }
495: return(0);
496: }
498: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
499: {
500: PetscInt cnt;
504: cnt = 0;
505: if (p) *p = 0;
506: *min = PETSC_MAX_REAL;
507: VecMin_Nest_Recursive(x,&cnt,p,min);
508: return(0);
509: }
511: /* supports nested blocks */
512: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
513: {
514: Vec_Nest *bx = (Vec_Nest*)x->data;
515: PetscBool isascii;
516: PetscInt i;
520: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
521: if (isascii) {
522: PetscViewerASCIIPushTab(viewer);
523: PetscViewerASCIIPrintf(viewer,"VecNest, rows=%D, structure: \n",bx->nb);
524: for (i=0; i<bx->nb; i++) {
525: VecType type;
526: char name[256] = "",prefix[256] = "";
527: PetscInt NR;
529: VecGetSize(bx->v[i],&NR);
530: VecGetType(bx->v[i],&type);
531: if (((PetscObject)bx->v[i])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);}
532: if (((PetscObject)bx->v[i])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);}
534: PetscViewerASCIIPrintf(viewer,"(%D) : %s%stype=%s, rows=%D \n",i,name,prefix,type,NR);
536: PetscViewerASCIIPushTab(viewer); /* push1 */
537: VecView(bx->v[i],viewer);
538: PetscViewerASCIIPopTab(viewer); /* pop1 */
539: }
540: PetscViewerASCIIPopTab(viewer); /* pop0 */
541: }
542: return(0);
543: }
545: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
546: {
547: Vec_Nest *bx;
548: PetscInt size,i,nr;
549: PetscBool isnest;
553: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
554: if (!isnest) {
555: /* Not nest */
556: if (globalsize) { VecGetSize(x,&size); }
557: else { VecGetLocalSize(x,&size); }
558: *L = *L + size;
559: return(0);
560: }
562: /* Otherwise we have a nest */
563: bx = (Vec_Nest*)x->data;
564: nr = bx->nb;
566: /* now descend recursively */
567: for (i=0; i<nr; i++) {
568: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
569: }
570: return(0);
571: }
573: /* Returns the sum of the global size of all the consituent vectors in the nest */
574: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
575: {
577: *N = x->map->N;
578: return(0);
579: }
581: /* Returns the sum of the local size of all the consituent vectors in the nest */
582: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
583: {
585: *n = x->map->n;
586: return(0);
587: }
589: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
590: {
591: Vec_Nest *bx = (Vec_Nest*)x->data;
592: Vec_Nest *by = (Vec_Nest*)y->data;
593: PetscInt i,nr;
594: PetscReal local_max,m;
598: VecNestCheckCompatible2(x,1,y,2);
599: nr = bx->nb;
600: m = 0.0;
601: for (i=0; i<nr; i++) {
602: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
603: if (local_max > m) m = local_max;
604: }
605: *max = m;
606: return(0);
607: }
609: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
610: {
611: Vec_Nest *bx = (Vec_Nest*)X->data;
612: PetscInt i;
616: *x = NULL;
617: for (i=0; i<bx->nb; i++) {
618: PetscBool issame = PETSC_FALSE;
619: ISEqual(is,bx->is[i],&issame);
620: if (issame) {
621: *x = bx->v[i];
622: PetscObjectReference((PetscObject)(*x));
623: break;
624: }
625: }
626: if (!*x) SETERRQ(PetscObjectComm((PetscObject)is),PETSC_ERR_ARG_OUTOFRANGE,"Index set not found in nested Vec");
627: return(0);
628: }
630: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
631: {
635: VecDestroy(x);
636: return(0);
637: }
639: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
640: {
641: Vec_Nest *bx = (Vec_Nest*)X->data;
642: PetscInt i,m,rstart,rend;
646: VecGetOwnershipRange(X,&rstart,&rend);
647: VecGetLocalSize(X,&m);
648: PetscMalloc1(m,x);
649: for (i=0; i<bx->nb; i++) {
650: Vec subvec = bx->v[i];
651: IS isy = bx->is[i];
652: PetscInt j,sm;
653: const PetscInt *ixy;
654: const PetscScalar *y;
655: VecGetLocalSize(subvec,&sm);
656: VecGetArrayRead(subvec,&y);
657: ISGetIndices(isy,&ixy);
658: for (j=0; j<sm; j++) {
659: PetscInt ix = ixy[j];
660: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
661: (*x)[ix-rstart] = y[j];
662: }
663: ISRestoreIndices(isy,&ixy);
664: VecRestoreArrayRead(subvec,&y);
665: }
666: return(0);
667: }
669: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
670: {
671: Vec_Nest *bx = (Vec_Nest*)X->data;
672: PetscInt i,m,rstart,rend;
676: VecGetOwnershipRange(X,&rstart,&rend);
677: VecGetLocalSize(X,&m);
678: for (i=0; i<bx->nb; i++) {
679: Vec subvec = bx->v[i];
680: IS isy = bx->is[i];
681: PetscInt j,sm;
682: const PetscInt *ixy;
683: PetscScalar *y;
684: VecGetLocalSize(subvec,&sm);
685: VecGetArray(subvec,&y);
686: ISGetIndices(isy,&ixy);
687: for (j=0; j<sm; j++) {
688: PetscInt ix = ixy[j];
689: if (ix < rstart || rend <= ix) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for getting array from nonlocal subvec");
690: y[j] = (*x)[ix-rstart];
691: }
692: ISRestoreIndices(isy,&ixy);
693: VecRestoreArray(subvec,&y);
694: }
695: PetscFree(*x);
696: return(0);
697: }
699: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X,const PetscScalar **x)
700: {
704: PetscFree(*x);
705: return(0);
706: }
708: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
709: {
711: if (nx > 0) SETERRQ(PetscObjectComm((PetscObject)(*X)), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
712: return(0);
713: }
715: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
716: {
718: ops->duplicate = VecDuplicate_Nest;
719: ops->duplicatevecs = VecDuplicateVecs_Default;
720: ops->destroyvecs = VecDestroyVecs_Default;
721: ops->dot = VecDot_Nest;
722: ops->mdot = VecMDot_Nest;
723: ops->norm = VecNorm_Nest;
724: ops->tdot = VecTDot_Nest;
725: ops->mtdot = VecMTDot_Nest;
726: ops->scale = VecScale_Nest;
727: ops->copy = VecCopy_Nest;
728: ops->set = VecSet_Nest;
729: ops->swap = VecSwap_Nest;
730: ops->axpy = VecAXPY_Nest;
731: ops->axpby = VecAXPBY_Nest;
732: ops->maxpy = VecMAXPY_Nest;
733: ops->aypx = VecAYPX_Nest;
734: ops->waxpy = VecWAXPY_Nest;
735: ops->axpbypcz = NULL;
736: ops->pointwisemult = VecPointwiseMult_Nest;
737: ops->pointwisedivide = VecPointwiseDivide_Nest;
738: ops->setvalues = NULL;
739: ops->assemblybegin = VecAssemblyBegin_Nest;
740: ops->assemblyend = VecAssemblyEnd_Nest;
741: ops->getarray = VecGetArray_Nest;
742: ops->getsize = VecGetSize_Nest;
743: ops->getlocalsize = VecGetLocalSize_Nest;
744: ops->restorearray = VecRestoreArray_Nest;
745: ops->restorearrayread = VecRestoreArrayRead_Nest;
746: ops->max = VecMax_Nest;
747: ops->min = VecMin_Nest;
748: ops->setrandom = NULL;
749: ops->setoption = NULL;
750: ops->setvaluesblocked = NULL;
751: ops->destroy = VecDestroy_Nest;
752: ops->view = VecView_Nest;
753: ops->placearray = NULL;
754: ops->replacearray = NULL;
755: ops->dot_local = VecDot_Nest;
756: ops->tdot_local = VecTDot_Nest;
757: ops->norm_local = VecNorm_Nest;
758: ops->mdot_local = VecMDot_Nest;
759: ops->mtdot_local = VecMTDot_Nest;
760: ops->load = NULL;
761: ops->reciprocal = VecReciprocal_Nest;
762: ops->conjugate = VecConjugate_Nest;
763: ops->setlocaltoglobalmapping = NULL;
764: ops->setvalueslocal = NULL;
765: ops->resetarray = NULL;
766: ops->setfromoptions = NULL;
767: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
768: ops->load = NULL;
769: ops->pointwisemax = NULL;
770: ops->pointwisemaxabs = NULL;
771: ops->pointwisemin = NULL;
772: ops->getvalues = NULL;
773: ops->sqrt = NULL;
774: ops->abs = NULL;
775: ops->exp = NULL;
776: ops->shift = NULL;
777: ops->create = NULL;
778: ops->stridegather = NULL;
779: ops->stridescatter = NULL;
780: ops->dotnorm2 = VecDotNorm2_Nest;
781: ops->getsubvector = VecGetSubVector_Nest;
782: ops->restoresubvector = VecRestoreSubVector_Nest;
783: ops->axpbypcz = VecAXPBYPCZ_Nest;
784: ops->concatenate = VecConcatenate_Nest;
785: return(0);
786: }
788: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
789: {
790: Vec_Nest *b = (Vec_Nest*)x->data;
791: PetscInt i;
792: PetscInt row;
795: if (!m) return(0);
796: for (i=0; i<m; i++) {
797: row = idxm[i];
798: if (row >= b->nb) SETERRQ2(PetscObjectComm((PetscObject)x),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,b->nb-1);
799: vec[i] = b->v[row];
800: }
801: return(0);
802: }
804: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
805: {
809: VecNestGetSubVecs_Private(X,1,&idxm,sx);
810: return(0);
811: }
813: /*@
814: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
816: Not collective
818: Input Parameters:
819: + X - nest vector
820: - idxm - index of the vector within the nest
822: Output Parameter:
823: . sx - vector at index idxm within the nest
825: Notes:
827: Level: developer
829: .seealso: VecNestGetSize(), VecNestGetSubVecs()
830: @*/
831: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
832: {
836: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
837: return(0);
838: }
840: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
841: {
842: Vec_Nest *b = (Vec_Nest*)X->data;
845: if (N) *N = b->nb;
846: if (sx) *sx = b->v;
847: return(0);
848: }
850: /*@C
851: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
853: Not collective
855: Input Parameter:
856: . X - nest vector
858: Output Parameters:
859: + N - number of nested vecs
860: - sx - array of vectors
862: Notes:
863: The user should not free the array sx.
865: Fortran Notes:
866: The caller must allocate the array to hold the subvectors.
868: Level: developer
870: .seealso: VecNestGetSize(), VecNestGetSubVec()
871: @*/
872: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
873: {
877: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
878: return(0);
879: }
881: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
882: {
883: Vec_Nest *bx = (Vec_Nest*)X->data;
884: PetscInt i,offset=0,n=0,bs;
885: IS is;
887: PetscBool issame = PETSC_FALSE;
888: PetscInt N=0;
890: /* check if idxm < bx->nb */
891: if (idxm >= bx->nb) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",idxm,bx->nb);
894: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
895: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
896: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
898: /* check if we need to update the IS for the block */
899: offset = X->map->rstart;
900: for (i=0; i<idxm; i++) {
901: n=0;
902: VecGetLocalSize(bx->v[i],&n);
903: offset += n;
904: }
906: /* get the local size and block size */
907: VecGetLocalSize(x,&n);
908: VecGetBlockSize(x,&bs);
910: /* create the new IS */
911: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
912: ISSetBlockSize(is,bs);
914: /* check if they are equal */
915: ISEqual(is,bx->is[idxm],&issame);
917: if (!issame) {
918: /* The IS of given vector has a different layout compared to the existing block vector.
919: Destroy the existing reference and update the IS. */
920: ISDestroy(&bx->is[idxm]);
921: ISDuplicate(is,&bx->is[idxm]);
922: ISCopy(is,bx->is[idxm]);
924: offset += n;
925: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
926: for (i=idxm+1; i<bx->nb; i++) {
927: /* get the local size and block size */
928: VecGetLocalSize(bx->v[i],&n);
929: VecGetBlockSize(bx->v[i],&bs);
931: /* destroy the old and create the new IS */
932: ISDestroy(&bx->is[i]);
933: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
934: ISSetBlockSize(bx->is[i],bs);
936: offset += n;
937: }
939: n=0;
940: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
941: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
942: PetscLayoutSetSize(X->map,N);
943: PetscLayoutSetLocalSize(X->map,n);
944: }
946: ISDestroy(&is);
947: return(0);
948: }
950: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
951: {
955: VecNestSetSubVec_Private(X,idxm,sx);
956: return(0);
957: }
959: /*@
960: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
962: Not collective
964: Input Parameters:
965: + X - nest vector
966: . idxm - index of the vector within the nest vector
967: - sx - vector at index idxm within the nest vector
969: Notes:
970: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
972: Level: developer
974: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
975: @*/
976: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
977: {
981: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
982: return(0);
983: }
985: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
986: {
987: PetscInt i;
991: for (i=0; i<N; i++) {
992: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
993: }
994: return(0);
995: }
997: /*@C
998: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1000: Not collective
1002: Input Parameters:
1003: + X - nest vector
1004: . N - number of component vecs in sx
1005: . idxm - indices of component vecs that are to be replaced
1006: - sx - array of vectors
1008: Notes:
1009: The components in the vector array sx do not have to be of the same size as corresponding
1010: components in X. The user can also free the array "sx" after the call.
1012: Level: developer
1014: .seealso: VecNestGetSize(), VecNestGetSubVec()
1015: @*/
1016: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
1017: {
1021: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
1022: return(0);
1023: }
1025: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
1026: {
1027: Vec_Nest *b = (Vec_Nest*)X->data;
1030: *N = b->nb;
1031: return(0);
1032: }
1034: /*@
1035: VecNestGetSize - Returns the size of the nest vector.
1037: Not collective
1039: Input Parameter:
1040: . X - nest vector
1042: Output Parameter:
1043: . N - number of nested vecs
1045: Notes:
1047: Level: developer
1049: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
1050: @*/
1051: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
1052: {
1058: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
1059: return(0);
1060: }
1062: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
1063: {
1064: Vec_Nest *ctx = (Vec_Nest*)V->data;
1065: PetscInt i;
1069: if (ctx->setup_called) return(0);
1071: ctx->nb = nb;
1072: if (ctx->nb < 0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONG,"Cannot create VECNEST with < 0 blocks.");
1074: /* Create space */
1075: PetscMalloc1(ctx->nb,&ctx->v);
1076: for (i=0; i<ctx->nb; i++) {
1077: ctx->v[i] = x[i];
1078: PetscObjectReference((PetscObject)x[i]);
1079: /* Do not allocate memory for internal sub blocks */
1080: }
1082: PetscMalloc1(ctx->nb,&ctx->is);
1084: ctx->setup_called = PETSC_TRUE;
1085: return(0);
1086: }
1088: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
1089: {
1090: Vec_Nest *ctx = (Vec_Nest*)V->data;
1091: PetscInt i,offset,m,n,M,N;
1095: if (is) { /* Do some consistency checks and reference the is */
1096: offset = V->map->rstart;
1097: for (i=0; i<ctx->nb; i++) {
1098: ISGetSize(is[i],&M);
1099: VecGetSize(ctx->v[i],&N);
1100: 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);
1101: ISGetLocalSize(is[i],&m);
1102: VecGetLocalSize(ctx->v[i],&n);
1103: 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);
1104: if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
1105: PetscInt start;
1106: PetscBool contiguous;
1107: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1108: if (!contiguous) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D is not contiguous with layout of matching vector",i);
1109: if (start != 0) SETERRQ1(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Index set %D introduces overlap or a hole",i);
1110: }
1111: PetscObjectReference((PetscObject)is[i]);
1112: ctx->is[i] = is[i];
1113: offset += n;
1114: }
1115: } else { /* Create a contiguous ISStride for each entry */
1116: offset = V->map->rstart;
1117: for (i=0; i<ctx->nb; i++) {
1118: PetscInt bs;
1119: VecGetLocalSize(ctx->v[i],&n);
1120: VecGetBlockSize(ctx->v[i],&bs);
1121: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1122: ISSetBlockSize(ctx->is[i],bs);
1123: offset += n;
1124: }
1125: }
1126: return(0);
1127: }
1129: /*@C
1130: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1132: Collective on Vec
1134: Input Parameters:
1135: + comm - Communicator for the new Vec
1136: . nb - number of nested blocks
1137: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1138: - x - array of nb sub-vectors
1140: Output Parameter:
1141: . Y - new vector
1143: Level: advanced
1145: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1146: @*/
1147: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1148: {
1149: Vec V;
1150: Vec_Nest *s;
1151: PetscInt n,N;
1155: VecCreate(comm,&V);
1156: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1158: /* allocate and set pointer for implememtation data */
1159: PetscNew(&s);
1160: V->data = (void*)s;
1161: s->setup_called = PETSC_FALSE;
1162: s->nb = -1;
1163: s->v = NULL;
1165: VecSetUp_Nest_Private(V,nb,x);
1167: n = N = 0;
1168: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1169: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1170: PetscLayoutSetSize(V->map,N);
1171: PetscLayoutSetLocalSize(V->map,n);
1172: PetscLayoutSetBlockSize(V->map,1);
1173: PetscLayoutSetUp(V->map);
1175: VecSetUp_NestIS_Private(V,nb,is);
1177: VecNestSetOps_Private(V->ops);
1178: 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*/