Actual source code: vecnest.c
1: #include <../src/vec/vec/impls/nest/vecnestimpl.h>
3: /* check all blocks are filled */
4: static PetscErrorCode VecAssemblyBegin_Nest(Vec v)
5: {
6: Vec_Nest *vs = (Vec_Nest *)v->data;
7: PetscInt i;
9: PetscFunctionBegin;
10: for (i = 0; i < vs->nb; i++) {
11: PetscCheck(vs->v[i], PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Nest vector cannot contain NULL blocks");
12: PetscCall(VecAssemblyBegin(vs->v[i]));
13: }
14: PetscFunctionReturn(PETSC_SUCCESS);
15: }
17: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
18: {
19: Vec_Nest *vs = (Vec_Nest *)v->data;
20: PetscInt i;
22: PetscFunctionBegin;
23: for (i = 0; i < vs->nb; i++) PetscCall(VecAssemblyEnd(vs->v[i]));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode VecDestroy_Nest(Vec v)
28: {
29: Vec_Nest *vs = (Vec_Nest *)v->data;
30: PetscInt i;
32: PetscFunctionBegin;
33: if (vs->v) {
34: for (i = 0; i < vs->nb; i++) PetscCall(VecDestroy(&vs->v[i]));
35: PetscCall(PetscFree(vs->v));
36: }
37: for (i = 0; i < vs->nb; i++) PetscCall(ISDestroy(&vs->is[i]));
38: PetscCall(PetscFree(vs->is));
39: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVec_C", NULL));
40: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecs_C", NULL));
41: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecsRead_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestRestoreSubVecsRead_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVec_C", NULL));
44: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL));
45: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL));
47: PetscCall(PetscFree(vs));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
52: {
53: Vec_Nest *bx = (Vec_Nest *)x->data;
54: Vec_Nest *by = (Vec_Nest *)y->data;
55: PetscInt i;
57: PetscFunctionBegin;
58: VecNestCheckCompatible2(x, 1, y, 2);
59: for (i = 0; i < bx->nb; i++) PetscCall(VecCopy(bx->v[i], by->v[i]));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
64: {
65: Vec_Nest *bx = (Vec_Nest *)x->data;
66: Vec Y;
67: Vec *sub;
68: PetscInt i;
70: PetscFunctionBegin;
71: PetscCall(PetscMalloc1(bx->nb, &sub));
72: for (i = 0; i < bx->nb; i++) PetscCall(VecDuplicate(bx->v[i], &sub[i]));
73: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y));
74: for (i = 0; i < bx->nb; i++) PetscCall(VecDestroy(&sub[i]));
75: PetscCall(PetscFree(sub));
76: *y = Y;
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
81: {
82: Vec_Nest *bx = (Vec_Nest *)x->data;
83: Vec_Nest *by = (Vec_Nest *)y->data;
84: PetscInt i, nr;
85: PetscScalar x_dot_y, _val;
87: PetscFunctionBegin;
88: VecNestCheckCompatible2(x, 1, y, 2);
89: nr = bx->nb;
90: _val = 0.0;
91: for (i = 0; i < nr; i++) {
92: PetscCall(VecDot(bx->v[i], by->v[i], &x_dot_y));
93: _val = _val + x_dot_y;
94: }
95: *val = _val;
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
100: {
101: Vec_Nest *bx = (Vec_Nest *)x->data;
102: Vec_Nest *by = (Vec_Nest *)y->data;
103: PetscInt i, nr;
104: PetscScalar x_dot_y, _val;
106: PetscFunctionBegin;
107: VecNestCheckCompatible2(x, 1, y, 2);
108: nr = bx->nb;
109: _val = 0.0;
110: for (i = 0; i < nr; i++) {
111: PetscCall(VecTDot(bx->v[i], by->v[i], &x_dot_y));
112: _val = _val + x_dot_y;
113: }
114: *val = _val;
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
119: {
120: Vec_Nest *bx = (Vec_Nest *)x->data;
121: Vec_Nest *by = (Vec_Nest *)y->data;
122: PetscInt i, nr;
123: PetscScalar x_dot_y, _dp, _nm;
124: PetscReal norm2_y;
126: PetscFunctionBegin;
127: VecNestCheckCompatible2(x, 1, y, 2);
128: nr = bx->nb;
129: _dp = 0.0;
130: _nm = 0.0;
131: for (i = 0; i < nr; i++) {
132: PetscCall(VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y));
133: _dp += x_dot_y;
134: _nm += norm2_y;
135: }
136: *dp = _dp;
137: *nm = _nm;
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode VecAXPY_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: PetscFunctionBegin;
148: VecNestCheckCompatible2(y, 1, x, 3);
149: nr = bx->nb;
150: for (i = 0; i < nr; i++) PetscCall(VecAXPY(by->v[i], alpha, bx->v[i]));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
155: {
156: Vec_Nest *bx = (Vec_Nest *)x->data;
157: Vec_Nest *by = (Vec_Nest *)y->data;
158: PetscInt i, nr;
160: PetscFunctionBegin;
161: VecNestCheckCompatible2(y, 1, x, 3);
162: nr = bx->nb;
163: for (i = 0; i < nr; i++) PetscCall(VecAYPX(by->v[i], alpha, bx->v[i]));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
168: {
169: Vec_Nest *bx = (Vec_Nest *)x->data;
170: Vec_Nest *by = (Vec_Nest *)y->data;
171: PetscInt i, nr;
173: PetscFunctionBegin;
174: VecNestCheckCompatible2(y, 1, x, 4);
175: nr = bx->nb;
176: for (i = 0; i < nr; i++) PetscCall(VecAXPBY(by->v[i], alpha, beta, bx->v[i]));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
181: {
182: Vec_Nest *bx = (Vec_Nest *)x->data;
183: Vec_Nest *by = (Vec_Nest *)y->data;
184: Vec_Nest *bz = (Vec_Nest *)z->data;
185: PetscInt i, nr;
187: PetscFunctionBegin;
188: VecNestCheckCompatible3(z, 1, x, 5, y, 6);
189: nr = bx->nb;
190: for (i = 0; i < nr; i++) PetscCall(VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
195: {
196: Vec_Nest *bx = (Vec_Nest *)x->data;
197: PetscInt i, nr;
199: PetscFunctionBegin;
200: nr = bx->nb;
201: for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
206: {
207: Vec_Nest *bx = (Vec_Nest *)x->data;
208: Vec_Nest *by = (Vec_Nest *)y->data;
209: Vec_Nest *bw = (Vec_Nest *)w->data;
210: PetscInt i, nr;
212: PetscFunctionBegin;
213: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
214: nr = bx->nb;
215: for (i = 0; i < nr; i++) PetscCall(VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
220: {
221: Vec_Nest *bx = (Vec_Nest *)x->data;
222: Vec_Nest *by = (Vec_Nest *)y->data;
223: Vec_Nest *bw = (Vec_Nest *)w->data;
224: PetscInt i, nr;
226: PetscFunctionBegin;
227: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
228: nr = bx->nb;
229: for (i = 0; i < nr; i++) PetscCall(VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]));
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: static PetscErrorCode VecReciprocal_Nest(Vec x)
234: {
235: Vec_Nest *bx = (Vec_Nest *)x->data;
236: PetscInt i, nr;
238: PetscFunctionBegin;
239: nr = bx->nb;
240: for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
245: {
246: Vec_Nest *bx = (Vec_Nest *)xin->data;
247: PetscInt i, nr;
248: PetscReal z_i;
249: PetscReal _z;
251: PetscFunctionBegin;
252: nr = bx->nb;
253: _z = 0.0;
255: if (type == NORM_2) {
256: PetscScalar dot;
257: PetscCall(VecDot(xin, xin, &dot));
258: _z = PetscAbsScalar(PetscSqrtScalar(dot));
259: } else if (type == NORM_1) {
260: for (i = 0; i < nr; i++) {
261: PetscCall(VecNorm(bx->v[i], type, &z_i));
262: _z = _z + z_i;
263: }
264: } else if (type == NORM_INFINITY) {
265: for (i = 0; i < nr; i++) {
266: PetscCall(VecNorm(bx->v[i], type, &z_i));
267: if (z_i > _z) _z = z_i;
268: }
269: }
271: *z = _z;
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
276: {
277: PetscFunctionBegin;
278: /* TODO: implement proper MAXPY. For now, do axpy on each vector */
279: for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
284: {
285: PetscFunctionBegin;
286: /* TODO: implement proper MDOT. For now, do dot on each vector */
287: for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
292: {
293: PetscFunctionBegin;
294: /* TODO: implement proper MTDOT. For now, do tdot on each vector */
295: for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
300: {
301: Vec_Nest *bx = (Vec_Nest *)x->data;
302: PetscInt i, nr;
304: PetscFunctionBegin;
305: nr = bx->nb;
306: for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode VecConjugate_Nest(Vec x)
311: {
312: Vec_Nest *bx = (Vec_Nest *)x->data;
313: PetscInt j, nr;
315: PetscFunctionBegin;
316: nr = bx->nb;
317: for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
322: {
323: Vec_Nest *bx = (Vec_Nest *)x->data;
324: Vec_Nest *by = (Vec_Nest *)y->data;
325: PetscInt i, nr;
327: PetscFunctionBegin;
328: VecNestCheckCompatible2(x, 1, y, 2);
329: nr = bx->nb;
330: for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
335: {
336: Vec_Nest *bx = (Vec_Nest *)x->data;
337: Vec_Nest *by = (Vec_Nest *)y->data;
338: Vec_Nest *bw = (Vec_Nest *)w->data;
339: PetscInt i, nr;
341: PetscFunctionBegin;
342: VecNestCheckCompatible3(w, 1, x, 3, y, 4);
343: nr = bx->nb;
344: for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
349: {
350: PetscInt i, lp = -1, lw = -1;
351: PetscReal lv;
352: Vec_Nest *bx = (Vec_Nest *)x->data;
354: PetscFunctionBegin;
355: *v = PETSC_MAX_REAL;
356: for (i = 0; i < bx->nb; i++) {
357: PetscInt tp;
358: PetscCall(VecMin(bx->v[i], &tp, &lv));
359: if (lv < *v) {
360: *v = lv;
361: lw = i;
362: lp = tp;
363: }
364: }
365: if (p && lw > -1) {
366: PetscInt st, en;
367: const PetscInt *idxs;
369: *p = -1;
370: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
371: if (lp >= st && lp < en) {
372: PetscCall(ISGetIndices(bx->is[lw], &idxs));
373: *p = idxs[lp - st];
374: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
375: }
376: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
377: }
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
382: {
383: PetscInt i, lp = -1, lw = -1;
384: PetscReal lv;
385: Vec_Nest *bx = (Vec_Nest *)x->data;
387: PetscFunctionBegin;
388: *v = PETSC_MIN_REAL;
389: for (i = 0; i < bx->nb; i++) {
390: PetscInt tp;
391: PetscCall(VecMax(bx->v[i], &tp, &lv));
392: if (lv > *v) {
393: *v = lv;
394: lw = i;
395: lp = tp;
396: }
397: }
398: if (p && lw > -1) {
399: PetscInt st, en;
400: const PetscInt *idxs;
402: *p = -1;
403: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
404: if (lp >= st && lp < en) {
405: PetscCall(ISGetIndices(bx->is[lw], &idxs));
406: *p = idxs[lp - st];
407: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
408: }
409: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
410: }
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
415: {
416: Vec_Nest *bx = (Vec_Nest *)x->data;
417: PetscBool isascii;
418: PetscInt i;
420: PetscFunctionBegin;
421: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
422: if (isascii) {
423: PetscCall(PetscViewerASCIIPushTab(viewer));
424: PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ", structure:\n", bx->nb));
425: for (i = 0; i < bx->nb; i++) {
426: VecType type;
427: char name[256] = "", prefix[256] = "";
428: PetscInt NR;
430: PetscCall(VecGetSize(bx->v[i], &NR));
431: PetscCall(VecGetType(bx->v[i], &type));
432: if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
433: if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));
435: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT "\n", i, name, prefix, type, NR));
437: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
438: PetscCall(VecView(bx->v[i], viewer));
439: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
440: }
441: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
442: }
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
447: {
448: Vec_Nest *bx;
449: PetscInt size, i, nr;
450: PetscBool isnest;
452: PetscFunctionBegin;
453: PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
454: if (!isnest) {
455: /* Not nest */
456: if (globalsize) PetscCall(VecGetSize(x, &size));
457: else PetscCall(VecGetLocalSize(x, &size));
458: *L = *L + size;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /* Otherwise we have a nest */
463: bx = (Vec_Nest *)x->data;
464: nr = bx->nb;
466: /* now descend recursively */
467: for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /* Returns the sum of the global size of all the constituent vectors in the nest */
472: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
473: {
474: PetscFunctionBegin;
475: *N = x->map->N;
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: /* Returns the sum of the local size of all the constituent vectors in the nest */
480: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
481: {
482: PetscFunctionBegin;
483: *n = x->map->n;
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
488: {
489: Vec_Nest *bx = (Vec_Nest *)x->data;
490: Vec_Nest *by = (Vec_Nest *)y->data;
491: PetscInt i, nr;
492: PetscReal local_max, m;
494: PetscFunctionBegin;
495: VecNestCheckCompatible2(x, 1, y, 2);
496: nr = bx->nb;
497: m = 0.0;
498: for (i = 0; i < nr; i++) {
499: PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
500: if (local_max > m) m = local_max;
501: }
502: *max = m;
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
507: {
508: Vec_Nest *bx = (Vec_Nest *)X->data;
509: PetscInt i;
511: PetscFunctionBegin;
512: *x = NULL;
513: for (i = 0; i < bx->nb; i++) {
514: PetscBool issame = PETSC_FALSE;
515: PetscCall(ISEqual(is, bx->is[i], &issame));
516: if (issame) {
517: *x = bx->v[i];
518: PetscCall(PetscObjectReference((PetscObject)*x));
519: break;
520: }
521: }
522: PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
527: {
528: PetscFunctionBegin;
529: PetscCall(PetscObjectStateIncrease((PetscObject)X));
530: PetscCall(VecDestroy(x));
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
535: {
536: Vec_Nest *bx = (Vec_Nest *)X->data;
537: PetscInt i, m, rstart, rend;
539: PetscFunctionBegin;
540: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
541: PetscCall(VecGetLocalSize(X, &m));
542: PetscCall(PetscMalloc1(m, x));
543: for (i = 0; i < bx->nb; i++) {
544: Vec subvec = bx->v[i];
545: IS isy = bx->is[i];
546: PetscInt j, sm;
547: const PetscInt *ixy;
548: const PetscScalar *y;
549: PetscCall(VecGetLocalSize(subvec, &sm));
550: PetscCall(VecGetArrayRead(subvec, &y));
551: PetscCall(ISGetIndices(isy, &ixy));
552: for (j = 0; j < sm; j++) {
553: PetscInt ix = ixy[j];
554: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
555: (*x)[ix - rstart] = y[j];
556: }
557: PetscCall(ISRestoreIndices(isy, &ixy));
558: PetscCall(VecRestoreArrayRead(subvec, &y));
559: }
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
564: {
565: Vec_Nest *bx = (Vec_Nest *)X->data;
566: PetscInt i, m, rstart, rend;
568: PetscFunctionBegin;
569: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
570: PetscCall(VecGetLocalSize(X, &m));
571: for (i = 0; i < bx->nb; i++) {
572: Vec subvec = bx->v[i];
573: IS isy = bx->is[i];
574: PetscInt j, sm;
575: const PetscInt *ixy;
576: PetscScalar *y;
577: PetscCall(VecGetLocalSize(subvec, &sm));
578: PetscCall(VecGetArray(subvec, &y));
579: PetscCall(ISGetIndices(isy, &ixy));
580: for (j = 0; j < sm; j++) {
581: PetscInt ix = ixy[j];
582: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
583: y[j] = (*x)[ix - rstart];
584: }
585: PetscCall(ISRestoreIndices(isy, &ixy));
586: PetscCall(VecRestoreArray(subvec, &y));
587: }
588: PetscCall(PetscFree(*x));
589: PetscCall(PetscObjectStateIncrease((PetscObject)X));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
594: {
595: PetscFunctionBegin;
596: PetscCall(PetscFree(*x));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
601: {
602: PetscFunctionBegin;
603: PetscCheck(nx <= 0, PetscObjectComm((PetscObject)*X), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
608: {
609: Vec *ww;
610: IS *wis;
611: Vec_Nest *bv = (Vec_Nest *)v->data;
612: PetscInt i;
613: PetscMPIInt size;
615: PetscFunctionBegin;
616: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
617: if (size == 1) {
618: PetscCall(VecDuplicate(v, w));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
621: PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
622: for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
623: for (i = 0; i < bv->nb; i++) {
624: PetscInt off;
626: PetscCall(VecGetOwnershipRange(v, &off, NULL));
627: PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
628: PetscCall(ISShift(wis[i], -off, wis[i]));
629: }
630: PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
631: for (i = 0; i < bv->nb; i++) {
632: PetscCall(VecDestroy(&ww[i]));
633: PetscCall(ISDestroy(&wis[i]));
634: }
635: PetscCall(PetscFree2(ww, wis));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
640: {
641: Vec_Nest *bv = (Vec_Nest *)v->data;
642: Vec_Nest *bw = (Vec_Nest *)w->data;
643: PetscInt i;
645: PetscFunctionBegin;
646: PetscCheckSameType(v, 1, w, 2);
647: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
648: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
653: {
654: Vec_Nest *bv = (Vec_Nest *)v->data;
655: Vec_Nest *bw = (Vec_Nest *)w->data;
656: PetscInt i;
658: PetscFunctionBegin;
659: PetscCheckSameType(v, 1, w, 2);
660: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
661: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
662: PetscCall(PetscObjectStateIncrease((PetscObject)v));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
667: {
668: Vec_Nest *bv = (Vec_Nest *)v->data;
669: Vec_Nest *bw = (Vec_Nest *)w->data;
670: PetscInt i;
672: PetscFunctionBegin;
673: PetscCheckSameType(v, 1, w, 2);
674: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
675: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }
679: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
680: {
681: Vec_Nest *bv = (Vec_Nest *)v->data;
682: Vec_Nest *bw = (Vec_Nest *)w->data;
683: PetscInt i;
685: PetscFunctionBegin;
686: PetscCheckSameType(v, 1, w, 2);
687: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
688: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
693: {
694: Vec_Nest *bv = (Vec_Nest *)v->data;
696: PetscFunctionBegin;
697: for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode VecErrorWeightedNorms_Nest(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
702: {
703: Vec_Nest *bu = (Vec_Nest *)U->data;
704: Vec_Nest *by = (Vec_Nest *)Y->data;
705: Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
706: Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
707: Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;
709: PetscFunctionBegin;
710: VecNestCheckCompatible2(U, 1, Y, 2);
711: if (E) VecNestCheckCompatible2(U, 1, E, 3);
712: if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
713: if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
714: *norm = 0.0;
715: *norma = 0.0;
716: *normr = 0.0;
717: *norm_loc = 0;
718: *norma_loc = 0;
719: *normr_loc = 0;
720: for (PetscInt i = 0; i < bu->nb; i++) {
721: PetscReal n, na, nr;
722: PetscInt n_loc, na_loc, nr_loc;
724: PetscCall(VecErrorWeightedNorms(bu->v[i], by->v[i], be ? be->v[i] : NULL, wnormtype, atol, ba ? ba->v[i] : NULL, rtol, br ? br->v[i] : NULL, ignore_max, &n, &n_loc, &na, &na_loc, &nr, &nr_loc));
725: if (wnormtype == NORM_INFINITY) {
726: *norm = PetscMax(*norm, n);
727: *norma = PetscMax(*norma, na);
728: *normr = PetscMax(*normr, nr);
729: } else {
730: *norm += PetscSqr(n);
731: *norma += PetscSqr(na);
732: *normr += PetscSqr(nr);
733: }
734: *norm_loc += n_loc;
735: *norma_loc += na_loc;
736: *normr_loc += nr_loc;
737: }
738: if (wnormtype == NORM_2) {
739: *norm = PetscSqrtReal(*norm);
740: *norma = PetscSqrtReal(*norma);
741: *normr = PetscSqrtReal(*normr);
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
747: {
748: PetscFunctionBegin;
749: ops->duplicate = VecDuplicate_Nest;
750: ops->duplicatevecs = VecDuplicateVecs_Default;
751: ops->destroyvecs = VecDestroyVecs_Default;
752: ops->dot = VecDot_Nest;
753: ops->mdot = VecMDot_Nest;
754: ops->norm = VecNorm_Nest;
755: ops->tdot = VecTDot_Nest;
756: ops->mtdot = VecMTDot_Nest;
757: ops->scale = VecScale_Nest;
758: ops->copy = VecCopy_Nest;
759: ops->set = VecSet_Nest;
760: ops->swap = VecSwap_Nest;
761: ops->axpy = VecAXPY_Nest;
762: ops->axpby = VecAXPBY_Nest;
763: ops->maxpy = VecMAXPY_Nest;
764: ops->aypx = VecAYPX_Nest;
765: ops->waxpy = VecWAXPY_Nest;
766: ops->axpbypcz = NULL;
767: ops->pointwisemult = VecPointwiseMult_Nest;
768: ops->pointwisedivide = VecPointwiseDivide_Nest;
769: ops->setvalues = NULL;
770: ops->assemblybegin = VecAssemblyBegin_Nest;
771: ops->assemblyend = VecAssemblyEnd_Nest;
772: ops->getarray = VecGetArray_Nest;
773: ops->getsize = VecGetSize_Nest;
774: ops->getlocalsize = VecGetLocalSize_Nest;
775: ops->restorearray = VecRestoreArray_Nest;
776: ops->restorearrayread = VecRestoreArrayRead_Nest;
777: ops->max = VecMax_Nest;
778: ops->min = VecMin_Nest;
779: ops->setrandom = NULL;
780: ops->setoption = NULL;
781: ops->setvaluesblocked = NULL;
782: ops->destroy = VecDestroy_Nest;
783: ops->view = VecView_Nest;
784: ops->placearray = NULL;
785: ops->replacearray = NULL;
786: ops->dot_local = VecDot_Nest;
787: ops->tdot_local = VecTDot_Nest;
788: ops->norm_local = VecNorm_Nest;
789: ops->mdot_local = VecMDot_Nest;
790: ops->mtdot_local = VecMTDot_Nest;
791: ops->load = NULL;
792: ops->reciprocal = VecReciprocal_Nest;
793: ops->conjugate = VecConjugate_Nest;
794: ops->setlocaltoglobalmapping = NULL;
795: ops->resetarray = NULL;
796: ops->setfromoptions = NULL;
797: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
798: ops->load = NULL;
799: ops->pointwisemax = NULL;
800: ops->pointwisemaxabs = NULL;
801: ops->pointwisemin = NULL;
802: ops->getvalues = NULL;
803: ops->sqrt = NULL;
804: ops->abs = NULL;
805: ops->exp = NULL;
806: ops->shift = NULL;
807: ops->create = NULL;
808: ops->stridegather = NULL;
809: ops->stridescatter = NULL;
810: ops->dotnorm2 = VecDotNorm2_Nest;
811: ops->getsubvector = VecGetSubVector_Nest;
812: ops->restoresubvector = VecRestoreSubVector_Nest;
813: ops->axpbypcz = VecAXPBYPCZ_Nest;
814: ops->concatenate = VecConcatenate_Nest;
815: ops->createlocalvector = VecCreateLocalVector_Nest;
816: ops->getlocalvector = VecGetLocalVector_Nest;
817: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
818: ops->restorelocalvector = VecRestoreLocalVector_Nest;
819: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
820: ops->setrandom = VecSetRandom_Nest;
821: ops->errorwnorm = VecErrorWeightedNorms_Nest;
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
826: {
827: Vec_Nest *b = (Vec_Nest *)x->data;
828: PetscInt i;
829: PetscInt row;
831: PetscFunctionBegin;
832: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
833: for (i = 0; i < m; i++) {
834: row = idxm[i];
835: PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
836: vec[i] = b->v[row];
837: }
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
842: {
843: PetscFunctionBegin;
844: PetscCall(PetscObjectStateIncrease((PetscObject)X));
845: PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: /*@
850: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
852: Not Collective
854: Input Parameters:
855: + X - nest vector
856: - idxm - index of the vector within the nest
858: Output Parameter:
859: . sx - vector at index `idxm` within the nest
861: Level: developer
863: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
864: @*/
865: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
866: {
867: PetscFunctionBegin;
868: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
873: {
874: Vec_Nest *b = (Vec_Nest *)X->data;
876: PetscFunctionBegin;
877: PetscCall(PetscObjectStateIncrease((PetscObject)X));
878: if (N) *N = b->nb;
879: if (sx) *sx = b->v;
880: PetscFunctionReturn(PETSC_SUCCESS);
881: }
883: /*@C
884: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
886: Not Collective
888: Input Parameter:
889: . X - nest vector
891: Output Parameters:
892: + N - number of nested vecs
893: - sx - array of vectors, can pass in `NULL`
895: Level: developer
897: Note:
898: The user should not free the array `sx`.
900: Fortran Notes:
901: The caller must allocate the array to hold the subvectors and pass it in.
903: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestGetSubVecsRead()`
904: @*/
905: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
906: {
907: PetscFunctionBegin;
908: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
913: {
914: Vec_Nest *bx = (Vec_Nest *)X->data;
915: PetscInt i, offset = 0, n = 0, bs;
916: IS is;
917: PetscBool issame = PETSC_FALSE;
918: PetscInt N = 0;
920: PetscFunctionBegin;
921: /* check if idxm < bx->nb */
922: PetscCheck(idxm < bx->nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT ", should be less than %" PetscInt_FMT, idxm, bx->nb);
924: PetscCall(PetscObjectReference((PetscObject)x));
925: PetscCall(VecDestroy(&bx->v[idxm]));
926: bx->v[idxm] = x;
928: /* check if we need to update the IS for the block */
929: offset = X->map->rstart;
930: for (i = 0; i < idxm; i++) {
931: n = 0;
932: PetscCall(VecGetLocalSize(bx->v[i], &n));
933: offset += n;
934: }
936: /* get the local size and block size */
937: PetscCall(VecGetLocalSize(x, &n));
938: PetscCall(VecGetBlockSize(x, &bs));
940: /* create the new IS */
941: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
942: PetscCall(ISSetBlockSize(is, bs));
944: /* check if they are equal */
945: PetscCall(ISEqual(is, bx->is[idxm], &issame));
947: if (!issame) {
948: /* The IS of given vector has a different layout compared to the existing block vector.
949: Destroy the existing reference and update the IS. */
950: PetscCall(ISDestroy(&bx->is[idxm]));
951: PetscCall(ISDuplicate(is, &bx->is[idxm]));
952: PetscCall(ISCopy(is, bx->is[idxm]));
954: offset += n;
955: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
956: for (i = idxm + 1; i < bx->nb; i++) {
957: /* get the local size and block size */
958: PetscCall(VecGetLocalSize(bx->v[i], &n));
959: PetscCall(VecGetBlockSize(bx->v[i], &bs));
961: /* destroy the old and create the new IS */
962: PetscCall(ISDestroy(&bx->is[i]));
963: PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
964: PetscCall(ISSetBlockSize(bx->is[i], bs));
966: offset += n;
967: }
969: n = 0;
970: PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
971: PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
972: PetscCall(PetscLayoutSetSize(X->map, N));
973: PetscCall(PetscLayoutSetLocalSize(X->map, n));
974: }
976: PetscCall(ISDestroy(&is));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
981: {
982: PetscFunctionBegin;
983: PetscCall(PetscObjectStateIncrease((PetscObject)X));
984: PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /*@
989: VecNestGetSubVecsRead - Access the subvecs of a `VECNEST` vector for read-only access
991: Logically collective
993: Input Parameter:
994: . X - nest vector
996: Output Parameters:
997: + N - number of nested vecs
998: - sx - array of read-locked vectors
1000: Level: advanced
1002: Notes:
1003: Each of the subvecs will be read locked (`VecLockReadPush()`), which is a logically collective operation.
1004: When access is complete, you must call `VecNestRestoreSubVecsRead()` to release the locks.
1006: Developer Note:
1007: This function does not increase the state of `X` (`PetscObjectStateIncrease()`).
1009: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestRestoreSubVecsRead()`
1010: @*/
1011: PetscErrorCode VecNestGetSubVecsRead(Vec X, PetscInt *N, Vec *sx[])
1012: {
1013: PetscFunctionBegin;
1014: PetscUseMethod(X, "VecNestGetSubVecsRead_C", (Vec, PetscInt *, Vec **), (X, N, sx));
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: static PetscErrorCode VecNestGetSubVecsRead_Nest(Vec X, PetscInt *N, Vec **sx)
1019: {
1020: Vec_Nest *b = (Vec_Nest *)X->data;
1022: PetscFunctionBegin;
1023: for (PetscInt i = 0; i < b->nb; i++) PetscCall(VecLockReadPush(b->v[i]));
1024: if (N) *N = b->nb;
1025: if (sx) *sx = b->v;
1026: PetscFunctionReturn(PETSC_SUCCESS);
1027: }
1029: /*@
1030: VecNestRestoreSubVecsRead - Restore access the subvecs of a `VECNEST` vector obtained with `VecNestGetSubVecsRead()`
1032: Logically collective
1034: Input Parameters:
1035: + X - nest vector
1036: . N - number of nested vecs
1037: - sx - array of read-locked vectors
1039: Level: advanced
1041: Note:
1042: The same arguments to `VecNestGetSubVecsRead()` should be the argument to `VecNestRestoreSubVecsRead()`.
1044: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestGetSubVecsRead()`
1045: @*/
1046: PetscErrorCode VecNestRestoreSubVecsRead(Vec X, PetscInt *N, Vec *sx[])
1047: {
1048: PetscFunctionBegin;
1049: PetscUseMethod(X, "VecNestRestoreSubVecsRead_C", (Vec, PetscInt *, Vec **), (X, N, sx));
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: static PetscErrorCode VecNestRestoreSubVecsRead_Nest(Vec X, PetscInt *N, Vec **sx)
1054: {
1055: Vec_Nest *b = (Vec_Nest *)X->data;
1057: PetscFunctionBegin;
1058: if (N) *N = 0;
1059: if (sx) {
1060: PetscCheck(*sx == b->v, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONG, "Restoring incorrect array of vectors");
1061: *sx = NULL;
1062: }
1063: for (PetscInt i = 0; i < b->nb; i++) PetscCall(VecLockReadPop(b->v[i]));
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: /*@
1068: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
1070: Not Collective
1072: Input Parameters:
1073: + X - nest vector
1074: . idxm - index of the vector within the nest vector
1075: - sx - vector at index `idxm` within the nest vector
1077: Level: developer
1079: Notes:
1080: The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1082: The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.
1084: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
1085: @*/
1086: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
1087: {
1088: PetscFunctionBegin;
1089: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1090: PetscFunctionReturn(PETSC_SUCCESS);
1091: }
1093: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1094: {
1095: PetscInt i;
1097: PetscFunctionBegin;
1098: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1099: for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: /*@
1104: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1106: Not Collective
1108: Input Parameters:
1109: + X - nest vector
1110: . N - number of component vecs in `sx`
1111: . idxm - indices of component vectors that are to be replaced
1112: - sx - array of vectors
1114: Level: developer
1116: Notes:
1117: The components in the vector array `sx` do not have to be of the same size as corresponding
1118: components in `X`. The user can also free the array `sx` after the call.
1120: The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.
1122: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1123: @*/
1124: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
1125: {
1126: PetscFunctionBegin;
1127: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1132: {
1133: Vec_Nest *b = (Vec_Nest *)X->data;
1135: PetscFunctionBegin;
1136: *N = b->nb;
1137: PetscFunctionReturn(PETSC_SUCCESS);
1138: }
1140: /*@
1141: VecNestGetSize - Returns the size of the nest vector.
1143: Not Collective
1145: Input Parameter:
1146: . X - nest vector
1148: Output Parameter:
1149: . N - number of nested vecs
1151: Level: developer
1153: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1154: @*/
1155: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1156: {
1157: PetscFunctionBegin;
1159: PetscAssertPointer(N, 2);
1160: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1165: {
1166: Vec_Nest *ctx = (Vec_Nest *)V->data;
1167: PetscInt i;
1169: PetscFunctionBegin;
1170: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1172: ctx->nb = nb;
1173: PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1175: /* Create space */
1176: PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1177: for (i = 0; i < ctx->nb; i++) {
1178: ctx->v[i] = x[i];
1179: PetscCall(PetscObjectReference((PetscObject)x[i]));
1180: /* Do not allocate memory for internal sub blocks */
1181: }
1183: PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1185: ctx->setup_called = PETSC_TRUE;
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1190: {
1191: Vec_Nest *ctx = (Vec_Nest *)V->data;
1192: PetscInt i, offset, m, n, M, N;
1194: PetscFunctionBegin;
1195: (void)nb;
1196: if (is) { /* Do some consistency checks and reference the is */
1197: offset = V->map->rstart;
1198: for (i = 0; i < ctx->nb; i++) {
1199: PetscCall(ISGetSize(is[i], &M));
1200: PetscCall(VecGetSize(ctx->v[i], &N));
1201: 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);
1202: PetscCall(ISGetLocalSize(is[i], &m));
1203: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1204: 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);
1205: PetscCall(PetscObjectReference((PetscObject)is[i]));
1206: ctx->is[i] = is[i];
1207: offset += n;
1208: }
1209: } else { /* Create a contiguous ISStride for each entry */
1210: offset = V->map->rstart;
1211: for (i = 0; i < ctx->nb; i++) {
1212: PetscInt bs;
1213: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1214: PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1215: PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1216: PetscCall(ISSetBlockSize(ctx->is[i], bs));
1217: offset += n;
1218: }
1219: }
1220: PetscFunctionReturn(PETSC_SUCCESS);
1221: }
1223: /*MC
1224: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1226: Level: intermediate
1228: Notes:
1229: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1230: It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1232: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecCreateNest()`, `MatCreateNest()`
1233: M*/
1235: /*@C
1236: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1238: Collective
1240: Input Parameters:
1241: + comm - Communicator for the new `Vec`
1242: . nb - number of nested blocks
1243: . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1244: - x - array of `nb` sub-vectors
1246: Output Parameter:
1247: . Y - new vector
1249: Level: advanced
1251: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1252: @*/
1253: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1254: {
1255: Vec V;
1256: Vec_Nest *s;
1257: PetscInt n, N;
1259: PetscFunctionBegin;
1260: PetscCall(VecCreate(comm, &V));
1261: PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1263: /* allocate and set pointer for implementation data */
1264: PetscCall(PetscNew(&s));
1265: V->data = (void *)s;
1266: s->setup_called = PETSC_FALSE;
1267: s->nb = -1;
1268: s->v = NULL;
1270: PetscCall(VecSetUp_Nest_Private(V, nb, x));
1272: n = N = 0;
1273: PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1274: PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1275: PetscCall(PetscLayoutSetSize(V->map, N));
1276: PetscCall(PetscLayoutSetLocalSize(V->map, n));
1277: PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1278: PetscCall(PetscLayoutSetUp(V->map));
1280: PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1282: PetscCall(VecNestSetOps_Private(V->ops));
1283: V->petscnative = PETSC_FALSE;
1285: /* expose block api's */
1286: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1287: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1288: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecsRead_C", VecNestGetSubVecsRead_Nest));
1289: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestRestoreSubVecsRead_C", VecNestRestoreSubVecsRead_Nest));
1290: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1291: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1292: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1294: *Y = V;
1295: PetscFunctionReturn(PETSC_SUCCESS);
1296: }