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++) VecAssemblyEnd(vs->v[i]);
23: return 0;
24: }
26: static PetscErrorCode VecDestroy_Nest(Vec v)
27: {
28: Vec_Nest *vs = (Vec_Nest *)v->data;
29: PetscInt i;
31: if (vs->v) {
32: for (i = 0; i < vs->nb; i++) VecDestroy(&vs->v[i]);
33: PetscFree(vs->v);
34: }
35: for (i = 0; i < vs->nb; i++) ISDestroy(&vs->is[i]);
36: PetscFree(vs->is);
37: PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVec_C", NULL);
38: PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecs_C", NULL);
39: PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVec_C", NULL);
40: PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL);
41: PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL);
43: PetscFree(vs);
44: return 0;
45: }
47: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
48: {
49: Vec_Nest *bx = (Vec_Nest *)x->data;
50: Vec_Nest *by = (Vec_Nest *)y->data;
51: PetscInt i;
54: VecNestCheckCompatible2(x, 1, y, 2);
55: for (i = 0; i < bx->nb; i++) VecCopy(bx->v[i], by->v[i]);
56: return 0;
57: }
59: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
60: {
61: Vec_Nest *bx = (Vec_Nest *)x->data;
62: Vec Y;
63: Vec *sub;
64: PetscInt i;
66: PetscMalloc1(bx->nb, &sub);
67: for (i = 0; i < bx->nb; i++) VecDuplicate(bx->v[i], &sub[i]);
68: VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y);
69: for (i = 0; i < bx->nb; i++) VecDestroy(&sub[i]);
70: PetscFree(sub);
71: *y = Y;
72: return 0;
73: }
75: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
76: {
77: Vec_Nest *bx = (Vec_Nest *)x->data;
78: Vec_Nest *by = (Vec_Nest *)y->data;
79: PetscInt i, nr;
80: PetscScalar x_dot_y, _val;
82: nr = bx->nb;
83: _val = 0.0;
84: for (i = 0; i < nr; i++) {
85: VecDot(bx->v[i], by->v[i], &x_dot_y);
86: _val = _val + x_dot_y;
87: }
88: *val = _val;
89: return 0;
90: }
92: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
93: {
94: Vec_Nest *bx = (Vec_Nest *)x->data;
95: Vec_Nest *by = (Vec_Nest *)y->data;
96: PetscInt i, nr;
97: PetscScalar x_dot_y, _val;
99: nr = bx->nb;
100: _val = 0.0;
101: for (i = 0; i < nr; i++) {
102: VecTDot(bx->v[i], by->v[i], &x_dot_y);
103: _val = _val + x_dot_y;
104: }
105: *val = _val;
106: return 0;
107: }
109: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
110: {
111: Vec_Nest *bx = (Vec_Nest *)x->data;
112: Vec_Nest *by = (Vec_Nest *)y->data;
113: PetscInt i, nr;
114: PetscScalar x_dot_y, _dp, _nm;
115: PetscReal norm2_y;
117: nr = bx->nb;
118: _dp = 0.0;
119: _nm = 0.0;
120: for (i = 0; i < nr; i++) {
121: VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y);
122: _dp += x_dot_y;
123: _nm += norm2_y;
124: }
125: *dp = _dp;
126: *nm = _nm;
127: return 0;
128: }
130: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
131: {
132: Vec_Nest *bx = (Vec_Nest *)x->data;
133: Vec_Nest *by = (Vec_Nest *)y->data;
134: PetscInt i, nr;
136: nr = bx->nb;
137: for (i = 0; i < nr; i++) VecAXPY(by->v[i], alpha, bx->v[i]);
138: return 0;
139: }
141: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
142: {
143: Vec_Nest *bx = (Vec_Nest *)x->data;
144: Vec_Nest *by = (Vec_Nest *)y->data;
145: PetscInt i, nr;
147: nr = bx->nb;
148: for (i = 0; i < nr; i++) VecAYPX(by->v[i], alpha, bx->v[i]);
149: return 0;
150: }
152: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
153: {
154: Vec_Nest *bx = (Vec_Nest *)x->data;
155: Vec_Nest *by = (Vec_Nest *)y->data;
156: PetscInt i, nr;
158: nr = bx->nb;
159: for (i = 0; i < nr; i++) VecAXPBY(by->v[i], alpha, beta, bx->v[i]);
160: return 0;
161: }
163: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
164: {
165: Vec_Nest *bx = (Vec_Nest *)x->data;
166: Vec_Nest *by = (Vec_Nest *)y->data;
167: Vec_Nest *bz = (Vec_Nest *)z->data;
168: PetscInt i, nr;
170: nr = bx->nb;
171: for (i = 0; i < nr; i++) VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]);
172: return 0;
173: }
175: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
176: {
177: Vec_Nest *bx = (Vec_Nest *)x->data;
178: PetscInt i, nr;
180: nr = bx->nb;
181: for (i = 0; i < nr; i++) VecScale(bx->v[i], alpha);
182: return 0;
183: }
185: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
186: {
187: Vec_Nest *bx = (Vec_Nest *)x->data;
188: Vec_Nest *by = (Vec_Nest *)y->data;
189: Vec_Nest *bw = (Vec_Nest *)w->data;
190: PetscInt i, nr;
192: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
193: nr = bx->nb;
194: for (i = 0; i < nr; i++) VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]);
195: return 0;
196: }
198: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
199: {
200: Vec_Nest *bx = (Vec_Nest *)x->data;
201: Vec_Nest *by = (Vec_Nest *)y->data;
202: Vec_Nest *bw = (Vec_Nest *)w->data;
203: PetscInt i, nr;
205: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
207: nr = bx->nb;
208: for (i = 0; i < nr; i++) VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]);
209: return 0;
210: }
212: static PetscErrorCode VecReciprocal_Nest(Vec x)
213: {
214: Vec_Nest *bx = (Vec_Nest *)x->data;
215: PetscInt i, nr;
217: nr = bx->nb;
218: for (i = 0; i < nr; i++) VecReciprocal(bx->v[i]);
219: return 0;
220: }
222: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
223: {
224: Vec_Nest *bx = (Vec_Nest *)xin->data;
225: PetscInt i, nr;
226: PetscReal z_i;
227: PetscReal _z;
229: nr = bx->nb;
230: _z = 0.0;
232: if (type == NORM_2) {
233: PetscScalar dot;
234: VecDot(xin, xin, &dot);
235: _z = PetscAbsScalar(PetscSqrtScalar(dot));
236: } else if (type == NORM_1) {
237: for (i = 0; i < nr; i++) {
238: VecNorm(bx->v[i], type, &z_i);
239: _z = _z + z_i;
240: }
241: } else if (type == NORM_INFINITY) {
242: for (i = 0; i < nr; i++) {
243: VecNorm(bx->v[i], type, &z_i);
244: if (z_i > _z) _z = z_i;
245: }
246: }
248: *z = _z;
249: return 0;
250: }
252: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
253: {
254: PetscInt v;
256: for (v = 0; v < nv; v++) {
257: /* Do axpy on each vector,v */
258: VecAXPY(y, alpha[v], x[v]);
259: }
260: return 0;
261: }
263: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
264: {
265: PetscInt j;
267: for (j = 0; j < nv; j++) VecDot(x, y[j], &val[j]);
268: return 0;
269: }
271: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
272: {
273: PetscInt j;
275: for (j = 0; j < nv; j++) VecTDot(x, y[j], &val[j]);
276: return 0;
277: }
279: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
280: {
281: Vec_Nest *bx = (Vec_Nest *)x->data;
282: PetscInt i, nr;
284: nr = bx->nb;
285: for (i = 0; i < nr; i++) VecSet(bx->v[i], alpha);
286: return 0;
287: }
289: static PetscErrorCode VecConjugate_Nest(Vec x)
290: {
291: Vec_Nest *bx = (Vec_Nest *)x->data;
292: PetscInt j, nr;
294: nr = bx->nb;
295: for (j = 0; j < nr; j++) VecConjugate(bx->v[j]);
296: return 0;
297: }
299: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
300: {
301: Vec_Nest *bx = (Vec_Nest *)x->data;
302: Vec_Nest *by = (Vec_Nest *)y->data;
303: PetscInt i, nr;
305: VecNestCheckCompatible2(x, 1, y, 2);
306: nr = bx->nb;
307: for (i = 0; i < nr; i++) VecSwap(bx->v[i], by->v[i]);
308: return 0;
309: }
311: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
312: {
313: Vec_Nest *bx = (Vec_Nest *)x->data;
314: Vec_Nest *by = (Vec_Nest *)y->data;
315: Vec_Nest *bw = (Vec_Nest *)w->data;
316: PetscInt i, nr;
318: VecNestCheckCompatible3(w, 1, x, 3, y, 4);
320: nr = bx->nb;
321: for (i = 0; i < nr; i++) VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]);
322: return 0;
323: }
325: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
326: {
327: PetscInt i, lp = -1, lw = -1;
328: PetscReal lv;
329: Vec_Nest *bx = (Vec_Nest *)x->data;
331: *v = PETSC_MAX_REAL;
332: for (i = 0; i < bx->nb; i++) {
333: PetscInt tp;
334: VecMin(bx->v[i], &tp, &lv);
335: if (lv < *v) {
336: *v = lv;
337: lw = i;
338: lp = tp;
339: }
340: }
341: if (p && lw > -1) {
342: PetscInt st, en;
343: const PetscInt *idxs;
345: *p = -1;
346: VecGetOwnershipRange(bx->v[lw], &st, &en);
347: if (lp >= st && lp < en) {
348: ISGetIndices(bx->is[lw], &idxs);
349: *p = idxs[lp - st];
350: ISRestoreIndices(bx->is[lw], &idxs);
351: }
352: MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x));
353: }
354: return 0;
355: }
357: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
358: {
359: PetscInt i, lp = -1, lw = -1;
360: PetscReal lv;
361: Vec_Nest *bx = (Vec_Nest *)x->data;
363: *v = PETSC_MIN_REAL;
364: for (i = 0; i < bx->nb; i++) {
365: PetscInt tp;
366: VecMax(bx->v[i], &tp, &lv);
367: if (lv > *v) {
368: *v = lv;
369: lw = i;
370: lp = tp;
371: }
372: }
373: if (p && lw > -1) {
374: PetscInt st, en;
375: const PetscInt *idxs;
377: *p = -1;
378: VecGetOwnershipRange(bx->v[lw], &st, &en);
379: if (lp >= st && lp < en) {
380: ISGetIndices(bx->is[lw], &idxs);
381: *p = idxs[lp - st];
382: ISRestoreIndices(bx->is[lw], &idxs);
383: }
384: MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x));
385: }
386: return 0;
387: }
389: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
390: {
391: Vec_Nest *bx = (Vec_Nest *)x->data;
392: PetscBool isascii;
393: PetscInt i;
395: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
396: if (isascii) {
397: PetscViewerASCIIPushTab(viewer);
398: PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ", structure: \n", bx->nb);
399: for (i = 0; i < bx->nb; i++) {
400: VecType type;
401: char name[256] = "", prefix[256] = "";
402: PetscInt NR;
404: VecGetSize(bx->v[i], &NR);
405: VecGetType(bx->v[i], &type);
406: if (((PetscObject)bx->v[i])->name) PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name);
407: if (((PetscObject)bx->v[i])->prefix) PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix);
409: PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n", i, name, prefix, type, NR);
411: PetscViewerASCIIPushTab(viewer); /* push1 */
412: VecView(bx->v[i], viewer);
413: PetscViewerASCIIPopTab(viewer); /* pop1 */
414: }
415: PetscViewerASCIIPopTab(viewer); /* pop0 */
416: }
417: return 0;
418: }
420: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
421: {
422: Vec_Nest *bx;
423: PetscInt size, i, nr;
424: PetscBool isnest;
426: PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest);
427: if (!isnest) {
428: /* Not nest */
429: if (globalsize) VecGetSize(x, &size);
430: else VecGetLocalSize(x, &size);
431: *L = *L + size;
432: return 0;
433: }
435: /* Otherwise we have a nest */
436: bx = (Vec_Nest *)x->data;
437: nr = bx->nb;
439: /* now descend recursively */
440: for (i = 0; i < nr; i++) VecSize_Nest_Recursive(bx->v[i], globalsize, L);
441: return 0;
442: }
444: /* Returns the sum of the global size of all the constituent vectors in the nest */
445: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
446: {
447: *N = x->map->N;
448: return 0;
449: }
451: /* Returns the sum of the local size of all the constituent vectors in the nest */
452: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
453: {
454: *n = x->map->n;
455: return 0;
456: }
458: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
459: {
460: Vec_Nest *bx = (Vec_Nest *)x->data;
461: Vec_Nest *by = (Vec_Nest *)y->data;
462: PetscInt i, nr;
463: PetscReal local_max, m;
465: VecNestCheckCompatible2(x, 1, y, 2);
466: nr = bx->nb;
467: m = 0.0;
468: for (i = 0; i < nr; i++) {
469: VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max);
470: if (local_max > m) m = local_max;
471: }
472: *max = m;
473: return 0;
474: }
476: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
477: {
478: Vec_Nest *bx = (Vec_Nest *)X->data;
479: PetscInt i;
481: *x = NULL;
482: for (i = 0; i < bx->nb; i++) {
483: PetscBool issame = PETSC_FALSE;
484: ISEqual(is, bx->is[i], &issame);
485: if (issame) {
486: *x = bx->v[i];
487: PetscObjectReference((PetscObject)(*x));
488: break;
489: }
490: }
492: return 0;
493: }
495: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
496: {
497: PetscObjectStateIncrease((PetscObject)X);
498: VecDestroy(x);
499: return 0;
500: }
502: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
503: {
504: Vec_Nest *bx = (Vec_Nest *)X->data;
505: PetscInt i, m, rstart, rend;
507: VecGetOwnershipRange(X, &rstart, &rend);
508: VecGetLocalSize(X, &m);
509: PetscMalloc1(m, x);
510: for (i = 0; i < bx->nb; i++) {
511: Vec subvec = bx->v[i];
512: IS isy = bx->is[i];
513: PetscInt j, sm;
514: const PetscInt *ixy;
515: const PetscScalar *y;
516: VecGetLocalSize(subvec, &sm);
517: VecGetArrayRead(subvec, &y);
518: ISGetIndices(isy, &ixy);
519: for (j = 0; j < sm; j++) {
520: PetscInt ix = ixy[j];
522: (*x)[ix - rstart] = y[j];
523: }
524: ISRestoreIndices(isy, &ixy);
525: VecRestoreArrayRead(subvec, &y);
526: }
527: return 0;
528: }
530: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
531: {
532: Vec_Nest *bx = (Vec_Nest *)X->data;
533: PetscInt i, m, rstart, rend;
535: VecGetOwnershipRange(X, &rstart, &rend);
536: VecGetLocalSize(X, &m);
537: for (i = 0; i < bx->nb; i++) {
538: Vec subvec = bx->v[i];
539: IS isy = bx->is[i];
540: PetscInt j, sm;
541: const PetscInt *ixy;
542: PetscScalar *y;
543: VecGetLocalSize(subvec, &sm);
544: VecGetArray(subvec, &y);
545: ISGetIndices(isy, &ixy);
546: for (j = 0; j < sm; j++) {
547: PetscInt ix = ixy[j];
549: y[j] = (*x)[ix - rstart];
550: }
551: ISRestoreIndices(isy, &ixy);
552: VecRestoreArray(subvec, &y);
553: }
554: PetscFree(*x);
555: PetscObjectStateIncrease((PetscObject)X);
556: return 0;
557: }
559: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
560: {
561: PetscFree(*x);
562: return 0;
563: }
565: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
566: {
568: return 0;
569: }
571: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
572: {
573: Vec *ww;
574: IS *wis;
575: Vec_Nest *bv = (Vec_Nest *)v->data;
576: PetscInt i;
578: PetscMalloc2(bv->nb, &ww, bv->nb, &wis);
579: for (i = 0; i < bv->nb; i++) VecCreateLocalVector(bv->v[i], &ww[i]);
580: for (i = 0; i < bv->nb; i++) {
581: PetscInt off;
583: VecGetOwnershipRange(v, &off, NULL);
584: ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]);
585: ISShift(wis[i], -off, wis[i]);
586: }
587: VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w);
588: for (i = 0; i < bv->nb; i++) {
589: VecDestroy(&ww[i]);
590: ISDestroy(&wis[i]);
591: }
592: PetscFree2(ww, wis);
593: return 0;
594: }
596: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
597: {
598: Vec_Nest *bv = (Vec_Nest *)v->data;
599: Vec_Nest *bw = (Vec_Nest *)w->data;
600: PetscInt i;
604: for (i = 0; i < bv->nb; i++) VecGetLocalVector(bv->v[i], bw->v[i]);
605: return 0;
606: }
608: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
609: {
610: Vec_Nest *bv = (Vec_Nest *)v->data;
611: Vec_Nest *bw = (Vec_Nest *)w->data;
612: PetscInt i;
616: for (i = 0; i < bv->nb; i++) VecRestoreLocalVector(bv->v[i], bw->v[i]);
617: PetscObjectStateIncrease((PetscObject)v);
618: return 0;
619: }
621: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
622: {
623: Vec_Nest *bv = (Vec_Nest *)v->data;
624: Vec_Nest *bw = (Vec_Nest *)w->data;
625: PetscInt i;
629: for (i = 0; i < bv->nb; i++) VecGetLocalVectorRead(bv->v[i], bw->v[i]);
630: return 0;
631: }
633: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
634: {
635: Vec_Nest *bv = (Vec_Nest *)v->data;
636: Vec_Nest *bw = (Vec_Nest *)w->data;
637: PetscInt i;
641: for (i = 0; i < bv->nb; i++) VecRestoreLocalVectorRead(bv->v[i], bw->v[i]);
642: return 0;
643: }
645: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
646: {
647: Vec_Nest *bv = (Vec_Nest *)v->data;
648: PetscInt i;
650: for (i = 0; i < bv->nb; i++) VecSetRandom(bv->v[i], r);
651: return 0;
652: }
654: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
655: {
656: ops->duplicate = VecDuplicate_Nest;
657: ops->duplicatevecs = VecDuplicateVecs_Default;
658: ops->destroyvecs = VecDestroyVecs_Default;
659: ops->dot = VecDot_Nest;
660: ops->mdot = VecMDot_Nest;
661: ops->norm = VecNorm_Nest;
662: ops->tdot = VecTDot_Nest;
663: ops->mtdot = VecMTDot_Nest;
664: ops->scale = VecScale_Nest;
665: ops->copy = VecCopy_Nest;
666: ops->set = VecSet_Nest;
667: ops->swap = VecSwap_Nest;
668: ops->axpy = VecAXPY_Nest;
669: ops->axpby = VecAXPBY_Nest;
670: ops->maxpy = VecMAXPY_Nest;
671: ops->aypx = VecAYPX_Nest;
672: ops->waxpy = VecWAXPY_Nest;
673: ops->axpbypcz = NULL;
674: ops->pointwisemult = VecPointwiseMult_Nest;
675: ops->pointwisedivide = VecPointwiseDivide_Nest;
676: ops->setvalues = NULL;
677: ops->assemblybegin = VecAssemblyBegin_Nest;
678: ops->assemblyend = VecAssemblyEnd_Nest;
679: ops->getarray = VecGetArray_Nest;
680: ops->getsize = VecGetSize_Nest;
681: ops->getlocalsize = VecGetLocalSize_Nest;
682: ops->restorearray = VecRestoreArray_Nest;
683: ops->restorearrayread = VecRestoreArrayRead_Nest;
684: ops->max = VecMax_Nest;
685: ops->min = VecMin_Nest;
686: ops->setrandom = NULL;
687: ops->setoption = NULL;
688: ops->setvaluesblocked = NULL;
689: ops->destroy = VecDestroy_Nest;
690: ops->view = VecView_Nest;
691: ops->placearray = NULL;
692: ops->replacearray = NULL;
693: ops->dot_local = VecDot_Nest;
694: ops->tdot_local = VecTDot_Nest;
695: ops->norm_local = VecNorm_Nest;
696: ops->mdot_local = VecMDot_Nest;
697: ops->mtdot_local = VecMTDot_Nest;
698: ops->load = NULL;
699: ops->reciprocal = VecReciprocal_Nest;
700: ops->conjugate = VecConjugate_Nest;
701: ops->setlocaltoglobalmapping = NULL;
702: ops->setvalueslocal = NULL;
703: ops->resetarray = NULL;
704: ops->setfromoptions = NULL;
705: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
706: ops->load = NULL;
707: ops->pointwisemax = NULL;
708: ops->pointwisemaxabs = NULL;
709: ops->pointwisemin = NULL;
710: ops->getvalues = NULL;
711: ops->sqrt = NULL;
712: ops->abs = NULL;
713: ops->exp = NULL;
714: ops->shift = NULL;
715: ops->create = NULL;
716: ops->stridegather = NULL;
717: ops->stridescatter = NULL;
718: ops->dotnorm2 = VecDotNorm2_Nest;
719: ops->getsubvector = VecGetSubVector_Nest;
720: ops->restoresubvector = VecRestoreSubVector_Nest;
721: ops->axpbypcz = VecAXPBYPCZ_Nest;
722: ops->concatenate = VecConcatenate_Nest;
723: ops->createlocalvector = VecCreateLocalVector_Nest;
724: ops->getlocalvector = VecGetLocalVector_Nest;
725: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
726: ops->restorelocalvector = VecRestoreLocalVector_Nest;
727: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
728: ops->setrandom = VecSetRandom_Nest;
729: return 0;
730: }
732: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
733: {
734: Vec_Nest *b = (Vec_Nest *)x->data;
735: PetscInt i;
736: PetscInt row;
738: if (!m) return 0;
739: for (i = 0; i < m; i++) {
740: row = idxm[i];
742: vec[i] = b->v[row];
743: }
744: return 0;
745: }
747: PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
748: {
749: PetscObjectStateIncrease((PetscObject)X);
750: VecNestGetSubVecs_Private(X, 1, &idxm, sx);
751: return 0;
752: }
754: /*@
755: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
757: Not collective
759: Input Parameters:
760: + X - nest vector
761: - idxm - index of the vector within the nest
763: Output Parameter:
764: . sx - vector at index idxm within the nest
766: Notes:
768: Level: developer
770: .seealso: `VecNestGetSize()`, `VecNestGetSubVecs()`
771: @*/
772: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
773: {
774: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
775: return 0;
776: }
778: PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
779: {
780: Vec_Nest *b = (Vec_Nest *)X->data;
782: PetscObjectStateIncrease((PetscObject)X);
783: if (N) *N = b->nb;
784: if (sx) *sx = b->v;
785: return 0;
786: }
788: /*@C
789: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
791: Not collective
793: Input Parameter:
794: . X - nest vector
796: Output Parameters:
797: + N - number of nested vecs
798: - sx - array of vectors
800: Notes:
801: The user should not free the array sx.
803: Fortran Notes:
804: The caller must allocate the array to hold the subvectors.
806: Level: developer
808: .seealso: `VecNestGetSize()`, `VecNestGetSubVec()`
809: @*/
810: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec **sx)
811: {
812: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
813: return 0;
814: }
816: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
817: {
818: Vec_Nest *bx = (Vec_Nest *)X->data;
819: PetscInt i, offset = 0, n = 0, bs;
820: IS is;
821: PetscBool issame = PETSC_FALSE;
822: PetscInt N = 0;
824: /* check if idxm < bx->nb */
827: VecDestroy(&bx->v[idxm]); /* destroy the existing vector */
828: VecDuplicate(x, &bx->v[idxm]); /* duplicate the layout of given vector */
829: VecCopy(x, bx->v[idxm]); /* copy the contents of the given vector */
831: /* check if we need to update the IS for the block */
832: offset = X->map->rstart;
833: for (i = 0; i < idxm; i++) {
834: n = 0;
835: VecGetLocalSize(bx->v[i], &n);
836: offset += n;
837: }
839: /* get the local size and block size */
840: VecGetLocalSize(x, &n);
841: VecGetBlockSize(x, &bs);
843: /* create the new IS */
844: ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is);
845: ISSetBlockSize(is, bs);
847: /* check if they are equal */
848: ISEqual(is, bx->is[idxm], &issame);
850: if (!issame) {
851: /* The IS of given vector has a different layout compared to the existing block vector.
852: Destroy the existing reference and update the IS. */
853: ISDestroy(&bx->is[idxm]);
854: ISDuplicate(is, &bx->is[idxm]);
855: ISCopy(is, bx->is[idxm]);
857: offset += n;
858: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
859: for (i = idxm + 1; i < bx->nb; i++) {
860: /* get the local size and block size */
861: VecGetLocalSize(bx->v[i], &n);
862: VecGetBlockSize(bx->v[i], &bs);
864: /* destroy the old and create the new IS */
865: ISDestroy(&bx->is[i]);
866: ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]);
867: ISSetBlockSize(bx->is[i], bs);
869: offset += n;
870: }
872: n = 0;
873: VecSize_Nest_Recursive(X, PETSC_TRUE, &N);
874: VecSize_Nest_Recursive(X, PETSC_FALSE, &n);
875: PetscLayoutSetSize(X->map, N);
876: PetscLayoutSetLocalSize(X->map, n);
877: }
879: ISDestroy(&is);
880: return 0;
881: }
883: PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
884: {
885: PetscObjectStateIncrease((PetscObject)X);
886: VecNestSetSubVec_Private(X, idxm, sx);
887: return 0;
888: }
890: /*@
891: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
893: Not collective
895: Input Parameters:
896: + X - nest vector
897: . idxm - index of the vector within the nest vector
898: - sx - vector at index idxm within the nest vector
900: Notes:
901: The new vector sx does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
903: Level: developer
905: .seealso: `VecNestSetSubVecs()`, `VecNestGetSubVec()`
906: @*/
907: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
908: {
909: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
910: return 0;
911: }
913: PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
914: {
915: PetscInt i;
917: PetscObjectStateIncrease((PetscObject)X);
918: for (i = 0; i < N; i++) VecNestSetSubVec_Private(X, idxm[i], sx[i]);
919: return 0;
920: }
922: /*@C
923: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
925: Not collective
927: Input Parameters:
928: + X - nest vector
929: . N - number of component vecs in sx
930: . idxm - indices of component vecs that are to be replaced
931: - sx - array of vectors
933: Notes:
934: The components in the vector array sx do not have to be of the same size as corresponding
935: components in X. The user can also free the array "sx" after the call.
937: Level: developer
939: .seealso: `VecNestGetSize()`, `VecNestGetSubVec()`
940: @*/
941: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
942: {
943: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
944: return 0;
945: }
947: PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
948: {
949: Vec_Nest *b = (Vec_Nest *)X->data;
951: *N = b->nb;
952: return 0;
953: }
955: /*@
956: VecNestGetSize - Returns the size of the nest vector.
958: Not collective
960: Input Parameter:
961: . X - nest vector
963: Output Parameter:
964: . N - number of nested vecs
966: Notes:
968: Level: developer
970: .seealso: `VecNestGetSubVec()`, `VecNestGetSubVecs()`
971: @*/
972: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
973: {
976: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
977: return 0;
978: }
980: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
981: {
982: Vec_Nest *ctx = (Vec_Nest *)V->data;
983: PetscInt i;
985: if (ctx->setup_called) return 0;
987: ctx->nb = nb;
990: /* Create space */
991: PetscMalloc1(ctx->nb, &ctx->v);
992: for (i = 0; i < ctx->nb; i++) {
993: ctx->v[i] = x[i];
994: PetscObjectReference((PetscObject)x[i]);
995: /* Do not allocate memory for internal sub blocks */
996: }
998: PetscMalloc1(ctx->nb, &ctx->is);
1000: ctx->setup_called = PETSC_TRUE;
1001: return 0;
1002: }
1004: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1005: {
1006: Vec_Nest *ctx = (Vec_Nest *)V->data;
1007: PetscInt i, offset, m, n, M, N;
1009: if (is) { /* Do some consistency checks and reference the is */
1010: offset = V->map->rstart;
1011: for (i = 0; i < ctx->nb; i++) {
1012: ISGetSize(is[i], &M);
1013: VecGetSize(ctx->v[i], &N);
1015: ISGetLocalSize(is[i], &m);
1016: VecGetLocalSize(ctx->v[i], &n);
1018: PetscObjectReference((PetscObject)is[i]);
1019: ctx->is[i] = is[i];
1020: offset += n;
1021: }
1022: } else { /* Create a contiguous ISStride for each entry */
1023: offset = V->map->rstart;
1024: for (i = 0; i < ctx->nb; i++) {
1025: PetscInt bs;
1026: VecGetLocalSize(ctx->v[i], &n);
1027: VecGetBlockSize(ctx->v[i], &bs);
1028: ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]);
1029: ISSetBlockSize(ctx->is[i], bs);
1030: offset += n;
1031: }
1032: }
1033: return 0;
1034: }
1036: /*@C
1037: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1039: Collective
1041: Input Parameters:
1042: + comm - Communicator for the new `Vec`
1043: . nb - number of nested blocks
1044: . is - array of nb index sets describing each nested block, or NULL to pack subvectors contiguously
1045: - x - array of nb sub-vectors
1047: Output Parameter:
1048: . Y - new vector
1050: Level: advanced
1052: .seealso: `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`, `VECNEST`
1053: @*/
1054: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1055: {
1056: Vec V;
1057: Vec_Nest *s;
1058: PetscInt n, N;
1060: VecCreate(comm, &V);
1061: PetscObjectChangeTypeName((PetscObject)V, VECNEST);
1063: /* allocate and set pointer for implementation data */
1064: PetscNew(&s);
1065: V->data = (void *)s;
1066: s->setup_called = PETSC_FALSE;
1067: s->nb = -1;
1068: s->v = NULL;
1070: VecSetUp_Nest_Private(V, nb, x);
1072: n = N = 0;
1073: VecSize_Nest_Recursive(V, PETSC_TRUE, &N);
1074: VecSize_Nest_Recursive(V, PETSC_FALSE, &n);
1075: PetscLayoutSetSize(V->map, N);
1076: PetscLayoutSetLocalSize(V->map, n);
1077: PetscLayoutSetBlockSize(V->map, 1);
1078: PetscLayoutSetUp(V->map);
1080: VecSetUp_NestIS_Private(V, nb, is);
1082: VecNestSetOps_Private(V->ops);
1083: V->petscnative = PETSC_FALSE;
1085: /* expose block api's */
1086: PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest);
1087: PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest);
1088: PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest);
1089: PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest);
1090: PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest);
1092: *Y = V;
1093: return 0;
1094: }
1096: /*MC
1097: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1099: Level: intermediate
1101: Notes:
1102: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1103: It is usually used with MATNEST and DMComposite via DMSetVecType().
1105: .seealso: `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1106: M*/