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: PetscFunctionBegin;
11: for (i = 0; i < vs->nb; i++) {
12: PetscCheck(vs->v[i], PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Nest vector cannot contain NULL blocks");
13: PetscCall(VecAssemblyBegin(vs->v[i]));
14: }
15: PetscFunctionReturn(PETSC_SUCCESS);
16: }
18: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
19: {
20: Vec_Nest *vs = (Vec_Nest *)v->data;
21: PetscInt i;
23: PetscFunctionBegin;
24: for (i = 0; i < vs->nb; i++) PetscCall(VecAssemblyEnd(vs->v[i]));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode VecDestroy_Nest(Vec v)
29: {
30: Vec_Nest *vs = (Vec_Nest *)v->data;
31: PetscInt i;
33: PetscFunctionBegin;
34: if (vs->v) {
35: for (i = 0; i < vs->nb; i++) PetscCall(VecDestroy(&vs->v[i]));
36: PetscCall(PetscFree(vs->v));
37: }
38: for (i = 0; i < vs->nb; i++) PetscCall(ISDestroy(&vs->is[i]));
39: PetscCall(PetscFree(vs->is));
40: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVec_C", NULL));
41: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecs_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVec_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL));
44: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL));
46: PetscCall(PetscFree(vs));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
51: {
52: Vec_Nest *bx = (Vec_Nest *)x->data;
53: Vec_Nest *by = (Vec_Nest *)y->data;
54: PetscInt i;
56: PetscFunctionBegin;
57: VecNestCheckCompatible2(x, 1, y, 2);
58: for (i = 0; i < bx->nb; i++) PetscCall(VecCopy(bx->v[i], by->v[i]));
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
63: {
64: Vec_Nest *bx = (Vec_Nest *)x->data;
65: Vec Y;
66: Vec *sub;
67: PetscInt i;
69: PetscFunctionBegin;
70: PetscCall(PetscMalloc1(bx->nb, &sub));
71: for (i = 0; i < bx->nb; i++) PetscCall(VecDuplicate(bx->v[i], &sub[i]));
72: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y));
73: for (i = 0; i < bx->nb; i++) PetscCall(VecDestroy(&sub[i]));
74: PetscCall(PetscFree(sub));
75: *y = Y;
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
80: {
81: Vec_Nest *bx = (Vec_Nest *)x->data;
82: Vec_Nest *by = (Vec_Nest *)y->data;
83: PetscInt i, nr;
84: PetscScalar x_dot_y, _val;
86: PetscFunctionBegin;
87: VecNestCheckCompatible2(x, 1, y, 2);
88: nr = bx->nb;
89: _val = 0.0;
90: for (i = 0; i < nr; i++) {
91: PetscCall(VecDot(bx->v[i], by->v[i], &x_dot_y));
92: _val = _val + x_dot_y;
93: }
94: *val = _val;
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
99: {
100: Vec_Nest *bx = (Vec_Nest *)x->data;
101: Vec_Nest *by = (Vec_Nest *)y->data;
102: PetscInt i, nr;
103: PetscScalar x_dot_y, _val;
105: PetscFunctionBegin;
106: VecNestCheckCompatible2(x, 1, y, 2);
107: nr = bx->nb;
108: _val = 0.0;
109: for (i = 0; i < nr; i++) {
110: PetscCall(VecTDot(bx->v[i], by->v[i], &x_dot_y));
111: _val = _val + x_dot_y;
112: }
113: *val = _val;
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
118: {
119: Vec_Nest *bx = (Vec_Nest *)x->data;
120: Vec_Nest *by = (Vec_Nest *)y->data;
121: PetscInt i, nr;
122: PetscScalar x_dot_y, _dp, _nm;
123: PetscReal norm2_y;
125: PetscFunctionBegin;
126: VecNestCheckCompatible2(x, 1, y, 2);
127: nr = bx->nb;
128: _dp = 0.0;
129: _nm = 0.0;
130: for (i = 0; i < nr; i++) {
131: PetscCall(VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y));
132: _dp += x_dot_y;
133: _nm += norm2_y;
134: }
135: *dp = _dp;
136: *nm = _nm;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
141: {
142: Vec_Nest *bx = (Vec_Nest *)x->data;
143: Vec_Nest *by = (Vec_Nest *)y->data;
144: PetscInt i, nr;
146: PetscFunctionBegin;
147: VecNestCheckCompatible2(y, 1, x, 3);
148: nr = bx->nb;
149: for (i = 0; i < nr; i++) PetscCall(VecAXPY(by->v[i], alpha, bx->v[i]));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
154: {
155: Vec_Nest *bx = (Vec_Nest *)x->data;
156: Vec_Nest *by = (Vec_Nest *)y->data;
157: PetscInt i, nr;
159: PetscFunctionBegin;
160: VecNestCheckCompatible2(y, 1, x, 3);
161: nr = bx->nb;
162: for (i = 0; i < nr; i++) PetscCall(VecAYPX(by->v[i], alpha, bx->v[i]));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
167: {
168: Vec_Nest *bx = (Vec_Nest *)x->data;
169: Vec_Nest *by = (Vec_Nest *)y->data;
170: PetscInt i, nr;
172: PetscFunctionBegin;
173: VecNestCheckCompatible2(y, 1, x, 4);
174: nr = bx->nb;
175: for (i = 0; i < nr; i++) PetscCall(VecAXPBY(by->v[i], alpha, beta, bx->v[i]));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
180: {
181: Vec_Nest *bx = (Vec_Nest *)x->data;
182: Vec_Nest *by = (Vec_Nest *)y->data;
183: Vec_Nest *bz = (Vec_Nest *)z->data;
184: PetscInt i, nr;
186: PetscFunctionBegin;
187: VecNestCheckCompatible3(z, 1, x, 5, y, 6);
188: nr = bx->nb;
189: for (i = 0; i < nr; i++) PetscCall(VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
194: {
195: Vec_Nest *bx = (Vec_Nest *)x->data;
196: PetscInt i, nr;
198: PetscFunctionBegin;
199: nr = bx->nb;
200: for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
205: {
206: Vec_Nest *bx = (Vec_Nest *)x->data;
207: Vec_Nest *by = (Vec_Nest *)y->data;
208: Vec_Nest *bw = (Vec_Nest *)w->data;
209: PetscInt i, nr;
211: PetscFunctionBegin;
212: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
213: nr = bx->nb;
214: for (i = 0; i < nr; i++) PetscCall(VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
219: {
220: Vec_Nest *bx = (Vec_Nest *)x->data;
221: Vec_Nest *by = (Vec_Nest *)y->data;
222: Vec_Nest *bw = (Vec_Nest *)w->data;
223: PetscInt i, nr;
225: PetscFunctionBegin;
226: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
227: nr = bx->nb;
228: for (i = 0; i < nr; i++) PetscCall(VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode VecReciprocal_Nest(Vec x)
233: {
234: Vec_Nest *bx = (Vec_Nest *)x->data;
235: PetscInt i, nr;
237: PetscFunctionBegin;
238: nr = bx->nb;
239: for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
244: {
245: Vec_Nest *bx = (Vec_Nest *)xin->data;
246: PetscInt i, nr;
247: PetscReal z_i;
248: PetscReal _z;
250: PetscFunctionBegin;
251: nr = bx->nb;
252: _z = 0.0;
254: if (type == NORM_2) {
255: PetscScalar dot;
256: PetscCall(VecDot(xin, xin, &dot));
257: _z = PetscAbsScalar(PetscSqrtScalar(dot));
258: } else if (type == NORM_1) {
259: for (i = 0; i < nr; i++) {
260: PetscCall(VecNorm(bx->v[i], type, &z_i));
261: _z = _z + z_i;
262: }
263: } else if (type == NORM_INFINITY) {
264: for (i = 0; i < nr; i++) {
265: PetscCall(VecNorm(bx->v[i], type, &z_i));
266: if (z_i > _z) _z = z_i;
267: }
268: }
270: *z = _z;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
275: {
276: PetscFunctionBegin;
277: /* TODO: implement proper MAXPY. For now, do axpy on each vector */
278: for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
283: {
284: PetscFunctionBegin;
285: /* TODO: implement proper MDOT. For now, do dot on each vector */
286: for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
291: {
292: PetscFunctionBegin;
293: /* TODO: implement proper MTDOT. For now, do tdot on each vector */
294: for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
299: {
300: Vec_Nest *bx = (Vec_Nest *)x->data;
301: PetscInt i, nr;
303: PetscFunctionBegin;
304: nr = bx->nb;
305: for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode VecConjugate_Nest(Vec x)
310: {
311: Vec_Nest *bx = (Vec_Nest *)x->data;
312: PetscInt j, nr;
314: PetscFunctionBegin;
315: nr = bx->nb;
316: for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
321: {
322: Vec_Nest *bx = (Vec_Nest *)x->data;
323: Vec_Nest *by = (Vec_Nest *)y->data;
324: PetscInt i, nr;
326: PetscFunctionBegin;
327: VecNestCheckCompatible2(x, 1, y, 2);
328: nr = bx->nb;
329: for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
334: {
335: Vec_Nest *bx = (Vec_Nest *)x->data;
336: Vec_Nest *by = (Vec_Nest *)y->data;
337: Vec_Nest *bw = (Vec_Nest *)w->data;
338: PetscInt i, nr;
340: PetscFunctionBegin;
341: VecNestCheckCompatible3(w, 1, x, 3, y, 4);
342: nr = bx->nb;
343: for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
348: {
349: PetscInt i, lp = -1, lw = -1;
350: PetscReal lv;
351: Vec_Nest *bx = (Vec_Nest *)x->data;
353: PetscFunctionBegin;
354: *v = PETSC_MAX_REAL;
355: for (i = 0; i < bx->nb; i++) {
356: PetscInt tp;
357: PetscCall(VecMin(bx->v[i], &tp, &lv));
358: if (lv < *v) {
359: *v = lv;
360: lw = i;
361: lp = tp;
362: }
363: }
364: if (p && lw > -1) {
365: PetscInt st, en;
366: const PetscInt *idxs;
368: *p = -1;
369: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
370: if (lp >= st && lp < en) {
371: PetscCall(ISGetIndices(bx->is[lw], &idxs));
372: *p = idxs[lp - st];
373: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
374: }
375: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
376: }
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
381: {
382: PetscInt i, lp = -1, lw = -1;
383: PetscReal lv;
384: Vec_Nest *bx = (Vec_Nest *)x->data;
386: PetscFunctionBegin;
387: *v = PETSC_MIN_REAL;
388: for (i = 0; i < bx->nb; i++) {
389: PetscInt tp;
390: PetscCall(VecMax(bx->v[i], &tp, &lv));
391: if (lv > *v) {
392: *v = lv;
393: lw = i;
394: lp = tp;
395: }
396: }
397: if (p && lw > -1) {
398: PetscInt st, en;
399: const PetscInt *idxs;
401: *p = -1;
402: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
403: if (lp >= st && lp < en) {
404: PetscCall(ISGetIndices(bx->is[lw], &idxs));
405: *p = idxs[lp - st];
406: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
407: }
408: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
409: }
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
414: {
415: Vec_Nest *bx = (Vec_Nest *)x->data;
416: PetscBool isascii;
417: PetscInt i;
419: PetscFunctionBegin;
420: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
421: if (isascii) {
422: PetscCall(PetscViewerASCIIPushTab(viewer));
423: PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ", structure: \n", bx->nb));
424: for (i = 0; i < bx->nb; i++) {
425: VecType type;
426: char name[256] = "", prefix[256] = "";
427: PetscInt NR;
429: PetscCall(VecGetSize(bx->v[i], &NR));
430: PetscCall(VecGetType(bx->v[i], &type));
431: if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
432: if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));
434: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n", i, name, prefix, type, NR));
436: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
437: PetscCall(VecView(bx->v[i], viewer));
438: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
439: }
440: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
446: {
447: Vec_Nest *bx;
448: PetscInt size, i, nr;
449: PetscBool isnest;
451: PetscFunctionBegin;
452: PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
453: if (!isnest) {
454: /* Not nest */
455: if (globalsize) PetscCall(VecGetSize(x, &size));
456: else PetscCall(VecGetLocalSize(x, &size));
457: *L = *L + size;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /* Otherwise we have a nest */
462: bx = (Vec_Nest *)x->data;
463: nr = bx->nb;
465: /* now descend recursively */
466: for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /* Returns the sum of the global size of all the constituent vectors in the nest */
471: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
472: {
473: PetscFunctionBegin;
474: *N = x->map->N;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: /* Returns the sum of the local size of all the constituent vectors in the nest */
479: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
480: {
481: PetscFunctionBegin;
482: *n = x->map->n;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
487: {
488: Vec_Nest *bx = (Vec_Nest *)x->data;
489: Vec_Nest *by = (Vec_Nest *)y->data;
490: PetscInt i, nr;
491: PetscReal local_max, m;
493: PetscFunctionBegin;
494: VecNestCheckCompatible2(x, 1, y, 2);
495: nr = bx->nb;
496: m = 0.0;
497: for (i = 0; i < nr; i++) {
498: PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
499: if (local_max > m) m = local_max;
500: }
501: *max = m;
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
506: {
507: Vec_Nest *bx = (Vec_Nest *)X->data;
508: PetscInt i;
510: PetscFunctionBegin;
511: *x = NULL;
512: for (i = 0; i < bx->nb; i++) {
513: PetscBool issame = PETSC_FALSE;
514: PetscCall(ISEqual(is, bx->is[i], &issame));
515: if (issame) {
516: *x = bx->v[i];
517: PetscCall(PetscObjectReference((PetscObject)(*x)));
518: break;
519: }
520: }
521: PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
526: {
527: PetscFunctionBegin;
528: PetscCall(PetscObjectStateIncrease((PetscObject)X));
529: PetscCall(VecDestroy(x));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
534: {
535: Vec_Nest *bx = (Vec_Nest *)X->data;
536: PetscInt i, m, rstart, rend;
538: PetscFunctionBegin;
539: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
540: PetscCall(VecGetLocalSize(X, &m));
541: PetscCall(PetscMalloc1(m, x));
542: for (i = 0; i < bx->nb; i++) {
543: Vec subvec = bx->v[i];
544: IS isy = bx->is[i];
545: PetscInt j, sm;
546: const PetscInt *ixy;
547: const PetscScalar *y;
548: PetscCall(VecGetLocalSize(subvec, &sm));
549: PetscCall(VecGetArrayRead(subvec, &y));
550: PetscCall(ISGetIndices(isy, &ixy));
551: for (j = 0; j < sm; j++) {
552: PetscInt ix = ixy[j];
553: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
554: (*x)[ix - rstart] = y[j];
555: }
556: PetscCall(ISRestoreIndices(isy, &ixy));
557: PetscCall(VecRestoreArrayRead(subvec, &y));
558: }
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
563: {
564: Vec_Nest *bx = (Vec_Nest *)X->data;
565: PetscInt i, m, rstart, rend;
567: PetscFunctionBegin;
568: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
569: PetscCall(VecGetLocalSize(X, &m));
570: for (i = 0; i < bx->nb; i++) {
571: Vec subvec = bx->v[i];
572: IS isy = bx->is[i];
573: PetscInt j, sm;
574: const PetscInt *ixy;
575: PetscScalar *y;
576: PetscCall(VecGetLocalSize(subvec, &sm));
577: PetscCall(VecGetArray(subvec, &y));
578: PetscCall(ISGetIndices(isy, &ixy));
579: for (j = 0; j < sm; j++) {
580: PetscInt ix = ixy[j];
581: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
582: y[j] = (*x)[ix - rstart];
583: }
584: PetscCall(ISRestoreIndices(isy, &ixy));
585: PetscCall(VecRestoreArray(subvec, &y));
586: }
587: PetscCall(PetscFree(*x));
588: PetscCall(PetscObjectStateIncrease((PetscObject)X));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
593: {
594: PetscFunctionBegin;
595: PetscCall(PetscFree(*x));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
600: {
601: PetscFunctionBegin;
602: PetscCheck(nx <= 0, PetscObjectComm((PetscObject)(*X)), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
607: {
608: Vec *ww;
609: IS *wis;
610: Vec_Nest *bv = (Vec_Nest *)v->data;
611: PetscInt i;
613: PetscFunctionBegin;
614: PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
615: for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
616: for (i = 0; i < bv->nb; i++) {
617: PetscInt off;
619: PetscCall(VecGetOwnershipRange(v, &off, NULL));
620: PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
621: PetscCall(ISShift(wis[i], -off, wis[i]));
622: }
623: PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
624: for (i = 0; i < bv->nb; i++) {
625: PetscCall(VecDestroy(&ww[i]));
626: PetscCall(ISDestroy(&wis[i]));
627: }
628: PetscCall(PetscFree2(ww, wis));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
633: {
634: Vec_Nest *bv = (Vec_Nest *)v->data;
635: Vec_Nest *bw = (Vec_Nest *)w->data;
636: PetscInt i;
638: PetscFunctionBegin;
639: PetscCheckSameType(v, 1, w, 2);
640: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
641: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
646: {
647: Vec_Nest *bv = (Vec_Nest *)v->data;
648: Vec_Nest *bw = (Vec_Nest *)w->data;
649: PetscInt i;
651: PetscFunctionBegin;
652: PetscCheckSameType(v, 1, w, 2);
653: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
654: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
655: PetscCall(PetscObjectStateIncrease((PetscObject)v));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
660: {
661: Vec_Nest *bv = (Vec_Nest *)v->data;
662: Vec_Nest *bw = (Vec_Nest *)w->data;
663: PetscInt i;
665: PetscFunctionBegin;
666: PetscCheckSameType(v, 1, w, 2);
667: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
668: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
673: {
674: Vec_Nest *bv = (Vec_Nest *)v->data;
675: Vec_Nest *bw = (Vec_Nest *)w->data;
676: PetscInt i;
678: PetscFunctionBegin;
679: PetscCheckSameType(v, 1, w, 2);
680: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
681: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
686: {
687: Vec_Nest *bv = (Vec_Nest *)v->data;
688: PetscInt i;
690: PetscFunctionBegin;
691: for (i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
696: {
697: PetscFunctionBegin;
698: ops->duplicate = VecDuplicate_Nest;
699: ops->duplicatevecs = VecDuplicateVecs_Default;
700: ops->destroyvecs = VecDestroyVecs_Default;
701: ops->dot = VecDot_Nest;
702: ops->mdot = VecMDot_Nest;
703: ops->norm = VecNorm_Nest;
704: ops->tdot = VecTDot_Nest;
705: ops->mtdot = VecMTDot_Nest;
706: ops->scale = VecScale_Nest;
707: ops->copy = VecCopy_Nest;
708: ops->set = VecSet_Nest;
709: ops->swap = VecSwap_Nest;
710: ops->axpy = VecAXPY_Nest;
711: ops->axpby = VecAXPBY_Nest;
712: ops->maxpy = VecMAXPY_Nest;
713: ops->aypx = VecAYPX_Nest;
714: ops->waxpy = VecWAXPY_Nest;
715: ops->axpbypcz = NULL;
716: ops->pointwisemult = VecPointwiseMult_Nest;
717: ops->pointwisedivide = VecPointwiseDivide_Nest;
718: ops->setvalues = NULL;
719: ops->assemblybegin = VecAssemblyBegin_Nest;
720: ops->assemblyend = VecAssemblyEnd_Nest;
721: ops->getarray = VecGetArray_Nest;
722: ops->getsize = VecGetSize_Nest;
723: ops->getlocalsize = VecGetLocalSize_Nest;
724: ops->restorearray = VecRestoreArray_Nest;
725: ops->restorearrayread = VecRestoreArrayRead_Nest;
726: ops->max = VecMax_Nest;
727: ops->min = VecMin_Nest;
728: ops->setrandom = NULL;
729: ops->setoption = NULL;
730: ops->setvaluesblocked = NULL;
731: ops->destroy = VecDestroy_Nest;
732: ops->view = VecView_Nest;
733: ops->placearray = NULL;
734: ops->replacearray = NULL;
735: ops->dot_local = VecDot_Nest;
736: ops->tdot_local = VecTDot_Nest;
737: ops->norm_local = VecNorm_Nest;
738: ops->mdot_local = VecMDot_Nest;
739: ops->mtdot_local = VecMTDot_Nest;
740: ops->load = NULL;
741: ops->reciprocal = VecReciprocal_Nest;
742: ops->conjugate = VecConjugate_Nest;
743: ops->setlocaltoglobalmapping = NULL;
744: ops->setvalueslocal = NULL;
745: ops->resetarray = NULL;
746: ops->setfromoptions = NULL;
747: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
748: ops->load = NULL;
749: ops->pointwisemax = NULL;
750: ops->pointwisemaxabs = NULL;
751: ops->pointwisemin = NULL;
752: ops->getvalues = NULL;
753: ops->sqrt = NULL;
754: ops->abs = NULL;
755: ops->exp = NULL;
756: ops->shift = NULL;
757: ops->create = NULL;
758: ops->stridegather = NULL;
759: ops->stridescatter = NULL;
760: ops->dotnorm2 = VecDotNorm2_Nest;
761: ops->getsubvector = VecGetSubVector_Nest;
762: ops->restoresubvector = VecRestoreSubVector_Nest;
763: ops->axpbypcz = VecAXPBYPCZ_Nest;
764: ops->concatenate = VecConcatenate_Nest;
765: ops->createlocalvector = VecCreateLocalVector_Nest;
766: ops->getlocalvector = VecGetLocalVector_Nest;
767: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
768: ops->restorelocalvector = VecRestoreLocalVector_Nest;
769: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
770: ops->setrandom = VecSetRandom_Nest;
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
775: {
776: Vec_Nest *b = (Vec_Nest *)x->data;
777: PetscInt i;
778: PetscInt row;
780: PetscFunctionBegin;
781: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
782: for (i = 0; i < m; i++) {
783: row = idxm[i];
784: PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
785: vec[i] = b->v[row];
786: }
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
791: {
792: PetscFunctionBegin;
793: PetscCall(PetscObjectStateIncrease((PetscObject)X));
794: PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
801: Not Collective
803: Input Parameters:
804: + X - nest vector
805: - idxm - index of the vector within the nest
807: Output Parameter:
808: . sx - vector at index `idxm` within the nest
810: Level: developer
812: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
813: @*/
814: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
815: {
816: PetscFunctionBegin;
817: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
822: {
823: Vec_Nest *b = (Vec_Nest *)X->data;
825: PetscFunctionBegin;
826: PetscCall(PetscObjectStateIncrease((PetscObject)X));
827: if (N) *N = b->nb;
828: if (sx) *sx = b->v;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*@C
833: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
835: Not Collective
837: Input Parameter:
838: . X - nest vector
840: Output Parameters:
841: + N - number of nested vecs
842: - sx - array of vectors
844: Level: developer
846: Note:
847: The user should not free the array `sx`.
849: Fortran Note:
850: The caller must allocate the array to hold the subvectors.
852: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
853: @*/
854: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec **sx)
855: {
856: PetscFunctionBegin;
857: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
858: PetscFunctionReturn(PETSC_SUCCESS);
859: }
861: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
862: {
863: Vec_Nest *bx = (Vec_Nest *)X->data;
864: PetscInt i, offset = 0, n = 0, bs;
865: IS is;
866: PetscBool issame = PETSC_FALSE;
867: PetscInt N = 0;
869: PetscFunctionBegin;
870: /* check if idxm < bx->nb */
871: PetscCheck(idxm < bx->nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT " maximum %" PetscInt_FMT, idxm, bx->nb);
873: PetscCall(VecDestroy(&bx->v[idxm])); /* destroy the existing vector */
874: PetscCall(VecDuplicate(x, &bx->v[idxm])); /* duplicate the layout of given vector */
875: PetscCall(VecCopy(x, bx->v[idxm])); /* copy the contents of the given vector */
877: /* check if we need to update the IS for the block */
878: offset = X->map->rstart;
879: for (i = 0; i < idxm; i++) {
880: n = 0;
881: PetscCall(VecGetLocalSize(bx->v[i], &n));
882: offset += n;
883: }
885: /* get the local size and block size */
886: PetscCall(VecGetLocalSize(x, &n));
887: PetscCall(VecGetBlockSize(x, &bs));
889: /* create the new IS */
890: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
891: PetscCall(ISSetBlockSize(is, bs));
893: /* check if they are equal */
894: PetscCall(ISEqual(is, bx->is[idxm], &issame));
896: if (!issame) {
897: /* The IS of given vector has a different layout compared to the existing block vector.
898: Destroy the existing reference and update the IS. */
899: PetscCall(ISDestroy(&bx->is[idxm]));
900: PetscCall(ISDuplicate(is, &bx->is[idxm]));
901: PetscCall(ISCopy(is, bx->is[idxm]));
903: offset += n;
904: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
905: for (i = idxm + 1; i < bx->nb; i++) {
906: /* get the local size and block size */
907: PetscCall(VecGetLocalSize(bx->v[i], &n));
908: PetscCall(VecGetBlockSize(bx->v[i], &bs));
910: /* destroy the old and create the new IS */
911: PetscCall(ISDestroy(&bx->is[i]));
912: PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
913: PetscCall(ISSetBlockSize(bx->is[i], bs));
915: offset += n;
916: }
918: n = 0;
919: PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
920: PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
921: PetscCall(PetscLayoutSetSize(X->map, N));
922: PetscCall(PetscLayoutSetLocalSize(X->map, n));
923: }
925: PetscCall(ISDestroy(&is));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
930: {
931: PetscFunctionBegin;
932: PetscCall(PetscObjectStateIncrease((PetscObject)X));
933: PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
934: PetscFunctionReturn(PETSC_SUCCESS);
935: }
937: /*@
938: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
940: Not Collective
942: Input Parameters:
943: + X - nest vector
944: . idxm - index of the vector within the nest vector
945: - sx - vector at index `idxm` within the nest vector
947: Level: developer
949: Note:
950: The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
952: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
953: @*/
954: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
955: {
956: PetscFunctionBegin;
957: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
962: {
963: PetscInt i;
965: PetscFunctionBegin;
966: PetscCall(PetscObjectStateIncrease((PetscObject)X));
967: for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@C
972: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
974: Not Collective
976: Input Parameters:
977: + X - nest vector
978: . N - number of component vecs in `sx`
979: . idxm - indices of component vectors that are to be replaced
980: - sx - array of vectors
982: Level: developer
984: Note:
985: The components in the vector array `sx` do not have to be of the same size as corresponding
986: components in `X`. The user can also free the array `sx` after the call.
988: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
989: @*/
990: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
991: {
992: PetscFunctionBegin;
993: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
998: {
999: Vec_Nest *b = (Vec_Nest *)X->data;
1001: PetscFunctionBegin;
1002: *N = b->nb;
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*@
1007: VecNestGetSize - Returns the size of the nest vector.
1009: Not Collective
1011: Input Parameter:
1012: . X - nest vector
1014: Output Parameter:
1015: . N - number of nested vecs
1017: Level: developer
1019: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1020: @*/
1021: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1022: {
1023: PetscFunctionBegin;
1026: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1031: {
1032: Vec_Nest *ctx = (Vec_Nest *)V->data;
1033: PetscInt i;
1035: PetscFunctionBegin;
1036: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1038: ctx->nb = nb;
1039: PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1041: /* Create space */
1042: PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1043: for (i = 0; i < ctx->nb; i++) {
1044: ctx->v[i] = x[i];
1045: PetscCall(PetscObjectReference((PetscObject)x[i]));
1046: /* Do not allocate memory for internal sub blocks */
1047: }
1049: PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1051: ctx->setup_called = PETSC_TRUE;
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1056: {
1057: Vec_Nest *ctx = (Vec_Nest *)V->data;
1058: PetscInt i, offset, m, n, M, N;
1060: PetscFunctionBegin;
1061: if (is) { /* Do some consistency checks and reference the is */
1062: offset = V->map->rstart;
1063: for (i = 0; i < ctx->nb; i++) {
1064: PetscCall(ISGetSize(is[i], &M));
1065: PetscCall(VecGetSize(ctx->v[i], &N));
1066: PetscCheck(M == N, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of size %" PetscInt_FMT " is not compatible with Vec of size %" PetscInt_FMT, i, M, N);
1067: PetscCall(ISGetLocalSize(is[i], &m));
1068: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1069: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of local size %" PetscInt_FMT " is not compatible with Vec of local size %" PetscInt_FMT, i, m, n);
1070: PetscCall(PetscObjectReference((PetscObject)is[i]));
1071: ctx->is[i] = is[i];
1072: offset += n;
1073: }
1074: } else { /* Create a contiguous ISStride for each entry */
1075: offset = V->map->rstart;
1076: for (i = 0; i < ctx->nb; i++) {
1077: PetscInt bs;
1078: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1079: PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1080: PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1081: PetscCall(ISSetBlockSize(ctx->is[i], bs));
1082: offset += n;
1083: }
1084: }
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: /*MC
1089: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1091: Level: intermediate
1093: Notes:
1094: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1095: It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1097: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1098: M*/
1100: /*@C
1101: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1103: Collective
1105: Input Parameters:
1106: + comm - Communicator for the new `Vec`
1107: . nb - number of nested blocks
1108: . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1109: - x - array of `nb` sub-vectors
1111: Output Parameter:
1112: . Y - new vector
1114: Level: advanced
1116: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`, `VECNEST`
1117: @*/
1118: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1119: {
1120: Vec V;
1121: Vec_Nest *s;
1122: PetscInt n, N;
1124: PetscFunctionBegin;
1125: PetscCall(VecCreate(comm, &V));
1126: PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1128: /* allocate and set pointer for implementation data */
1129: PetscCall(PetscNew(&s));
1130: V->data = (void *)s;
1131: s->setup_called = PETSC_FALSE;
1132: s->nb = -1;
1133: s->v = NULL;
1135: PetscCall(VecSetUp_Nest_Private(V, nb, x));
1137: n = N = 0;
1138: PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1139: PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1140: PetscCall(PetscLayoutSetSize(V->map, N));
1141: PetscCall(PetscLayoutSetLocalSize(V->map, n));
1142: PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1143: PetscCall(PetscLayoutSetUp(V->map));
1145: PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1147: PetscCall(VecNestSetOps_Private(V->ops));
1148: V->petscnative = PETSC_FALSE;
1150: /* expose block api's */
1151: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1152: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1153: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1154: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1155: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1157: *Y = V;
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }