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;
10: for (i=0; i<vs->nb; i++) {
12: VecAssemblyBegin(vs->v[i]);
13: }
14: return 0;
15: }
17: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
18: {
19: Vec_Nest *vs = (Vec_Nest*)v->data;
20: PetscInt i;
22: for (i=0; i<vs->nb; i++) {
23: VecAssemblyEnd(vs->v[i]);
24: }
25: return 0;
26: }
28: static PetscErrorCode VecDestroy_Nest(Vec v)
29: {
30: Vec_Nest *vs = (Vec_Nest*)v->data;
31: PetscInt i;
33: if (vs->v) {
34: for (i=0; i<vs->nb; i++) {
35: VecDestroy(&vs->v[i]);
36: }
37: PetscFree(vs->v);
38: }
39: for (i=0; i<vs->nb; i++) {
40: ISDestroy(&vs->is[i]);
41: }
42: PetscFree(vs->is);
43: PetscObjectComposeFunction((PetscObject)v,"",NULL);
44: PetscObjectComposeFunction((PetscObject)v,"",NULL);
45: PetscObjectComposeFunction((PetscObject)v,"",NULL);
46: PetscObjectComposeFunction((PetscObject)v,"",NULL);
47: PetscObjectComposeFunction((PetscObject)v,"",NULL);
49: PetscFree(vs);
50: return 0;
51: }
53: /* supports nested blocks */
54: static PetscErrorCode VecCopy_Nest(Vec x,Vec y)
55: {
56: Vec_Nest *bx = (Vec_Nest*)x->data;
57: Vec_Nest *by = (Vec_Nest*)y->data;
58: PetscInt i;
61: VecNestCheckCompatible2(x,1,y,2);
62: for (i=0; i<bx->nb; i++) {
63: VecCopy(bx->v[i],by->v[i]);
64: }
65: return 0;
66: }
68: /* supports nested blocks */
69: static PetscErrorCode VecDuplicate_Nest(Vec x,Vec *y)
70: {
71: Vec_Nest *bx = (Vec_Nest*)x->data;
72: Vec Y;
73: Vec *sub;
74: PetscInt i;
76: PetscMalloc1(bx->nb,&sub);
77: for (i=0; i<bx->nb; i++) {
78: VecDuplicate(bx->v[i],&sub[i]);
79: }
80: VecCreateNest(PetscObjectComm((PetscObject)x),bx->nb,bx->is,sub,&Y);
81: for (i=0; i<bx->nb; i++) {
82: VecDestroy(&sub[i]);
83: }
84: PetscFree(sub);
85: *y = Y;
86: return 0;
87: }
89: /* supports nested blocks */
90: static PetscErrorCode VecDot_Nest(Vec x,Vec y,PetscScalar *val)
91: {
92: Vec_Nest *bx = (Vec_Nest*)x->data;
93: Vec_Nest *by = (Vec_Nest*)y->data;
94: PetscInt i,nr;
95: PetscScalar x_dot_y,_val;
97: nr = bx->nb;
98: _val = 0.0;
99: for (i=0; i<nr; i++) {
100: VecDot(bx->v[i],by->v[i],&x_dot_y);
101: _val = _val + x_dot_y;
102: }
103: *val = _val;
104: return 0;
105: }
107: /* supports nested blocks */
108: static PetscErrorCode VecTDot_Nest(Vec x,Vec y,PetscScalar *val)
109: {
110: Vec_Nest *bx = (Vec_Nest*)x->data;
111: Vec_Nest *by = (Vec_Nest*)y->data;
112: PetscInt i,nr;
113: PetscScalar x_dot_y,_val;
115: nr = bx->nb;
116: _val = 0.0;
117: for (i=0; i<nr; i++) {
118: VecTDot(bx->v[i],by->v[i],&x_dot_y);
119: _val = _val + x_dot_y;
120: }
121: *val = _val;
122: return 0;
123: }
125: static PetscErrorCode VecDotNorm2_Nest(Vec x,Vec y,PetscScalar *dp, PetscScalar *nm)
126: {
127: Vec_Nest *bx = (Vec_Nest*)x->data;
128: Vec_Nest *by = (Vec_Nest*)y->data;
129: PetscInt i,nr;
130: PetscScalar x_dot_y,_dp,_nm;
131: PetscReal norm2_y;
133: nr = bx->nb;
134: _dp = 0.0;
135: _nm = 0.0;
136: for (i=0; i<nr; i++) {
137: VecDotNorm2(bx->v[i],by->v[i],&x_dot_y,&norm2_y);
138: _dp += x_dot_y;
139: _nm += norm2_y;
140: }
141: *dp = _dp;
142: *nm = _nm;
143: return 0;
144: }
146: static PetscErrorCode VecAXPY_Nest(Vec y,PetscScalar alpha,Vec x)
147: {
148: Vec_Nest *bx = (Vec_Nest*)x->data;
149: Vec_Nest *by = (Vec_Nest*)y->data;
150: PetscInt i,nr;
152: nr = bx->nb;
153: for (i=0; i<nr; i++) {
154: VecAXPY(by->v[i],alpha,bx->v[i]);
155: }
156: return 0;
157: }
159: static PetscErrorCode VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)
160: {
161: Vec_Nest *bx = (Vec_Nest*)x->data;
162: Vec_Nest *by = (Vec_Nest*)y->data;
163: PetscInt i,nr;
165: nr = bx->nb;
166: for (i=0; i<nr; i++) {
167: VecAYPX(by->v[i],alpha,bx->v[i]);
168: }
169: return 0;
170: }
172: static PetscErrorCode VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
173: {
174: Vec_Nest *bx = (Vec_Nest*)x->data;
175: Vec_Nest *by = (Vec_Nest*)y->data;
176: PetscInt i,nr;
178: nr = bx->nb;
179: for (i=0; i<nr; i++) {
180: VecAXPBY(by->v[i],alpha,beta,bx->v[i]);
181: }
182: return 0;
183: }
185: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)
186: {
187: Vec_Nest *bx = (Vec_Nest*)x->data;
188: Vec_Nest *by = (Vec_Nest*)y->data;
189: Vec_Nest *bz = (Vec_Nest*)z->data;
190: PetscInt i,nr;
192: nr = bx->nb;
193: for (i=0; i<nr; i++) {
194: VecAXPBYPCZ(bz->v[i],alpha,beta,gamma,bx->v[i],by->v[i]);
195: }
196: return 0;
197: }
199: static PetscErrorCode VecScale_Nest(Vec x,PetscScalar alpha)
200: {
201: Vec_Nest *bx = (Vec_Nest*)x->data;
202: PetscInt i,nr;
204: nr = bx->nb;
205: for (i=0; i<nr; i++) {
206: VecScale(bx->v[i],alpha);
207: }
208: return 0;
209: }
211: static PetscErrorCode VecPointwiseMult_Nest(Vec w,Vec x,Vec y)
212: {
213: Vec_Nest *bx = (Vec_Nest*)x->data;
214: Vec_Nest *by = (Vec_Nest*)y->data;
215: Vec_Nest *bw = (Vec_Nest*)w->data;
216: PetscInt i,nr;
218: VecNestCheckCompatible3(w,1,x,2,y,3);
219: nr = bx->nb;
220: for (i=0; i<nr; i++) {
221: VecPointwiseMult(bw->v[i],bx->v[i],by->v[i]);
222: }
223: return 0;
224: }
226: static PetscErrorCode VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)
227: {
228: Vec_Nest *bx = (Vec_Nest*)x->data;
229: Vec_Nest *by = (Vec_Nest*)y->data;
230: Vec_Nest *bw = (Vec_Nest*)w->data;
231: PetscInt i,nr;
233: VecNestCheckCompatible3(w,1,x,2,y,3);
235: nr = bx->nb;
236: for (i=0; i<nr; i++) {
237: VecPointwiseDivide(bw->v[i],bx->v[i],by->v[i]);
238: }
239: return 0;
240: }
242: static PetscErrorCode VecReciprocal_Nest(Vec x)
243: {
244: Vec_Nest *bx = (Vec_Nest*)x->data;
245: PetscInt i,nr;
247: nr = bx->nb;
248: for (i=0; i<nr; i++) {
249: VecReciprocal(bx->v[i]);
250: }
251: return 0;
252: }
254: static PetscErrorCode VecNorm_Nest(Vec xin,NormType type,PetscReal *z)
255: {
256: Vec_Nest *bx = (Vec_Nest*)xin->data;
257: PetscInt i,nr;
258: PetscReal z_i;
259: PetscReal _z;
261: nr = bx->nb;
262: _z = 0.0;
264: if (type == NORM_2) {
265: PetscScalar dot;
266: VecDot(xin,xin,&dot);
267: _z = PetscAbsScalar(PetscSqrtScalar(dot));
268: } else if (type == NORM_1) {
269: for (i=0; i<nr; i++) {
270: VecNorm(bx->v[i],type,&z_i);
271: _z = _z + z_i;
272: }
273: } else if (type == NORM_INFINITY) {
274: for (i=0; i<nr; i++) {
275: VecNorm(bx->v[i],type,&z_i);
276: if (z_i > _z) _z = z_i;
277: }
278: }
280: *z = _z;
281: return 0;
282: }
284: static PetscErrorCode VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
285: {
286: PetscInt v;
288: for (v=0; v<nv; v++) {
289: /* Do axpy on each vector,v */
290: VecAXPY(y,alpha[v],x[v]);
291: }
292: return 0;
293: }
295: static PetscErrorCode VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
296: {
297: PetscInt j;
299: for (j=0; j<nv; j++) {
300: VecDot(x,y[j],&val[j]);
301: }
302: return 0;
303: }
305: static PetscErrorCode VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
306: {
307: PetscInt j;
309: for (j=0; j<nv; j++) {
310: VecTDot(x,y[j],&val[j]);
311: }
312: return 0;
313: }
315: static PetscErrorCode VecSet_Nest(Vec x,PetscScalar alpha)
316: {
317: Vec_Nest *bx = (Vec_Nest*)x->data;
318: PetscInt i,nr;
320: nr = bx->nb;
321: for (i=0; i<nr; i++) {
322: VecSet(bx->v[i],alpha);
323: }
324: return 0;
325: }
327: static PetscErrorCode VecConjugate_Nest(Vec x)
328: {
329: Vec_Nest *bx = (Vec_Nest*)x->data;
330: PetscInt j,nr;
332: nr = bx->nb;
333: for (j=0; j<nr; j++) {
334: VecConjugate(bx->v[j]);
335: }
336: return 0;
337: }
339: static PetscErrorCode VecSwap_Nest(Vec x,Vec y)
340: {
341: Vec_Nest *bx = (Vec_Nest*)x->data;
342: Vec_Nest *by = (Vec_Nest*)y->data;
343: PetscInt i,nr;
345: VecNestCheckCompatible2(x,1,y,2);
346: nr = bx->nb;
347: for (i=0; i<nr; i++) {
348: VecSwap(bx->v[i],by->v[i]);
349: }
350: return 0;
351: }
353: static PetscErrorCode VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)
354: {
355: Vec_Nest *bx = (Vec_Nest*)x->data;
356: Vec_Nest *by = (Vec_Nest*)y->data;
357: Vec_Nest *bw = (Vec_Nest*)w->data;
358: PetscInt i,nr;
360: VecNestCheckCompatible3(w,1,x,3,y,4);
362: nr = bx->nb;
363: for (i=0; i<nr; i++) {
364: VecWAXPY(bw->v[i],alpha,bx->v[i],by->v[i]);
365: }
366: return 0;
367: }
369: static PetscErrorCode VecMax_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *max)
370: {
371: Vec_Nest *bx;
372: PetscInt i,nr;
373: PetscBool isnest;
374: PetscInt L;
375: PetscInt _entry_loc;
376: PetscReal _entry_val;
378: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
379: if (!isnest) {
380: /* Not nest */
381: VecMax(x,&_entry_loc,&_entry_val);
382: if (_entry_val > *max) {
383: *max = _entry_val;
384: if (p) *p = _entry_loc + *cnt;
385: }
386: VecGetSize(x,&L);
387: *cnt = *cnt + L;
388: return 0;
389: }
391: /* Otherwise we have a nest */
392: bx = (Vec_Nest*)x->data;
393: nr = bx->nb;
395: /* now descend recursively */
396: for (i=0; i<nr; i++) {
397: VecMax_Nest_Recursive(bx->v[i],cnt,p,max);
398: }
399: return 0;
400: }
402: /* supports nested blocks */
403: static PetscErrorCode VecMax_Nest(Vec x,PetscInt *p,PetscReal *max)
404: {
405: PetscInt cnt;
407: cnt = 0;
408: if (p) *p = 0;
409: *max = PETSC_MIN_REAL;
410: VecMax_Nest_Recursive(x,&cnt,p,max);
411: return 0;
412: }
414: static PetscErrorCode VecMin_Nest_Recursive(Vec x,PetscInt *cnt,PetscInt *p,PetscReal *min)
415: {
416: Vec_Nest *bx = (Vec_Nest*)x->data;
417: PetscInt i,nr,L,_entry_loc;
418: PetscBool isnest;
419: PetscReal _entry_val;
421: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
422: if (!isnest) {
423: /* Not nest */
424: VecMin(x,&_entry_loc,&_entry_val);
425: if (_entry_val < *min) {
426: *min = _entry_val;
427: if (p) *p = _entry_loc + *cnt;
428: }
429: VecGetSize(x,&L);
430: *cnt = *cnt + L;
431: return 0;
432: }
434: /* Otherwise we have a nest */
435: nr = bx->nb;
437: /* now descend recursively */
438: for (i=0; i<nr; i++) {
439: VecMin_Nest_Recursive(bx->v[i],cnt,p,min);
440: }
441: return 0;
442: }
444: static PetscErrorCode VecMin_Nest(Vec x,PetscInt *p,PetscReal *min)
445: {
446: PetscInt cnt;
448: cnt = 0;
449: if (p) *p = 0;
450: *min = PETSC_MAX_REAL;
451: VecMin_Nest_Recursive(x,&cnt,p,min);
452: return 0;
453: }
455: /* supports nested blocks */
456: static PetscErrorCode VecView_Nest(Vec x,PetscViewer viewer)
457: {
458: Vec_Nest *bx = (Vec_Nest*)x->data;
459: PetscBool isascii;
460: PetscInt i;
462: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
463: if (isascii) {
464: PetscViewerASCIIPushTab(viewer);
465: PetscViewerASCIIPrintf(viewer,"VecNest, rows=%" PetscInt_FMT ", structure: \n",bx->nb);
466: for (i=0; i<bx->nb; i++) {
467: VecType type;
468: char name[256] = "",prefix[256] = "";
469: PetscInt NR;
471: VecGetSize(bx->v[i],&NR);
472: VecGetType(bx->v[i],&type);
473: if (((PetscObject)bx->v[i])->name) PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bx->v[i])->name);
474: if (((PetscObject)bx->v[i])->prefix) PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bx->v[i])->prefix);
476: PetscViewerASCIIPrintf(viewer,"(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n",i,name,prefix,type,NR);
478: PetscViewerASCIIPushTab(viewer); /* push1 */
479: VecView(bx->v[i],viewer);
480: PetscViewerASCIIPopTab(viewer); /* pop1 */
481: }
482: PetscViewerASCIIPopTab(viewer); /* pop0 */
483: }
484: return 0;
485: }
487: static PetscErrorCode VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt *L)
488: {
489: Vec_Nest *bx;
490: PetscInt size,i,nr;
491: PetscBool isnest;
493: PetscObjectTypeCompare((PetscObject)x,VECNEST,&isnest);
494: if (!isnest) {
495: /* Not nest */
496: if (globalsize) VecGetSize(x,&size);
497: else VecGetLocalSize(x,&size);
498: *L = *L + size;
499: return 0;
500: }
502: /* Otherwise we have a nest */
503: bx = (Vec_Nest*)x->data;
504: nr = bx->nb;
506: /* now descend recursively */
507: for (i=0; i<nr; i++) {
508: VecSize_Nest_Recursive(bx->v[i],globalsize,L);
509: }
510: return 0;
511: }
513: /* Returns the sum of the global size of all the consituent vectors in the nest */
514: static PetscErrorCode VecGetSize_Nest(Vec x,PetscInt *N)
515: {
516: *N = x->map->N;
517: return 0;
518: }
520: /* Returns the sum of the local size of all the consituent vectors in the nest */
521: static PetscErrorCode VecGetLocalSize_Nest(Vec x,PetscInt *n)
522: {
523: *n = x->map->n;
524: return 0;
525: }
527: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal *max)
528: {
529: Vec_Nest *bx = (Vec_Nest*)x->data;
530: Vec_Nest *by = (Vec_Nest*)y->data;
531: PetscInt i,nr;
532: PetscReal local_max,m;
534: VecNestCheckCompatible2(x,1,y,2);
535: nr = bx->nb;
536: m = 0.0;
537: for (i=0; i<nr; i++) {
538: VecMaxPointwiseDivide(bx->v[i],by->v[i],&local_max);
539: if (local_max > m) m = local_max;
540: }
541: *max = m;
542: return 0;
543: }
545: static PetscErrorCode VecGetSubVector_Nest(Vec X,IS is,Vec *x)
546: {
547: Vec_Nest *bx = (Vec_Nest*)X->data;
548: PetscInt i;
550: *x = NULL;
551: for (i=0; i<bx->nb; i++) {
552: PetscBool issame = PETSC_FALSE;
553: ISEqual(is,bx->is[i],&issame);
554: if (issame) {
555: *x = bx->v[i];
556: PetscObjectReference((PetscObject)(*x));
557: break;
558: }
559: }
561: return 0;
562: }
564: static PetscErrorCode VecRestoreSubVector_Nest(Vec X,IS is,Vec *x)
565: {
566: VecDestroy(x);
567: return 0;
568: }
570: static PetscErrorCode VecGetArray_Nest(Vec X,PetscScalar **x)
571: {
572: Vec_Nest *bx = (Vec_Nest*)X->data;
573: PetscInt i,m,rstart,rend;
575: VecGetOwnershipRange(X,&rstart,&rend);
576: VecGetLocalSize(X,&m);
577: PetscMalloc1(m,x);
578: for (i=0; i<bx->nb; i++) {
579: Vec subvec = bx->v[i];
580: IS isy = bx->is[i];
581: PetscInt j,sm;
582: const PetscInt *ixy;
583: const PetscScalar *y;
584: VecGetLocalSize(subvec,&sm);
585: VecGetArrayRead(subvec,&y);
586: ISGetIndices(isy,&ixy);
587: for (j=0; j<sm; j++) {
588: PetscInt ix = ixy[j];
590: (*x)[ix-rstart] = y[j];
591: }
592: ISRestoreIndices(isy,&ixy);
593: VecRestoreArrayRead(subvec,&y);
594: }
595: return 0;
596: }
598: static PetscErrorCode VecRestoreArray_Nest(Vec X,PetscScalar **x)
599: {
600: Vec_Nest *bx = (Vec_Nest*)X->data;
601: PetscInt i,m,rstart,rend;
603: VecGetOwnershipRange(X,&rstart,&rend);
604: VecGetLocalSize(X,&m);
605: for (i=0; i<bx->nb; i++) {
606: Vec subvec = bx->v[i];
607: IS isy = bx->is[i];
608: PetscInt j,sm;
609: const PetscInt *ixy;
610: PetscScalar *y;
611: VecGetLocalSize(subvec,&sm);
612: VecGetArray(subvec,&y);
613: ISGetIndices(isy,&ixy);
614: for (j=0; j<sm; j++) {
615: PetscInt ix = ixy[j];
617: y[j] = (*x)[ix-rstart];
618: }
619: ISRestoreIndices(isy,&ixy);
620: VecRestoreArray(subvec,&y);
621: }
622: PetscFree(*x);
623: return 0;
624: }
626: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X,const PetscScalar **x)
627: {
628: PetscFree(*x);
629: return 0;
630: }
632: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
633: {
635: return 0;
636: }
638: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
639: {
640: ops->duplicate = VecDuplicate_Nest;
641: ops->duplicatevecs = VecDuplicateVecs_Default;
642: ops->destroyvecs = VecDestroyVecs_Default;
643: ops->dot = VecDot_Nest;
644: ops->mdot = VecMDot_Nest;
645: ops->norm = VecNorm_Nest;
646: ops->tdot = VecTDot_Nest;
647: ops->mtdot = VecMTDot_Nest;
648: ops->scale = VecScale_Nest;
649: ops->copy = VecCopy_Nest;
650: ops->set = VecSet_Nest;
651: ops->swap = VecSwap_Nest;
652: ops->axpy = VecAXPY_Nest;
653: ops->axpby = VecAXPBY_Nest;
654: ops->maxpy = VecMAXPY_Nest;
655: ops->aypx = VecAYPX_Nest;
656: ops->waxpy = VecWAXPY_Nest;
657: ops->axpbypcz = NULL;
658: ops->pointwisemult = VecPointwiseMult_Nest;
659: ops->pointwisedivide = VecPointwiseDivide_Nest;
660: ops->setvalues = NULL;
661: ops->assemblybegin = VecAssemblyBegin_Nest;
662: ops->assemblyend = VecAssemblyEnd_Nest;
663: ops->getarray = VecGetArray_Nest;
664: ops->getsize = VecGetSize_Nest;
665: ops->getlocalsize = VecGetLocalSize_Nest;
666: ops->restorearray = VecRestoreArray_Nest;
667: ops->restorearrayread = VecRestoreArrayRead_Nest;
668: ops->max = VecMax_Nest;
669: ops->min = VecMin_Nest;
670: ops->setrandom = NULL;
671: ops->setoption = NULL;
672: ops->setvaluesblocked = NULL;
673: ops->destroy = VecDestroy_Nest;
674: ops->view = VecView_Nest;
675: ops->placearray = NULL;
676: ops->replacearray = NULL;
677: ops->dot_local = VecDot_Nest;
678: ops->tdot_local = VecTDot_Nest;
679: ops->norm_local = VecNorm_Nest;
680: ops->mdot_local = VecMDot_Nest;
681: ops->mtdot_local = VecMTDot_Nest;
682: ops->load = NULL;
683: ops->reciprocal = VecReciprocal_Nest;
684: ops->conjugate = VecConjugate_Nest;
685: ops->setlocaltoglobalmapping = NULL;
686: ops->setvalueslocal = NULL;
687: ops->resetarray = NULL;
688: ops->setfromoptions = NULL;
689: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
690: ops->load = NULL;
691: ops->pointwisemax = NULL;
692: ops->pointwisemaxabs = NULL;
693: ops->pointwisemin = NULL;
694: ops->getvalues = NULL;
695: ops->sqrt = NULL;
696: ops->abs = NULL;
697: ops->exp = NULL;
698: ops->shift = NULL;
699: ops->create = NULL;
700: ops->stridegather = NULL;
701: ops->stridescatter = NULL;
702: ops->dotnorm2 = VecDotNorm2_Nest;
703: ops->getsubvector = VecGetSubVector_Nest;
704: ops->restoresubvector = VecRestoreSubVector_Nest;
705: ops->axpbypcz = VecAXPBYPCZ_Nest;
706: ops->concatenate = VecConcatenate_Nest;
707: return 0;
708: }
710: static PetscErrorCode VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])
711: {
712: Vec_Nest *b = (Vec_Nest*)x->data;
713: PetscInt i;
714: PetscInt row;
716: if (!m) return 0;
717: for (i=0; i<m; i++) {
718: row = idxm[i];
720: vec[i] = b->v[row];
721: }
722: return 0;
723: }
725: PetscErrorCode VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec *sx)
726: {
727: VecNestGetSubVecs_Private(X,1,&idxm,sx);
728: return 0;
729: }
731: /*@
732: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
734: Not collective
736: Input Parameters:
737: + X - nest vector
738: - idxm - index of the vector within the nest
740: Output Parameter:
741: . sx - vector at index idxm within the nest
743: Notes:
745: Level: developer
747: .seealso: VecNestGetSize(), VecNestGetSubVecs()
748: @*/
749: PetscErrorCode VecNestGetSubVec(Vec X,PetscInt idxm,Vec *sx)
750: {
751: PetscUseMethod(X,"VecNestGetSubVec_C",(Vec,PetscInt,Vec*),(X,idxm,sx));
752: return 0;
753: }
755: PetscErrorCode VecNestGetSubVecs_Nest(Vec X,PetscInt *N,Vec **sx)
756: {
757: Vec_Nest *b = (Vec_Nest*)X->data;
759: if (N) *N = b->nb;
760: if (sx) *sx = b->v;
761: return 0;
762: }
764: /*@C
765: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
767: Not collective
769: Input Parameter:
770: . X - nest vector
772: Output Parameters:
773: + N - number of nested vecs
774: - sx - array of vectors
776: Notes:
777: The user should not free the array sx.
779: Fortran Notes:
780: The caller must allocate the array to hold the subvectors.
782: Level: developer
784: .seealso: VecNestGetSize(), VecNestGetSubVec()
785: @*/
786: PetscErrorCode VecNestGetSubVecs(Vec X,PetscInt *N,Vec **sx)
787: {
788: PetscUseMethod(X,"VecNestGetSubVecs_C",(Vec,PetscInt*,Vec**),(X,N,sx));
789: return 0;
790: }
792: static PetscErrorCode VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)
793: {
794: Vec_Nest *bx = (Vec_Nest*)X->data;
795: PetscInt i,offset=0,n=0,bs;
796: IS is;
797: PetscBool issame = PETSC_FALSE;
798: PetscInt N=0;
800: /* check if idxm < bx->nb */
803: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
804: VecDuplicate(x,&bx->v[idxm]); /* duplicate the layout of given vector */
805: VecCopy(x,bx->v[idxm]); /* copy the contents of the given vector */
807: /* check if we need to update the IS for the block */
808: offset = X->map->rstart;
809: for (i=0; i<idxm; i++) {
810: n=0;
811: VecGetLocalSize(bx->v[i],&n);
812: offset += n;
813: }
815: /* get the local size and block size */
816: VecGetLocalSize(x,&n);
817: VecGetBlockSize(x,&bs);
819: /* create the new IS */
820: ISCreateStride(PetscObjectComm((PetscObject)x),n,offset,1,&is);
821: ISSetBlockSize(is,bs);
823: /* check if they are equal */
824: ISEqual(is,bx->is[idxm],&issame);
826: if (!issame) {
827: /* The IS of given vector has a different layout compared to the existing block vector.
828: Destroy the existing reference and update the IS. */
829: ISDestroy(&bx->is[idxm]);
830: ISDuplicate(is,&bx->is[idxm]);
831: ISCopy(is,bx->is[idxm]);
833: offset += n;
834: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
835: for (i=idxm+1; i<bx->nb; i++) {
836: /* get the local size and block size */
837: VecGetLocalSize(bx->v[i],&n);
838: VecGetBlockSize(bx->v[i],&bs);
840: /* destroy the old and create the new IS */
841: ISDestroy(&bx->is[i]);
842: ISCreateStride(((PetscObject)bx->v[i])->comm,n,offset,1,&bx->is[i]);
843: ISSetBlockSize(bx->is[i],bs);
845: offset += n;
846: }
848: n=0;
849: VecSize_Nest_Recursive(X,PETSC_TRUE,&N);
850: VecSize_Nest_Recursive(X,PETSC_FALSE,&n);
851: PetscLayoutSetSize(X->map,N);
852: PetscLayoutSetLocalSize(X->map,n);
853: }
855: ISDestroy(&is);
856: return 0;
857: }
859: PetscErrorCode VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)
860: {
861: VecNestSetSubVec_Private(X,idxm,sx);
862: return 0;
863: }
865: /*@
866: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
868: Not collective
870: Input Parameters:
871: + X - nest vector
872: . idxm - index of the vector within the nest vector
873: - sx - vector at index idxm within the nest vector
875: Notes:
876: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
878: Level: developer
880: .seealso: VecNestSetSubVecs(), VecNestGetSubVec()
881: @*/
882: PetscErrorCode VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)
883: {
884: PetscUseMethod(X,"VecNestSetSubVec_C",(Vec,PetscInt,Vec),(X,idxm,sx));
885: return 0;
886: }
888: PetscErrorCode VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
889: {
890: PetscInt i;
892: for (i=0; i<N; i++) {
893: VecNestSetSubVec_Private(X,idxm[i],sx[i]);
894: }
895: return 0;
896: }
898: /*@C
899: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
901: Not collective
903: Input Parameters:
904: + X - nest vector
905: . N - number of component vecs in sx
906: . idxm - indices of component vecs that are to be replaced
907: - sx - array of vectors
909: Notes:
910: The components in the vector array sx do not have to be of the same size as corresponding
911: components in X. The user can also free the array "sx" after the call.
913: Level: developer
915: .seealso: VecNestGetSize(), VecNestGetSubVec()
916: @*/
917: PetscErrorCode VecNestSetSubVecs(Vec X,PetscInt N,PetscInt *idxm,Vec *sx)
918: {
919: PetscUseMethod(X,"VecNestSetSubVecs_C",(Vec,PetscInt,PetscInt*,Vec*),(X,N,idxm,sx));
920: return 0;
921: }
923: PetscErrorCode VecNestGetSize_Nest(Vec X,PetscInt *N)
924: {
925: Vec_Nest *b = (Vec_Nest*)X->data;
927: *N = b->nb;
928: return 0;
929: }
931: /*@
932: VecNestGetSize - Returns the size of the nest vector.
934: Not collective
936: Input Parameter:
937: . X - nest vector
939: Output Parameter:
940: . N - number of nested vecs
942: Notes:
944: Level: developer
946: .seealso: VecNestGetSubVec(), VecNestGetSubVecs()
947: @*/
948: PetscErrorCode VecNestGetSize(Vec X,PetscInt *N)
949: {
952: PetscUseMethod(X,"VecNestGetSize_C",(Vec,PetscInt*),(X,N));
953: return 0;
954: }
956: static PetscErrorCode VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])
957: {
958: Vec_Nest *ctx = (Vec_Nest*)V->data;
959: PetscInt i;
961: if (ctx->setup_called) return 0;
963: ctx->nb = nb;
966: /* Create space */
967: PetscMalloc1(ctx->nb,&ctx->v);
968: for (i=0; i<ctx->nb; i++) {
969: ctx->v[i] = x[i];
970: PetscObjectReference((PetscObject)x[i]);
971: /* Do not allocate memory for internal sub blocks */
972: }
974: PetscMalloc1(ctx->nb,&ctx->is);
976: ctx->setup_called = PETSC_TRUE;
977: return 0;
978: }
980: static PetscErrorCode VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])
981: {
982: Vec_Nest *ctx = (Vec_Nest*)V->data;
983: PetscInt i,offset,m,n,M,N;
985: if (is) { /* Do some consistency checks and reference the is */
986: offset = V->map->rstart;
987: for (i=0; i<ctx->nb; i++) {
988: ISGetSize(is[i],&M);
989: VecGetSize(ctx->v[i],&N);
991: ISGetLocalSize(is[i],&m);
992: VecGetLocalSize(ctx->v[i],&n);
994: if (PetscDefined(USE_DEBUG)) { /* This test can be expensive */
995: PetscInt start;
996: PetscBool contiguous;
997: ISContiguousLocal(is[i],offset,offset+n,&start,&contiguous);
1000: }
1001: PetscObjectReference((PetscObject)is[i]);
1002: ctx->is[i] = is[i];
1003: offset += n;
1004: }
1005: } else { /* Create a contiguous ISStride for each entry */
1006: offset = V->map->rstart;
1007: for (i=0; i<ctx->nb; i++) {
1008: PetscInt bs;
1009: VecGetLocalSize(ctx->v[i],&n);
1010: VecGetBlockSize(ctx->v[i],&bs);
1011: ISCreateStride(((PetscObject)ctx->v[i])->comm,n,offset,1,&ctx->is[i]);
1012: ISSetBlockSize(ctx->is[i],bs);
1013: offset += n;
1014: }
1015: }
1016: return 0;
1017: }
1019: /*@C
1020: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1022: Collective on Vec
1024: Input Parameters:
1025: + comm - Communicator for the new Vec
1026: . nb - number of nested blocks
1027: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1028: - x - array of nb sub-vectors
1030: Output Parameter:
1031: . Y - new vector
1033: Level: advanced
1035: .seealso: VecCreate(), MatCreateNest(), DMSetVecType(), VECNEST
1036: @*/
1037: PetscErrorCode VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec *Y)
1038: {
1039: Vec V;
1040: Vec_Nest *s;
1041: PetscInt n,N;
1043: VecCreate(comm,&V);
1044: PetscObjectChangeTypeName((PetscObject)V,VECNEST);
1046: /* allocate and set pointer for implememtation data */
1047: PetscNew(&s);
1048: V->data = (void*)s;
1049: s->setup_called = PETSC_FALSE;
1050: s->nb = -1;
1051: s->v = NULL;
1053: VecSetUp_Nest_Private(V,nb,x);
1055: n = N = 0;
1056: VecSize_Nest_Recursive(V,PETSC_TRUE,&N);
1057: VecSize_Nest_Recursive(V,PETSC_FALSE,&n);
1058: PetscLayoutSetSize(V->map,N);
1059: PetscLayoutSetLocalSize(V->map,n);
1060: PetscLayoutSetBlockSize(V->map,1);
1061: PetscLayoutSetUp(V->map);
1063: VecSetUp_NestIS_Private(V,nb,is);
1065: VecNestSetOps_Private(V->ops);
1066: V->petscnative = PETSC_FALSE;
1068: /* expose block api's */
1069: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVec_C",VecNestGetSubVec_Nest);
1070: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSubVecs_C",VecNestGetSubVecs_Nest);
1071: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVec_C",VecNestSetSubVec_Nest);
1072: PetscObjectComposeFunction((PetscObject)V,"VecNestSetSubVecs_C",VecNestSetSubVecs_Nest);
1073: PetscObjectComposeFunction((PetscObject)V,"VecNestGetSize_C",VecNestGetSize_Nest);
1075: *Y = V;
1076: return 0;
1077: }
1079: /*MC
1080: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1082: Level: intermediate
1084: Notes:
1085: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1086: It is usually used with MATNEST and DMComposite via DMSetVecType().
1088: .seealso: VecCreate(), VecType, VecCreateNest(), MatCreateNest()
1089: M*/