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 || type == NORM_FROBENIUS) {
256: PetscScalar dot;
258: PetscCall(VecDot(xin, xin, &dot));
259: _z = PetscAbsScalar(PetscSqrtScalar(dot));
260: } else if (type == NORM_1) {
261: for (i = 0; i < nr; i++) {
262: PetscCall(VecNorm(bx->v[i], type, &z_i));
263: _z = _z + z_i;
264: }
265: } else if (type == NORM_INFINITY) {
266: for (i = 0; i < nr; i++) {
267: PetscCall(VecNorm(bx->v[i], type, &z_i));
268: if (z_i > _z) _z = z_i;
269: }
270: } else if (type == NORM_1_AND_2) {
271: PetscCall(VecNorm_Nest(xin, NORM_1, &_z));
272: PetscCall(VecNorm_Nest(xin, NORM_2, z + 1));
273: }
275: *z = _z;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
280: {
281: PetscFunctionBegin;
282: /* TODO: implement proper MAXPY. For now, do axpy on each vector */
283: for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
288: {
289: PetscFunctionBegin;
290: /* TODO: implement proper MDOT. For now, do dot on each vector */
291: for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
296: {
297: PetscFunctionBegin;
298: /* TODO: implement proper MTDOT. For now, do tdot on each vector */
299: for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
304: {
305: Vec_Nest *bx = (Vec_Nest *)x->data;
306: PetscInt i, nr;
308: PetscFunctionBegin;
309: nr = bx->nb;
310: for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: static PetscErrorCode VecConjugate_Nest(Vec x)
315: {
316: Vec_Nest *bx = (Vec_Nest *)x->data;
317: PetscInt j, nr;
319: PetscFunctionBegin;
320: nr = bx->nb;
321: for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
326: {
327: Vec_Nest *bx = (Vec_Nest *)x->data;
328: Vec_Nest *by = (Vec_Nest *)y->data;
329: PetscInt i, nr;
331: PetscFunctionBegin;
332: VecNestCheckCompatible2(x, 1, y, 2);
333: nr = bx->nb;
334: for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
339: {
340: Vec_Nest *bx = (Vec_Nest *)x->data;
341: Vec_Nest *by = (Vec_Nest *)y->data;
342: Vec_Nest *bw = (Vec_Nest *)w->data;
343: PetscInt i, nr;
345: PetscFunctionBegin;
346: VecNestCheckCompatible3(w, 1, x, 3, y, 4);
347: nr = bx->nb;
348: for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
353: {
354: PetscInt i, lp = -1, lw = -1;
355: PetscReal lv;
356: Vec_Nest *bx = (Vec_Nest *)x->data;
358: PetscFunctionBegin;
359: *v = PETSC_MAX_REAL;
360: for (i = 0; i < bx->nb; i++) {
361: PetscInt tp;
362: PetscCall(VecMin(bx->v[i], &tp, &lv));
363: if (lv < *v) {
364: *v = lv;
365: lw = i;
366: lp = tp;
367: }
368: }
369: if (p && lw > -1) {
370: PetscInt st, en;
371: const PetscInt *idxs;
373: *p = -1;
374: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
375: if (lp >= st && lp < en) {
376: PetscCall(ISGetIndices(bx->is[lw], &idxs));
377: *p = idxs[lp - st];
378: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
379: }
380: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
386: {
387: PetscInt i, lp = -1, lw = -1;
388: PetscReal lv;
389: Vec_Nest *bx = (Vec_Nest *)x->data;
391: PetscFunctionBegin;
392: *v = PETSC_MIN_REAL;
393: for (i = 0; i < bx->nb; i++) {
394: PetscInt tp;
395: PetscCall(VecMax(bx->v[i], &tp, &lv));
396: if (lv > *v) {
397: *v = lv;
398: lw = i;
399: lp = tp;
400: }
401: }
402: if (p && lw > -1) {
403: PetscInt st, en;
404: const PetscInt *idxs;
406: *p = -1;
407: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
408: if (lp >= st && lp < en) {
409: PetscCall(ISGetIndices(bx->is[lw], &idxs));
410: *p = idxs[lp - st];
411: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
412: }
413: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
414: }
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
419: {
420: Vec_Nest *bx = (Vec_Nest *)x->data;
421: PetscBool isascii;
422: PetscInt i;
424: PetscFunctionBegin;
425: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
426: if (isascii) {
427: PetscCall(PetscViewerASCIIPushTab(viewer));
428: PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ", structure:\n", bx->nb));
429: for (i = 0; i < bx->nb; i++) {
430: VecType type;
431: char name[256] = "", prefix[256] = "";
432: PetscInt NR;
434: PetscCall(VecGetSize(bx->v[i], &NR));
435: PetscCall(VecGetType(bx->v[i], &type));
436: if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
437: if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));
439: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT "\n", i, name, prefix, type, NR));
441: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
442: PetscCall(VecView(bx->v[i], viewer));
443: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
444: }
445: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
451: {
452: Vec_Nest *bx;
453: PetscInt size, i, nr;
454: PetscBool isnest;
456: PetscFunctionBegin;
457: PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
458: if (!isnest) {
459: /* Not nest */
460: if (globalsize) PetscCall(VecGetSize(x, &size));
461: else PetscCall(VecGetLocalSize(x, &size));
462: *L = *L + size;
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /* Otherwise we have a nest */
467: bx = (Vec_Nest *)x->data;
468: nr = bx->nb;
470: /* now descend recursively */
471: for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /* Returns the sum of the global size of all the constituent vectors in the nest */
476: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
477: {
478: PetscFunctionBegin;
479: *N = x->map->N;
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: /* Returns the sum of the local size of all the constituent vectors in the nest */
484: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
485: {
486: PetscFunctionBegin;
487: *n = x->map->n;
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
492: {
493: Vec_Nest *bx = (Vec_Nest *)x->data;
494: Vec_Nest *by = (Vec_Nest *)y->data;
495: PetscInt i, nr;
496: PetscReal local_max, m;
498: PetscFunctionBegin;
499: VecNestCheckCompatible2(x, 1, y, 2);
500: nr = bx->nb;
501: m = 0.0;
502: for (i = 0; i < nr; i++) {
503: PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
504: if (local_max > m) m = local_max;
505: }
506: *max = m;
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
511: {
512: Vec_Nest *bx = (Vec_Nest *)X->data;
513: PetscInt i;
515: PetscFunctionBegin;
516: *x = NULL;
517: for (i = 0; i < bx->nb; i++) {
518: PetscBool issame = PETSC_FALSE;
519: PetscCall(ISEqual(is, bx->is[i], &issame));
520: if (issame) {
521: *x = bx->v[i];
522: PetscCall(PetscObjectReference((PetscObject)*x));
523: break;
524: }
525: }
526: PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
531: {
532: PetscFunctionBegin;
533: PetscCall(PetscObjectStateIncrease((PetscObject)X));
534: PetscCall(VecDestroy(x));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
539: {
540: Vec_Nest *bx = (Vec_Nest *)X->data;
541: PetscInt i, m, rstart, rend;
543: PetscFunctionBegin;
544: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
545: PetscCall(VecGetLocalSize(X, &m));
546: PetscCall(PetscMalloc1(m, x));
547: for (i = 0; i < bx->nb; i++) {
548: Vec subvec = bx->v[i];
549: IS isy = bx->is[i];
550: PetscInt sm;
551: const PetscInt *ixy;
552: const PetscScalar *y;
553: PetscCall(VecGetLocalSize(subvec, &sm));
554: PetscCall(VecGetArrayRead(subvec, &y));
555: PetscCall(ISGetIndices(isy, &ixy));
556: for (PetscInt j = 0; j < sm; j++) {
557: PetscInt ix = ixy[j];
558: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
559: (*x)[ix - rstart] = y[j];
560: }
561: PetscCall(ISRestoreIndices(isy, &ixy));
562: PetscCall(VecRestoreArrayRead(subvec, &y));
563: }
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
568: {
569: Vec_Nest *bx = (Vec_Nest *)X->data;
570: PetscInt i, m, rstart, rend;
572: PetscFunctionBegin;
573: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
574: PetscCall(VecGetLocalSize(X, &m));
575: for (i = 0; i < bx->nb; i++) {
576: Vec subvec = bx->v[i];
577: IS isy = bx->is[i];
578: PetscInt sm;
579: const PetscInt *ixy;
580: PetscScalar *y;
581: PetscCall(VecGetLocalSize(subvec, &sm));
582: PetscCall(VecGetArray(subvec, &y));
583: PetscCall(ISGetIndices(isy, &ixy));
584: for (PetscInt j = 0; j < sm; j++) {
585: PetscInt ix = ixy[j];
586: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
587: y[j] = (*x)[ix - rstart];
588: }
589: PetscCall(ISRestoreIndices(isy, &ixy));
590: PetscCall(VecRestoreArray(subvec, &y));
591: }
592: PetscCall(PetscFree(*x));
593: PetscCall(PetscObjectStateIncrease((PetscObject)X));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
598: {
599: PetscFunctionBegin;
600: PetscCall(PetscFree(*x));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
605: {
606: PetscFunctionBegin;
607: PetscCheck(nx <= 0, PetscObjectComm((PetscObject)*X), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
611: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
612: {
613: Vec *ww;
614: IS *wis;
615: Vec_Nest *bv = (Vec_Nest *)v->data;
616: PetscInt i;
617: PetscMPIInt size;
619: PetscFunctionBegin;
620: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
621: if (size == 1) {
622: PetscCall(VecDuplicate(v, w));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
625: PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
626: for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
627: for (i = 0; i < bv->nb; i++) {
628: PetscInt off;
630: PetscCall(VecGetOwnershipRange(v, &off, NULL));
631: PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
632: PetscCall(ISShift(wis[i], -off, wis[i]));
633: }
634: PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
635: for (i = 0; i < bv->nb; i++) {
636: PetscCall(VecDestroy(&ww[i]));
637: PetscCall(ISDestroy(&wis[i]));
638: }
639: PetscCall(PetscFree2(ww, wis));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
644: {
645: Vec_Nest *bv = (Vec_Nest *)v->data;
646: Vec_Nest *bw = (Vec_Nest *)w->data;
647: PetscInt i;
649: PetscFunctionBegin;
650: PetscCheckTypeName(v, VECNEST);
651: PetscCheckTypeName(w, VECNEST);
652: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
653: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
658: {
659: Vec_Nest *bv = (Vec_Nest *)v->data;
660: Vec_Nest *bw = (Vec_Nest *)w->data;
661: PetscInt i;
663: PetscFunctionBegin;
664: PetscCheckTypeName(v, VECNEST);
665: PetscCheckTypeName(w, VECNEST);
666: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
667: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
668: PetscCall(PetscObjectStateIncrease((PetscObject)v));
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode VecGetLocalVectorRead_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: PetscCheckTypeName(v, VECNEST);
680: PetscCheckTypeName(w, VECNEST);
681: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
682: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
687: {
688: Vec_Nest *bv = (Vec_Nest *)v->data;
689: Vec_Nest *bw = (Vec_Nest *)w->data;
690: PetscInt i;
692: PetscFunctionBegin;
693: PetscCheckTypeName(v, VECNEST);
694: PetscCheckTypeName(w, VECNEST);
695: PetscCheck(bv->nb == bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
696: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
701: {
702: Vec_Nest *bv = (Vec_Nest *)v->data;
704: PetscFunctionBegin;
705: for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: 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)
710: {
711: Vec_Nest *bu = (Vec_Nest *)U->data;
712: Vec_Nest *by = (Vec_Nest *)Y->data;
713: Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
714: Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
715: Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;
717: PetscFunctionBegin;
718: VecNestCheckCompatible2(U, 1, Y, 2);
719: if (E) VecNestCheckCompatible2(U, 1, E, 3);
720: if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
721: if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
722: *norm = 0.0;
723: *norma = 0.0;
724: *normr = 0.0;
725: *norm_loc = 0;
726: *norma_loc = 0;
727: *normr_loc = 0;
728: for (PetscInt i = 0; i < bu->nb; i++) {
729: PetscReal n, na, nr;
730: PetscInt n_loc, na_loc, nr_loc;
732: 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));
733: if (wnormtype == NORM_INFINITY) {
734: *norm = PetscMax(*norm, n);
735: *norma = PetscMax(*norma, na);
736: *normr = PetscMax(*normr, nr);
737: } else {
738: *norm += PetscSqr(n);
739: *norma += PetscSqr(na);
740: *normr += PetscSqr(nr);
741: }
742: *norm_loc += n_loc;
743: *norma_loc += na_loc;
744: *normr_loc += nr_loc;
745: }
746: if (wnormtype == NORM_2) {
747: *norm = PetscSqrtReal(*norm);
748: *norma = PetscSqrtReal(*norma);
749: *normr = PetscSqrtReal(*normr);
750: }
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
755: {
756: PetscFunctionBegin;
757: ops->duplicate = VecDuplicate_Nest;
758: ops->duplicatevecs = VecDuplicateVecs_Default;
759: ops->destroyvecs = VecDestroyVecs_Default;
760: ops->dot = VecDot_Nest;
761: ops->mdot = VecMDot_Nest;
762: ops->norm = VecNorm_Nest;
763: ops->tdot = VecTDot_Nest;
764: ops->mtdot = VecMTDot_Nest;
765: ops->scale = VecScale_Nest;
766: ops->copy = VecCopy_Nest;
767: ops->set = VecSet_Nest;
768: ops->swap = VecSwap_Nest;
769: ops->axpy = VecAXPY_Nest;
770: ops->axpby = VecAXPBY_Nest;
771: ops->maxpy = VecMAXPY_Nest;
772: ops->aypx = VecAYPX_Nest;
773: ops->waxpy = VecWAXPY_Nest;
774: ops->axpbypcz = NULL;
775: ops->pointwisemult = VecPointwiseMult_Nest;
776: ops->pointwisedivide = VecPointwiseDivide_Nest;
777: ops->setvalues = NULL;
778: ops->assemblybegin = VecAssemblyBegin_Nest;
779: ops->assemblyend = VecAssemblyEnd_Nest;
780: ops->getarray = VecGetArray_Nest;
781: ops->getsize = VecGetSize_Nest;
782: ops->getlocalsize = VecGetLocalSize_Nest;
783: ops->restorearray = VecRestoreArray_Nest;
784: ops->restorearrayread = VecRestoreArrayRead_Nest;
785: ops->max = VecMax_Nest;
786: ops->min = VecMin_Nest;
787: ops->setrandom = NULL;
788: ops->setoption = NULL;
789: ops->setvaluesblocked = NULL;
790: ops->destroy = VecDestroy_Nest;
791: ops->view = VecView_Nest;
792: ops->placearray = NULL;
793: ops->replacearray = NULL;
794: ops->dot_local = VecDot_Nest;
795: ops->tdot_local = VecTDot_Nest;
796: ops->norm_local = VecNorm_Nest;
797: ops->mdot_local = VecMDot_Nest;
798: ops->mtdot_local = VecMTDot_Nest;
799: ops->load = NULL;
800: ops->reciprocal = VecReciprocal_Nest;
801: ops->conjugate = VecConjugate_Nest;
802: ops->setlocaltoglobalmapping = NULL;
803: ops->resetarray = NULL;
804: ops->setfromoptions = NULL;
805: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
806: ops->load = NULL;
807: ops->pointwisemax = NULL;
808: ops->pointwisemaxabs = NULL;
809: ops->pointwisemin = NULL;
810: ops->getvalues = NULL;
811: ops->sqrt = NULL;
812: ops->abs = NULL;
813: ops->exp = NULL;
814: ops->shift = NULL;
815: ops->create = NULL;
816: ops->stridegather = NULL;
817: ops->stridescatter = NULL;
818: ops->dotnorm2 = VecDotNorm2_Nest;
819: ops->getsubvector = VecGetSubVector_Nest;
820: ops->restoresubvector = VecRestoreSubVector_Nest;
821: ops->axpbypcz = VecAXPBYPCZ_Nest;
822: ops->concatenate = VecConcatenate_Nest;
823: ops->createlocalvector = VecCreateLocalVector_Nest;
824: ops->getlocalvector = VecGetLocalVector_Nest;
825: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
826: ops->restorelocalvector = VecRestoreLocalVector_Nest;
827: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
828: ops->setrandom = VecSetRandom_Nest;
829: ops->errorwnorm = VecErrorWeightedNorms_Nest;
830: ops->setstdbasis = NULL;
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
835: {
836: Vec_Nest *b = (Vec_Nest *)x->data;
837: PetscInt i;
838: PetscInt row;
840: PetscFunctionBegin;
841: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
842: for (i = 0; i < m; i++) {
843: row = idxm[i];
844: PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
845: vec[i] = b->v[row];
846: }
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
851: {
852: PetscFunctionBegin;
853: PetscCall(PetscObjectStateIncrease((PetscObject)X));
854: PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: /*@
859: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
861: Not Collective
863: Input Parameters:
864: + X - nest vector
865: - idxm - index of the vector within the nest
867: Output Parameter:
868: . sx - vector at index `idxm` within the nest
870: Level: developer
872: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
873: @*/
874: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
875: {
876: PetscFunctionBegin;
877: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
882: {
883: Vec_Nest *b = (Vec_Nest *)X->data;
885: PetscFunctionBegin;
886: PetscCall(PetscObjectStateIncrease((PetscObject)X));
887: if (N) *N = b->nb;
888: if (sx) *sx = b->v;
889: PetscFunctionReturn(PETSC_SUCCESS);
890: }
892: /*@
893: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
895: Not Collective
897: Input Parameter:
898: . X - nest vector
900: Output Parameters:
901: + N - number of nested vecs
902: - sx - array of vectors, can pass in `NULL`
904: Level: developer
906: Note:
907: The user should not free the array `sx`.
909: Fortran Notes:
910: The caller must allocate the array to hold the subvectors and pass it in.
912: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestGetSubVecsRead()`
913: @*/
914: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
915: {
916: PetscFunctionBegin;
917: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
922: {
923: Vec_Nest *bx = (Vec_Nest *)X->data;
924: PetscInt i, offset = 0, n = 0, bs;
925: IS is;
926: PetscBool issame = PETSC_FALSE;
927: PetscInt N = 0;
929: PetscFunctionBegin;
930: /* check if idxm < bx->nb */
931: 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);
933: PetscCall(PetscObjectReference((PetscObject)x));
934: PetscCall(VecDestroy(&bx->v[idxm]));
935: bx->v[idxm] = x;
937: /* check if we need to update the IS for the block */
938: offset = X->map->rstart;
939: for (i = 0; i < idxm; i++) {
940: n = 0;
941: PetscCall(VecGetLocalSize(bx->v[i], &n));
942: offset += n;
943: }
945: /* get the local size and block size */
946: PetscCall(VecGetLocalSize(x, &n));
947: PetscCall(VecGetBlockSize(x, &bs));
949: /* create the new IS */
950: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
951: PetscCall(ISSetBlockSize(is, bs));
953: /* check if they are equal */
954: PetscCall(ISEqual(is, bx->is[idxm], &issame));
956: if (!issame) {
957: /* The IS of given vector has a different layout compared to the existing block vector.
958: Destroy the existing reference and update the IS. */
959: PetscCall(ISDestroy(&bx->is[idxm]));
960: PetscCall(ISDuplicate(is, &bx->is[idxm]));
961: PetscCall(ISCopy(is, bx->is[idxm]));
963: offset += n;
964: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
965: for (i = idxm + 1; i < bx->nb; i++) {
966: /* get the local size and block size */
967: PetscCall(VecGetLocalSize(bx->v[i], &n));
968: PetscCall(VecGetBlockSize(bx->v[i], &bs));
970: /* destroy the old and create the new IS */
971: PetscCall(ISDestroy(&bx->is[i]));
972: PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
973: PetscCall(ISSetBlockSize(bx->is[i], bs));
975: offset += n;
976: }
978: n = 0;
979: PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
980: PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
981: PetscCall(PetscLayoutSetSize(X->map, N));
982: PetscCall(PetscLayoutSetLocalSize(X->map, n));
983: }
985: PetscCall(ISDestroy(&is));
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
990: {
991: PetscFunctionBegin;
992: PetscCall(PetscObjectStateIncrease((PetscObject)X));
993: PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: /*@
998: VecNestGetSubVecsRead - Access the subvecs of a `VECNEST` vector for read-only access
1000: Logically collective
1002: Input Parameter:
1003: . X - nest vector
1005: Output Parameters:
1006: + N - number of nested vecs
1007: - sx - array of read-locked vectors
1009: Level: advanced
1011: Notes:
1012: Each of the subvecs will be read locked (`VecLockReadPush()`), which is a logically collective operation.
1013: When access is complete, you must call `VecNestRestoreSubVecsRead()` to release the locks.
1015: Developer Note:
1016: This function does not increase the state of `X` (`PetscObjectStateIncrease()`).
1018: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestRestoreSubVecsRead()`
1019: @*/
1020: PetscErrorCode VecNestGetSubVecsRead(Vec X, PetscInt *N, Vec *sx[])
1021: {
1022: PetscFunctionBegin;
1023: PetscUseMethod(X, "VecNestGetSubVecsRead_C", (Vec, PetscInt *, Vec **), (X, N, sx));
1024: PetscFunctionReturn(PETSC_SUCCESS);
1025: }
1027: static PetscErrorCode VecNestGetSubVecsRead_Nest(Vec X, PetscInt *N, Vec **sx)
1028: {
1029: Vec_Nest *b = (Vec_Nest *)X->data;
1031: PetscFunctionBegin;
1032: for (PetscInt i = 0; i < b->nb; i++) PetscCall(VecLockReadPush(b->v[i]));
1033: if (N) *N = b->nb;
1034: if (sx) *sx = b->v;
1035: PetscFunctionReturn(PETSC_SUCCESS);
1036: }
1038: /*@
1039: VecNestRestoreSubVecsRead - Restore access the subvecs of a `VECNEST` vector obtained with `VecNestGetSubVecsRead()`
1041: Logically collective
1043: Input Parameters:
1044: + X - nest vector
1045: . N - number of nested vecs
1046: - sx - array of read-locked vectors
1048: Level: advanced
1050: Note:
1051: The same arguments to `VecNestGetSubVecsRead()` should be the argument to `VecNestRestoreSubVecsRead()`.
1053: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`, `VecNestGetSubVecsRead()`
1054: @*/
1055: PetscErrorCode VecNestRestoreSubVecsRead(Vec X, PetscInt *N, Vec *sx[])
1056: {
1057: PetscFunctionBegin;
1058: PetscUseMethod(X, "VecNestRestoreSubVecsRead_C", (Vec, PetscInt *, Vec **), (X, N, sx));
1059: PetscFunctionReturn(PETSC_SUCCESS);
1060: }
1062: static PetscErrorCode VecNestRestoreSubVecsRead_Nest(Vec X, PetscInt *N, Vec **sx)
1063: {
1064: Vec_Nest *b = (Vec_Nest *)X->data;
1066: PetscFunctionBegin;
1067: if (N) *N = 0;
1068: if (sx) {
1069: PetscCheck(*sx == b->v, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONG, "Restoring incorrect array of vectors");
1070: *sx = NULL;
1071: }
1072: for (PetscInt i = 0; i < b->nb; i++) PetscCall(VecLockReadPop(b->v[i]));
1073: PetscFunctionReturn(PETSC_SUCCESS);
1074: }
1076: /*@
1077: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
1079: Not Collective
1081: Input Parameters:
1082: + X - nest vector
1083: . idxm - index of the vector within the nest vector
1084: - sx - vector at index `idxm` within the nest vector
1086: Level: developer
1088: Notes:
1089: The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1091: The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.
1093: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
1094: @*/
1095: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
1096: {
1097: PetscFunctionBegin;
1098: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1103: {
1104: PetscInt i;
1106: PetscFunctionBegin;
1107: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1108: for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: /*@
1113: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1115: Not Collective
1117: Input Parameters:
1118: + X - nest vector
1119: . N - number of component vecs in `sx`
1120: . idxm - indices of component vectors that are to be replaced
1121: - sx - array of vectors
1123: Level: developer
1125: Notes:
1126: The components in the vector array `sx` do not have to be of the same size as corresponding
1127: components in `X`. The user can also free the array `sx` after the call.
1129: The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.
1131: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1132: @*/
1133: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
1134: {
1135: PetscFunctionBegin;
1136: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1137: PetscFunctionReturn(PETSC_SUCCESS);
1138: }
1140: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1141: {
1142: Vec_Nest *b = (Vec_Nest *)X->data;
1144: PetscFunctionBegin;
1145: *N = b->nb;
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: /*@
1150: VecNestGetSize - Returns the size of the nest vector.
1152: Not Collective
1154: Input Parameter:
1155: . X - nest vector
1157: Output Parameter:
1158: . N - number of nested vecs
1160: Level: developer
1162: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1163: @*/
1164: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1165: {
1166: PetscFunctionBegin;
1168: PetscAssertPointer(N, 2);
1169: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1170: PetscFunctionReturn(PETSC_SUCCESS);
1171: }
1173: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1174: {
1175: Vec_Nest *ctx = (Vec_Nest *)V->data;
1176: PetscInt i;
1178: PetscFunctionBegin;
1179: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1181: ctx->nb = nb;
1182: PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1184: /* Create space */
1185: PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1186: for (i = 0; i < ctx->nb; i++) {
1187: ctx->v[i] = x[i];
1188: PetscCall(PetscObjectReference((PetscObject)x[i]));
1189: /* Do not allocate memory for internal sub blocks */
1190: }
1192: PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1194: ctx->setup_called = PETSC_TRUE;
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1199: {
1200: Vec_Nest *ctx = (Vec_Nest *)V->data;
1201: PetscInt i, offset, m, n, M, N;
1203: PetscFunctionBegin;
1204: (void)nb;
1205: if (is) { /* Do some consistency checks and reference the is */
1206: offset = V->map->rstart;
1207: for (i = 0; i < ctx->nb; i++) {
1208: PetscCall(ISGetSize(is[i], &M));
1209: PetscCall(VecGetSize(ctx->v[i], &N));
1210: 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);
1211: PetscCall(ISGetLocalSize(is[i], &m));
1212: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1213: 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);
1214: PetscCall(PetscObjectReference((PetscObject)is[i]));
1215: ctx->is[i] = is[i];
1216: offset += n;
1217: }
1218: } else { /* Create a contiguous ISStride for each entry */
1219: offset = V->map->rstart;
1220: for (i = 0; i < ctx->nb; i++) {
1221: PetscInt bs;
1222: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1223: PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1224: PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1225: PetscCall(ISSetBlockSize(ctx->is[i], bs));
1226: offset += n;
1227: }
1228: }
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: /*MC
1233: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1235: Level: intermediate
1237: Notes:
1238: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1239: It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1241: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecCreateNest()`, `MatCreateNest()`
1242: M*/
1244: /*@
1245: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1247: Collective
1249: Input Parameters:
1250: + comm - Communicator for the new `Vec`
1251: . nb - number of nested blocks
1252: . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1253: - x - array of `nb` sub-vectors
1255: Output Parameter:
1256: . Y - new vector
1258: Level: advanced
1260: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1261: @*/
1262: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1263: {
1264: Vec V;
1265: Vec_Nest *s;
1266: PetscInt n, N;
1268: PetscFunctionBegin;
1269: PetscCall(VecCreate(comm, &V));
1270: PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1272: /* allocate and set pointer for implementation data */
1273: PetscCall(PetscNew(&s));
1274: V->data = (void *)s;
1275: s->setup_called = PETSC_FALSE;
1276: s->nb = -1;
1277: s->v = NULL;
1279: PetscCall(VecSetUp_Nest_Private(V, nb, x));
1281: n = N = 0;
1282: PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1283: PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1284: PetscCall(PetscLayoutSetSize(V->map, N));
1285: PetscCall(PetscLayoutSetLocalSize(V->map, n));
1286: PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1287: PetscCall(PetscLayoutSetUp(V->map));
1289: PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1291: PetscCall(VecNestSetOps_Private(V->ops));
1292: V->petscnative = PETSC_FALSE;
1294: /* expose block api's */
1295: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1296: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1297: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecsRead_C", VecNestGetSubVecsRead_Nest));
1298: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestRestoreSubVecsRead_C", VecNestRestoreSubVecsRead_Nest));
1299: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1300: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1301: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1303: *Y = V;
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }