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, "VecNestSetSubVec_C", NULL));
42: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL));
43: PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL));
45: PetscCall(PetscFree(vs));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
50: {
51: Vec_Nest *bx = (Vec_Nest *)x->data;
52: Vec_Nest *by = (Vec_Nest *)y->data;
53: PetscInt i;
55: PetscFunctionBegin;
56: VecNestCheckCompatible2(x, 1, y, 2);
57: for (i = 0; i < bx->nb; i++) PetscCall(VecCopy(bx->v[i], by->v[i]));
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
62: {
63: Vec_Nest *bx = (Vec_Nest *)x->data;
64: Vec Y;
65: Vec *sub;
66: PetscInt i;
68: PetscFunctionBegin;
69: PetscCall(PetscMalloc1(bx->nb, &sub));
70: for (i = 0; i < bx->nb; i++) PetscCall(VecDuplicate(bx->v[i], &sub[i]));
71: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y));
72: for (i = 0; i < bx->nb; i++) PetscCall(VecDestroy(&sub[i]));
73: PetscCall(PetscFree(sub));
74: *y = Y;
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
79: {
80: Vec_Nest *bx = (Vec_Nest *)x->data;
81: Vec_Nest *by = (Vec_Nest *)y->data;
82: PetscInt i, nr;
83: PetscScalar x_dot_y, _val;
85: PetscFunctionBegin;
86: VecNestCheckCompatible2(x, 1, y, 2);
87: nr = bx->nb;
88: _val = 0.0;
89: for (i = 0; i < nr; i++) {
90: PetscCall(VecDot(bx->v[i], by->v[i], &x_dot_y));
91: _val = _val + x_dot_y;
92: }
93: *val = _val;
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
98: {
99: Vec_Nest *bx = (Vec_Nest *)x->data;
100: Vec_Nest *by = (Vec_Nest *)y->data;
101: PetscInt i, nr;
102: PetscScalar x_dot_y, _val;
104: PetscFunctionBegin;
105: VecNestCheckCompatible2(x, 1, y, 2);
106: nr = bx->nb;
107: _val = 0.0;
108: for (i = 0; i < nr; i++) {
109: PetscCall(VecTDot(bx->v[i], by->v[i], &x_dot_y));
110: _val = _val + x_dot_y;
111: }
112: *val = _val;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
117: {
118: Vec_Nest *bx = (Vec_Nest *)x->data;
119: Vec_Nest *by = (Vec_Nest *)y->data;
120: PetscInt i, nr;
121: PetscScalar x_dot_y, _dp, _nm;
122: PetscReal norm2_y;
124: PetscFunctionBegin;
125: VecNestCheckCompatible2(x, 1, y, 2);
126: nr = bx->nb;
127: _dp = 0.0;
128: _nm = 0.0;
129: for (i = 0; i < nr; i++) {
130: PetscCall(VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y));
131: _dp += x_dot_y;
132: _nm += norm2_y;
133: }
134: *dp = _dp;
135: *nm = _nm;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
140: {
141: Vec_Nest *bx = (Vec_Nest *)x->data;
142: Vec_Nest *by = (Vec_Nest *)y->data;
143: PetscInt i, nr;
145: PetscFunctionBegin;
146: VecNestCheckCompatible2(y, 1, x, 3);
147: nr = bx->nb;
148: for (i = 0; i < nr; i++) PetscCall(VecAXPY(by->v[i], alpha, bx->v[i]));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
153: {
154: Vec_Nest *bx = (Vec_Nest *)x->data;
155: Vec_Nest *by = (Vec_Nest *)y->data;
156: PetscInt i, nr;
158: PetscFunctionBegin;
159: VecNestCheckCompatible2(y, 1, x, 3);
160: nr = bx->nb;
161: for (i = 0; i < nr; i++) PetscCall(VecAYPX(by->v[i], alpha, bx->v[i]));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
166: {
167: Vec_Nest *bx = (Vec_Nest *)x->data;
168: Vec_Nest *by = (Vec_Nest *)y->data;
169: PetscInt i, nr;
171: PetscFunctionBegin;
172: VecNestCheckCompatible2(y, 1, x, 4);
173: nr = bx->nb;
174: for (i = 0; i < nr; i++) PetscCall(VecAXPBY(by->v[i], alpha, beta, bx->v[i]));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
179: {
180: Vec_Nest *bx = (Vec_Nest *)x->data;
181: Vec_Nest *by = (Vec_Nest *)y->data;
182: Vec_Nest *bz = (Vec_Nest *)z->data;
183: PetscInt i, nr;
185: PetscFunctionBegin;
186: VecNestCheckCompatible3(z, 1, x, 5, y, 6);
187: nr = bx->nb;
188: for (i = 0; i < nr; i++) PetscCall(VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
193: {
194: Vec_Nest *bx = (Vec_Nest *)x->data;
195: PetscInt i, nr;
197: PetscFunctionBegin;
198: nr = bx->nb;
199: for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
204: {
205: Vec_Nest *bx = (Vec_Nest *)x->data;
206: Vec_Nest *by = (Vec_Nest *)y->data;
207: Vec_Nest *bw = (Vec_Nest *)w->data;
208: PetscInt i, nr;
210: PetscFunctionBegin;
211: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
212: nr = bx->nb;
213: for (i = 0; i < nr; i++) PetscCall(VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
218: {
219: Vec_Nest *bx = (Vec_Nest *)x->data;
220: Vec_Nest *by = (Vec_Nest *)y->data;
221: Vec_Nest *bw = (Vec_Nest *)w->data;
222: PetscInt i, nr;
224: PetscFunctionBegin;
225: VecNestCheckCompatible3(w, 1, x, 2, y, 3);
226: nr = bx->nb;
227: for (i = 0; i < nr; i++) PetscCall(VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: static PetscErrorCode VecReciprocal_Nest(Vec x)
232: {
233: Vec_Nest *bx = (Vec_Nest *)x->data;
234: PetscInt i, nr;
236: PetscFunctionBegin;
237: nr = bx->nb;
238: for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
243: {
244: Vec_Nest *bx = (Vec_Nest *)xin->data;
245: PetscInt i, nr;
246: PetscReal z_i;
247: PetscReal _z;
249: PetscFunctionBegin;
250: nr = bx->nb;
251: _z = 0.0;
253: if (type == NORM_2) {
254: PetscScalar dot;
255: PetscCall(VecDot(xin, xin, &dot));
256: _z = PetscAbsScalar(PetscSqrtScalar(dot));
257: } else if (type == NORM_1) {
258: for (i = 0; i < nr; i++) {
259: PetscCall(VecNorm(bx->v[i], type, &z_i));
260: _z = _z + z_i;
261: }
262: } else if (type == NORM_INFINITY) {
263: for (i = 0; i < nr; i++) {
264: PetscCall(VecNorm(bx->v[i], type, &z_i));
265: if (z_i > _z) _z = z_i;
266: }
267: }
269: *z = _z;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
274: {
275: PetscFunctionBegin;
276: /* TODO: implement proper MAXPY. For now, do axpy on each vector */
277: for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
282: {
283: PetscFunctionBegin;
284: /* TODO: implement proper MDOT. For now, do dot on each vector */
285: for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
290: {
291: PetscFunctionBegin;
292: /* TODO: implement proper MTDOT. For now, do tdot on each vector */
293: for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
298: {
299: Vec_Nest *bx = (Vec_Nest *)x->data;
300: PetscInt i, nr;
302: PetscFunctionBegin;
303: nr = bx->nb;
304: for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode VecConjugate_Nest(Vec x)
309: {
310: Vec_Nest *bx = (Vec_Nest *)x->data;
311: PetscInt j, nr;
313: PetscFunctionBegin;
314: nr = bx->nb;
315: for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
320: {
321: Vec_Nest *bx = (Vec_Nest *)x->data;
322: Vec_Nest *by = (Vec_Nest *)y->data;
323: PetscInt i, nr;
325: PetscFunctionBegin;
326: VecNestCheckCompatible2(x, 1, y, 2);
327: nr = bx->nb;
328: for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
333: {
334: Vec_Nest *bx = (Vec_Nest *)x->data;
335: Vec_Nest *by = (Vec_Nest *)y->data;
336: Vec_Nest *bw = (Vec_Nest *)w->data;
337: PetscInt i, nr;
339: PetscFunctionBegin;
340: VecNestCheckCompatible3(w, 1, x, 3, y, 4);
341: nr = bx->nb;
342: for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
347: {
348: PetscInt i, lp = -1, lw = -1;
349: PetscReal lv;
350: Vec_Nest *bx = (Vec_Nest *)x->data;
352: PetscFunctionBegin;
353: *v = PETSC_MAX_REAL;
354: for (i = 0; i < bx->nb; i++) {
355: PetscInt tp;
356: PetscCall(VecMin(bx->v[i], &tp, &lv));
357: if (lv < *v) {
358: *v = lv;
359: lw = i;
360: lp = tp;
361: }
362: }
363: if (p && lw > -1) {
364: PetscInt st, en;
365: const PetscInt *idxs;
367: *p = -1;
368: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
369: if (lp >= st && lp < en) {
370: PetscCall(ISGetIndices(bx->is[lw], &idxs));
371: *p = idxs[lp - st];
372: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
373: }
374: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
375: }
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
380: {
381: PetscInt i, lp = -1, lw = -1;
382: PetscReal lv;
383: Vec_Nest *bx = (Vec_Nest *)x->data;
385: PetscFunctionBegin;
386: *v = PETSC_MIN_REAL;
387: for (i = 0; i < bx->nb; i++) {
388: PetscInt tp;
389: PetscCall(VecMax(bx->v[i], &tp, &lv));
390: if (lv > *v) {
391: *v = lv;
392: lw = i;
393: lp = tp;
394: }
395: }
396: if (p && lw > -1) {
397: PetscInt st, en;
398: const PetscInt *idxs;
400: *p = -1;
401: PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
402: if (lp >= st && lp < en) {
403: PetscCall(ISGetIndices(bx->is[lw], &idxs));
404: *p = idxs[lp - st];
405: PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
406: }
407: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
408: }
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
413: {
414: Vec_Nest *bx = (Vec_Nest *)x->data;
415: PetscBool isascii;
416: PetscInt i;
418: PetscFunctionBegin;
419: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
420: if (isascii) {
421: PetscCall(PetscViewerASCIIPushTab(viewer));
422: PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ", structure: \n", bx->nb));
423: for (i = 0; i < bx->nb; i++) {
424: VecType type;
425: char name[256] = "", prefix[256] = "";
426: PetscInt NR;
428: PetscCall(VecGetSize(bx->v[i], &NR));
429: PetscCall(VecGetType(bx->v[i], &type));
430: if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
431: if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));
433: PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n", i, name, prefix, type, NR));
435: PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
436: PetscCall(VecView(bx->v[i], viewer));
437: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
438: }
439: PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
440: }
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
445: {
446: Vec_Nest *bx;
447: PetscInt size, i, nr;
448: PetscBool isnest;
450: PetscFunctionBegin;
451: PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
452: if (!isnest) {
453: /* Not nest */
454: if (globalsize) PetscCall(VecGetSize(x, &size));
455: else PetscCall(VecGetLocalSize(x, &size));
456: *L = *L + size;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /* Otherwise we have a nest */
461: bx = (Vec_Nest *)x->data;
462: nr = bx->nb;
464: /* now descend recursively */
465: for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /* Returns the sum of the global size of all the constituent vectors in the nest */
470: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
471: {
472: PetscFunctionBegin;
473: *N = x->map->N;
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /* Returns the sum of the local size of all the constituent vectors in the nest */
478: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
479: {
480: PetscFunctionBegin;
481: *n = x->map->n;
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
486: {
487: Vec_Nest *bx = (Vec_Nest *)x->data;
488: Vec_Nest *by = (Vec_Nest *)y->data;
489: PetscInt i, nr;
490: PetscReal local_max, m;
492: PetscFunctionBegin;
493: VecNestCheckCompatible2(x, 1, y, 2);
494: nr = bx->nb;
495: m = 0.0;
496: for (i = 0; i < nr; i++) {
497: PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
498: if (local_max > m) m = local_max;
499: }
500: *max = m;
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
505: {
506: Vec_Nest *bx = (Vec_Nest *)X->data;
507: PetscInt i;
509: PetscFunctionBegin;
510: *x = NULL;
511: for (i = 0; i < bx->nb; i++) {
512: PetscBool issame = PETSC_FALSE;
513: PetscCall(ISEqual(is, bx->is[i], &issame));
514: if (issame) {
515: *x = bx->v[i];
516: PetscCall(PetscObjectReference((PetscObject)*x));
517: break;
518: }
519: }
520: PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
525: {
526: PetscFunctionBegin;
527: PetscCall(PetscObjectStateIncrease((PetscObject)X));
528: PetscCall(VecDestroy(x));
529: PetscFunctionReturn(PETSC_SUCCESS);
530: }
532: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
533: {
534: Vec_Nest *bx = (Vec_Nest *)X->data;
535: PetscInt i, m, rstart, rend;
537: PetscFunctionBegin;
538: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
539: PetscCall(VecGetLocalSize(X, &m));
540: PetscCall(PetscMalloc1(m, x));
541: for (i = 0; i < bx->nb; i++) {
542: Vec subvec = bx->v[i];
543: IS isy = bx->is[i];
544: PetscInt j, sm;
545: const PetscInt *ixy;
546: const PetscScalar *y;
547: PetscCall(VecGetLocalSize(subvec, &sm));
548: PetscCall(VecGetArrayRead(subvec, &y));
549: PetscCall(ISGetIndices(isy, &ixy));
550: for (j = 0; j < sm; j++) {
551: PetscInt ix = ixy[j];
552: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
553: (*x)[ix - rstart] = y[j];
554: }
555: PetscCall(ISRestoreIndices(isy, &ixy));
556: PetscCall(VecRestoreArrayRead(subvec, &y));
557: }
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
562: {
563: Vec_Nest *bx = (Vec_Nest *)X->data;
564: PetscInt i, m, rstart, rend;
566: PetscFunctionBegin;
567: PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
568: PetscCall(VecGetLocalSize(X, &m));
569: for (i = 0; i < bx->nb; i++) {
570: Vec subvec = bx->v[i];
571: IS isy = bx->is[i];
572: PetscInt j, sm;
573: const PetscInt *ixy;
574: PetscScalar *y;
575: PetscCall(VecGetLocalSize(subvec, &sm));
576: PetscCall(VecGetArray(subvec, &y));
577: PetscCall(ISGetIndices(isy, &ixy));
578: for (j = 0; j < sm; j++) {
579: PetscInt ix = ixy[j];
580: PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
581: y[j] = (*x)[ix - rstart];
582: }
583: PetscCall(ISRestoreIndices(isy, &ixy));
584: PetscCall(VecRestoreArray(subvec, &y));
585: }
586: PetscCall(PetscFree(*x));
587: PetscCall(PetscObjectStateIncrease((PetscObject)X));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
592: {
593: PetscFunctionBegin;
594: PetscCall(PetscFree(*x));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
599: {
600: PetscFunctionBegin;
601: PetscCheck(nx <= 0, PetscObjectComm((PetscObject)*X), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
606: {
607: Vec *ww;
608: IS *wis;
609: Vec_Nest *bv = (Vec_Nest *)v->data;
610: PetscInt i;
611: PetscMPIInt size;
613: PetscFunctionBegin;
614: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
615: if (size == 1) {
616: PetscCall(VecDuplicate(v, w));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
619: PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
620: for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
621: for (i = 0; i < bv->nb; i++) {
622: PetscInt off;
624: PetscCall(VecGetOwnershipRange(v, &off, NULL));
625: PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
626: PetscCall(ISShift(wis[i], -off, wis[i]));
627: }
628: PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
629: for (i = 0; i < bv->nb; i++) {
630: PetscCall(VecDestroy(&ww[i]));
631: PetscCall(ISDestroy(&wis[i]));
632: }
633: PetscCall(PetscFree2(ww, wis));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
638: {
639: Vec_Nest *bv = (Vec_Nest *)v->data;
640: Vec_Nest *bw = (Vec_Nest *)w->data;
641: PetscInt i;
643: PetscFunctionBegin;
644: PetscCheckSameType(v, 1, w, 2);
645: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
646: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
651: {
652: Vec_Nest *bv = (Vec_Nest *)v->data;
653: Vec_Nest *bw = (Vec_Nest *)w->data;
654: PetscInt i;
656: PetscFunctionBegin;
657: PetscCheckSameType(v, 1, w, 2);
658: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
659: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
660: PetscCall(PetscObjectStateIncrease((PetscObject)v));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
665: {
666: Vec_Nest *bv = (Vec_Nest *)v->data;
667: Vec_Nest *bw = (Vec_Nest *)w->data;
668: PetscInt i;
670: PetscFunctionBegin;
671: PetscCheckSameType(v, 1, w, 2);
672: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
673: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
678: {
679: Vec_Nest *bv = (Vec_Nest *)v->data;
680: Vec_Nest *bw = (Vec_Nest *)w->data;
681: PetscInt i;
683: PetscFunctionBegin;
684: PetscCheckSameType(v, 1, w, 2);
685: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
686: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
691: {
692: Vec_Nest *bv = (Vec_Nest *)v->data;
694: PetscFunctionBegin;
695: for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: 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)
700: {
701: Vec_Nest *bu = (Vec_Nest *)U->data;
702: Vec_Nest *by = (Vec_Nest *)Y->data;
703: Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
704: Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
705: Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;
707: PetscFunctionBegin;
708: VecNestCheckCompatible2(U, 1, Y, 2);
709: if (E) VecNestCheckCompatible2(U, 1, E, 3);
710: if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
711: if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
712: *norm = 0.0;
713: *norma = 0.0;
714: *normr = 0.0;
715: *norm_loc = 0;
716: *norma_loc = 0;
717: *normr_loc = 0;
718: for (PetscInt i = 0; i < bu->nb; i++) {
719: PetscReal n, na, nr;
720: PetscInt n_loc, na_loc, nr_loc;
722: 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));
723: if (wnormtype == NORM_INFINITY) {
724: *norm = PetscMax(*norm, n);
725: *norma = PetscMax(*norma, na);
726: *normr = PetscMax(*normr, nr);
727: } else {
728: *norm += PetscSqr(n);
729: *norma += PetscSqr(na);
730: *normr += PetscSqr(nr);
731: }
732: *norm_loc += n_loc;
733: *norma_loc += na_loc;
734: *normr_loc += nr_loc;
735: }
736: if (wnormtype == NORM_2) {
737: *norm = PetscSqrtReal(*norm);
738: *norma = PetscSqrtReal(*norma);
739: *normr = PetscSqrtReal(*normr);
740: }
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
745: {
746: PetscFunctionBegin;
747: ops->duplicate = VecDuplicate_Nest;
748: ops->duplicatevecs = VecDuplicateVecs_Default;
749: ops->destroyvecs = VecDestroyVecs_Default;
750: ops->dot = VecDot_Nest;
751: ops->mdot = VecMDot_Nest;
752: ops->norm = VecNorm_Nest;
753: ops->tdot = VecTDot_Nest;
754: ops->mtdot = VecMTDot_Nest;
755: ops->scale = VecScale_Nest;
756: ops->copy = VecCopy_Nest;
757: ops->set = VecSet_Nest;
758: ops->swap = VecSwap_Nest;
759: ops->axpy = VecAXPY_Nest;
760: ops->axpby = VecAXPBY_Nest;
761: ops->maxpy = VecMAXPY_Nest;
762: ops->aypx = VecAYPX_Nest;
763: ops->waxpy = VecWAXPY_Nest;
764: ops->axpbypcz = NULL;
765: ops->pointwisemult = VecPointwiseMult_Nest;
766: ops->pointwisedivide = VecPointwiseDivide_Nest;
767: ops->setvalues = NULL;
768: ops->assemblybegin = VecAssemblyBegin_Nest;
769: ops->assemblyend = VecAssemblyEnd_Nest;
770: ops->getarray = VecGetArray_Nest;
771: ops->getsize = VecGetSize_Nest;
772: ops->getlocalsize = VecGetLocalSize_Nest;
773: ops->restorearray = VecRestoreArray_Nest;
774: ops->restorearrayread = VecRestoreArrayRead_Nest;
775: ops->max = VecMax_Nest;
776: ops->min = VecMin_Nest;
777: ops->setrandom = NULL;
778: ops->setoption = NULL;
779: ops->setvaluesblocked = NULL;
780: ops->destroy = VecDestroy_Nest;
781: ops->view = VecView_Nest;
782: ops->placearray = NULL;
783: ops->replacearray = NULL;
784: ops->dot_local = VecDot_Nest;
785: ops->tdot_local = VecTDot_Nest;
786: ops->norm_local = VecNorm_Nest;
787: ops->mdot_local = VecMDot_Nest;
788: ops->mtdot_local = VecMTDot_Nest;
789: ops->load = NULL;
790: ops->reciprocal = VecReciprocal_Nest;
791: ops->conjugate = VecConjugate_Nest;
792: ops->setlocaltoglobalmapping = NULL;
793: ops->setvalueslocal = NULL;
794: ops->resetarray = NULL;
795: ops->setfromoptions = NULL;
796: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
797: ops->load = NULL;
798: ops->pointwisemax = NULL;
799: ops->pointwisemaxabs = NULL;
800: ops->pointwisemin = NULL;
801: ops->getvalues = NULL;
802: ops->sqrt = NULL;
803: ops->abs = NULL;
804: ops->exp = NULL;
805: ops->shift = NULL;
806: ops->create = NULL;
807: ops->stridegather = NULL;
808: ops->stridescatter = NULL;
809: ops->dotnorm2 = VecDotNorm2_Nest;
810: ops->getsubvector = VecGetSubVector_Nest;
811: ops->restoresubvector = VecRestoreSubVector_Nest;
812: ops->axpbypcz = VecAXPBYPCZ_Nest;
813: ops->concatenate = VecConcatenate_Nest;
814: ops->createlocalvector = VecCreateLocalVector_Nest;
815: ops->getlocalvector = VecGetLocalVector_Nest;
816: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
817: ops->restorelocalvector = VecRestoreLocalVector_Nest;
818: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
819: ops->setrandom = VecSetRandom_Nest;
820: ops->errorwnorm = VecErrorWeightedNorms_Nest;
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
824: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
825: {
826: Vec_Nest *b = (Vec_Nest *)x->data;
827: PetscInt i;
828: PetscInt row;
830: PetscFunctionBegin;
831: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
832: for (i = 0; i < m; i++) {
833: row = idxm[i];
834: PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
835: vec[i] = b->v[row];
836: }
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
841: {
842: PetscFunctionBegin;
843: PetscCall(PetscObjectStateIncrease((PetscObject)X));
844: PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: /*@
849: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
851: Not Collective
853: Input Parameters:
854: + X - nest vector
855: - idxm - index of the vector within the nest
857: Output Parameter:
858: . sx - vector at index `idxm` within the nest
860: Level: developer
862: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
863: @*/
864: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
865: {
866: PetscFunctionBegin;
867: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
872: {
873: Vec_Nest *b = (Vec_Nest *)X->data;
875: PetscFunctionBegin;
876: PetscCall(PetscObjectStateIncrease((PetscObject)X));
877: if (N) *N = b->nb;
878: if (sx) *sx = b->v;
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: /*@C
883: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
885: Not Collective
887: Input Parameter:
888: . X - nest vector
890: Output Parameters:
891: + N - number of nested vecs
892: - sx - array of vectors, can pass in `NULL`
894: Level: developer
896: Note:
897: The user should not free the array `sx`.
899: Fortran Notes:
900: The caller must allocate the array to hold the subvectors and pass it in.
902: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
903: @*/
904: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
905: {
906: PetscFunctionBegin;
907: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
912: {
913: Vec_Nest *bx = (Vec_Nest *)X->data;
914: PetscInt i, offset = 0, n = 0, bs;
915: IS is;
916: PetscBool issame = PETSC_FALSE;
917: PetscInt N = 0;
919: PetscFunctionBegin;
920: /* check if idxm < bx->nb */
921: 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);
923: PetscCall(PetscObjectReference((PetscObject)x));
924: PetscCall(VecDestroy(&bx->v[idxm]));
925: bx->v[idxm] = x;
927: /* check if we need to update the IS for the block */
928: offset = X->map->rstart;
929: for (i = 0; i < idxm; i++) {
930: n = 0;
931: PetscCall(VecGetLocalSize(bx->v[i], &n));
932: offset += n;
933: }
935: /* get the local size and block size */
936: PetscCall(VecGetLocalSize(x, &n));
937: PetscCall(VecGetBlockSize(x, &bs));
939: /* create the new IS */
940: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
941: PetscCall(ISSetBlockSize(is, bs));
943: /* check if they are equal */
944: PetscCall(ISEqual(is, bx->is[idxm], &issame));
946: if (!issame) {
947: /* The IS of given vector has a different layout compared to the existing block vector.
948: Destroy the existing reference and update the IS. */
949: PetscCall(ISDestroy(&bx->is[idxm]));
950: PetscCall(ISDuplicate(is, &bx->is[idxm]));
951: PetscCall(ISCopy(is, bx->is[idxm]));
953: offset += n;
954: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
955: for (i = idxm + 1; i < bx->nb; i++) {
956: /* get the local size and block size */
957: PetscCall(VecGetLocalSize(bx->v[i], &n));
958: PetscCall(VecGetBlockSize(bx->v[i], &bs));
960: /* destroy the old and create the new IS */
961: PetscCall(ISDestroy(&bx->is[i]));
962: PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
963: PetscCall(ISSetBlockSize(bx->is[i], bs));
965: offset += n;
966: }
968: n = 0;
969: PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
970: PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
971: PetscCall(PetscLayoutSetSize(X->map, N));
972: PetscCall(PetscLayoutSetLocalSize(X->map, n));
973: }
975: PetscCall(ISDestroy(&is));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
980: {
981: PetscFunctionBegin;
982: PetscCall(PetscObjectStateIncrease((PetscObject)X));
983: PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: /*@
988: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
990: Not Collective
992: Input Parameters:
993: + X - nest vector
994: . idxm - index of the vector within the nest vector
995: - sx - vector at index `idxm` within the nest vector
997: Level: developer
999: Notes:
1000: The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1002: The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.
1004: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
1005: @*/
1006: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
1007: {
1008: PetscFunctionBegin;
1009: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1014: {
1015: PetscInt i;
1017: PetscFunctionBegin;
1018: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1019: for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1020: PetscFunctionReturn(PETSC_SUCCESS);
1021: }
1023: /*@
1024: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1026: Not Collective
1028: Input Parameters:
1029: + X - nest vector
1030: . N - number of component vecs in `sx`
1031: . idxm - indices of component vectors that are to be replaced
1032: - sx - array of vectors
1034: Level: developer
1036: Notes:
1037: The components in the vector array `sx` do not have to be of the same size as corresponding
1038: components in `X`. The user can also free the array `sx` after the call.
1040: The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.
1042: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1043: @*/
1044: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
1045: {
1046: PetscFunctionBegin;
1047: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1052: {
1053: Vec_Nest *b = (Vec_Nest *)X->data;
1055: PetscFunctionBegin;
1056: *N = b->nb;
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: /*@
1061: VecNestGetSize - Returns the size of the nest vector.
1063: Not Collective
1065: Input Parameter:
1066: . X - nest vector
1068: Output Parameter:
1069: . N - number of nested vecs
1071: Level: developer
1073: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1074: @*/
1075: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1076: {
1077: PetscFunctionBegin;
1079: PetscAssertPointer(N, 2);
1080: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1085: {
1086: Vec_Nest *ctx = (Vec_Nest *)V->data;
1087: PetscInt i;
1089: PetscFunctionBegin;
1090: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1092: ctx->nb = nb;
1093: PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1095: /* Create space */
1096: PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1097: for (i = 0; i < ctx->nb; i++) {
1098: ctx->v[i] = x[i];
1099: PetscCall(PetscObjectReference((PetscObject)x[i]));
1100: /* Do not allocate memory for internal sub blocks */
1101: }
1103: PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1105: ctx->setup_called = PETSC_TRUE;
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1110: {
1111: Vec_Nest *ctx = (Vec_Nest *)V->data;
1112: PetscInt i, offset, m, n, M, N;
1114: PetscFunctionBegin;
1115: (void)nb;
1116: if (is) { /* Do some consistency checks and reference the is */
1117: offset = V->map->rstart;
1118: for (i = 0; i < ctx->nb; i++) {
1119: PetscCall(ISGetSize(is[i], &M));
1120: PetscCall(VecGetSize(ctx->v[i], &N));
1121: 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);
1122: PetscCall(ISGetLocalSize(is[i], &m));
1123: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1124: 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);
1125: PetscCall(PetscObjectReference((PetscObject)is[i]));
1126: ctx->is[i] = is[i];
1127: offset += n;
1128: }
1129: } else { /* Create a contiguous ISStride for each entry */
1130: offset = V->map->rstart;
1131: for (i = 0; i < ctx->nb; i++) {
1132: PetscInt bs;
1133: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1134: PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1135: PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1136: PetscCall(ISSetBlockSize(ctx->is[i], bs));
1137: offset += n;
1138: }
1139: }
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: /*MC
1144: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1146: Level: intermediate
1148: Notes:
1149: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1150: It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1152: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1153: M*/
1155: /*@C
1156: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1158: Collective
1160: Input Parameters:
1161: + comm - Communicator for the new `Vec`
1162: . nb - number of nested blocks
1163: . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1164: - x - array of `nb` sub-vectors
1166: Output Parameter:
1167: . Y - new vector
1169: Level: advanced
1171: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1172: @*/
1173: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1174: {
1175: Vec V;
1176: Vec_Nest *s;
1177: PetscInt n, N;
1179: PetscFunctionBegin;
1180: PetscCall(VecCreate(comm, &V));
1181: PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1183: /* allocate and set pointer for implementation data */
1184: PetscCall(PetscNew(&s));
1185: V->data = (void *)s;
1186: s->setup_called = PETSC_FALSE;
1187: s->nb = -1;
1188: s->v = NULL;
1190: PetscCall(VecSetUp_Nest_Private(V, nb, x));
1192: n = N = 0;
1193: PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1194: PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1195: PetscCall(PetscLayoutSetSize(V->map, N));
1196: PetscCall(PetscLayoutSetLocalSize(V->map, n));
1197: PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1198: PetscCall(PetscLayoutSetUp(V->map));
1200: PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1202: PetscCall(VecNestSetOps_Private(V->ops));
1203: V->petscnative = PETSC_FALSE;
1205: /* expose block api's */
1206: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1207: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1208: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1209: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1210: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1212: *Y = V;
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }