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;
612: PetscFunctionBegin;
613: PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
614: for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
615: for (i = 0; i < bv->nb; i++) {
616: PetscInt off;
618: PetscCall(VecGetOwnershipRange(v, &off, NULL));
619: PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
620: PetscCall(ISShift(wis[i], -off, wis[i]));
621: }
622: PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
623: for (i = 0; i < bv->nb; i++) {
624: PetscCall(VecDestroy(&ww[i]));
625: PetscCall(ISDestroy(&wis[i]));
626: }
627: PetscCall(PetscFree2(ww, wis));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
632: {
633: Vec_Nest *bv = (Vec_Nest *)v->data;
634: Vec_Nest *bw = (Vec_Nest *)w->data;
635: PetscInt i;
637: PetscFunctionBegin;
638: PetscCheckSameType(v, 1, w, 2);
639: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
640: for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
645: {
646: Vec_Nest *bv = (Vec_Nest *)v->data;
647: Vec_Nest *bw = (Vec_Nest *)w->data;
648: PetscInt i;
650: PetscFunctionBegin;
651: PetscCheckSameType(v, 1, w, 2);
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(VecRestoreLocalVector(bv->v[i], bw->v[i]));
654: PetscCall(PetscObjectStateIncrease((PetscObject)v));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
659: {
660: Vec_Nest *bv = (Vec_Nest *)v->data;
661: Vec_Nest *bw = (Vec_Nest *)w->data;
662: PetscInt i;
664: PetscFunctionBegin;
665: PetscCheckSameType(v, 1, w, 2);
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(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
672: {
673: Vec_Nest *bv = (Vec_Nest *)v->data;
674: Vec_Nest *bw = (Vec_Nest *)w->data;
675: PetscInt i;
677: PetscFunctionBegin;
678: PetscCheckSameType(v, 1, w, 2);
679: PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
680: for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
685: {
686: Vec_Nest *bv = (Vec_Nest *)v->data;
688: PetscFunctionBegin;
689: for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: 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)
694: {
695: Vec_Nest *bu = (Vec_Nest *)U->data;
696: Vec_Nest *by = (Vec_Nest *)Y->data;
697: Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
698: Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
699: Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;
701: PetscFunctionBegin;
702: VecNestCheckCompatible2(U, 1, Y, 2);
703: if (E) VecNestCheckCompatible2(U, 1, E, 3);
704: if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
705: if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
706: *norm = 0.0;
707: *norma = 0.0;
708: *normr = 0.0;
709: *norm_loc = 0;
710: *norma_loc = 0;
711: *normr_loc = 0;
712: for (PetscInt i = 0; i < bu->nb; i++) {
713: PetscReal n, na, nr;
714: PetscInt n_loc, na_loc, nr_loc;
716: 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));
717: if (wnormtype == NORM_INFINITY) {
718: *norm = PetscMax(*norm, n);
719: *norma = PetscMax(*norma, na);
720: *normr = PetscMax(*normr, nr);
721: } else {
722: *norm += PetscSqr(n);
723: *norma += PetscSqr(na);
724: *normr += PetscSqr(nr);
725: }
726: *norm_loc += n_loc;
727: *norma_loc += na_loc;
728: *normr_loc += nr_loc;
729: }
730: if (wnormtype == NORM_2) {
731: *norm = PetscSqrtReal(*norm);
732: *norma = PetscSqrtReal(*norma);
733: *normr = PetscSqrtReal(*normr);
734: }
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
739: {
740: PetscFunctionBegin;
741: ops->duplicate = VecDuplicate_Nest;
742: ops->duplicatevecs = VecDuplicateVecs_Default;
743: ops->destroyvecs = VecDestroyVecs_Default;
744: ops->dot = VecDot_Nest;
745: ops->mdot = VecMDot_Nest;
746: ops->norm = VecNorm_Nest;
747: ops->tdot = VecTDot_Nest;
748: ops->mtdot = VecMTDot_Nest;
749: ops->scale = VecScale_Nest;
750: ops->copy = VecCopy_Nest;
751: ops->set = VecSet_Nest;
752: ops->swap = VecSwap_Nest;
753: ops->axpy = VecAXPY_Nest;
754: ops->axpby = VecAXPBY_Nest;
755: ops->maxpy = VecMAXPY_Nest;
756: ops->aypx = VecAYPX_Nest;
757: ops->waxpy = VecWAXPY_Nest;
758: ops->axpbypcz = NULL;
759: ops->pointwisemult = VecPointwiseMult_Nest;
760: ops->pointwisedivide = VecPointwiseDivide_Nest;
761: ops->setvalues = NULL;
762: ops->assemblybegin = VecAssemblyBegin_Nest;
763: ops->assemblyend = VecAssemblyEnd_Nest;
764: ops->getarray = VecGetArray_Nest;
765: ops->getsize = VecGetSize_Nest;
766: ops->getlocalsize = VecGetLocalSize_Nest;
767: ops->restorearray = VecRestoreArray_Nest;
768: ops->restorearrayread = VecRestoreArrayRead_Nest;
769: ops->max = VecMax_Nest;
770: ops->min = VecMin_Nest;
771: ops->setrandom = NULL;
772: ops->setoption = NULL;
773: ops->setvaluesblocked = NULL;
774: ops->destroy = VecDestroy_Nest;
775: ops->view = VecView_Nest;
776: ops->placearray = NULL;
777: ops->replacearray = NULL;
778: ops->dot_local = VecDot_Nest;
779: ops->tdot_local = VecTDot_Nest;
780: ops->norm_local = VecNorm_Nest;
781: ops->mdot_local = VecMDot_Nest;
782: ops->mtdot_local = VecMTDot_Nest;
783: ops->load = NULL;
784: ops->reciprocal = VecReciprocal_Nest;
785: ops->conjugate = VecConjugate_Nest;
786: ops->setlocaltoglobalmapping = NULL;
787: ops->setvalueslocal = NULL;
788: ops->resetarray = NULL;
789: ops->setfromoptions = NULL;
790: ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
791: ops->load = NULL;
792: ops->pointwisemax = NULL;
793: ops->pointwisemaxabs = NULL;
794: ops->pointwisemin = NULL;
795: ops->getvalues = NULL;
796: ops->sqrt = NULL;
797: ops->abs = NULL;
798: ops->exp = NULL;
799: ops->shift = NULL;
800: ops->create = NULL;
801: ops->stridegather = NULL;
802: ops->stridescatter = NULL;
803: ops->dotnorm2 = VecDotNorm2_Nest;
804: ops->getsubvector = VecGetSubVector_Nest;
805: ops->restoresubvector = VecRestoreSubVector_Nest;
806: ops->axpbypcz = VecAXPBYPCZ_Nest;
807: ops->concatenate = VecConcatenate_Nest;
808: ops->createlocalvector = VecCreateLocalVector_Nest;
809: ops->getlocalvector = VecGetLocalVector_Nest;
810: ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
811: ops->restorelocalvector = VecRestoreLocalVector_Nest;
812: ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
813: ops->setrandom = VecSetRandom_Nest;
814: ops->errorwnorm = VecErrorWeightedNorms_Nest;
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
819: {
820: Vec_Nest *b = (Vec_Nest *)x->data;
821: PetscInt i;
822: PetscInt row;
824: PetscFunctionBegin;
825: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
826: for (i = 0; i < m; i++) {
827: row = idxm[i];
828: PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
829: vec[i] = b->v[row];
830: }
831: PetscFunctionReturn(PETSC_SUCCESS);
832: }
834: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
835: {
836: PetscFunctionBegin;
837: PetscCall(PetscObjectStateIncrease((PetscObject)X));
838: PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: /*@
843: VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
845: Not Collective
847: Input Parameters:
848: + X - nest vector
849: - idxm - index of the vector within the nest
851: Output Parameter:
852: . sx - vector at index `idxm` within the nest
854: Level: developer
856: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
857: @*/
858: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
859: {
860: PetscFunctionBegin;
861: PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
865: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
866: {
867: Vec_Nest *b = (Vec_Nest *)X->data;
869: PetscFunctionBegin;
870: PetscCall(PetscObjectStateIncrease((PetscObject)X));
871: if (N) *N = b->nb;
872: if (sx) *sx = b->v;
873: PetscFunctionReturn(PETSC_SUCCESS);
874: }
876: /*@C
877: VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
879: Not Collective
881: Input Parameter:
882: . X - nest vector
884: Output Parameters:
885: + N - number of nested vecs
886: - sx - array of vectors, can pass in `NULL`
888: Level: developer
890: Note:
891: The user should not free the array `sx`.
893: Fortran Notes:
894: The caller must allocate the array to hold the subvectors and pass it in.
896: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
897: @*/
898: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
899: {
900: PetscFunctionBegin;
901: PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
906: {
907: Vec_Nest *bx = (Vec_Nest *)X->data;
908: PetscInt i, offset = 0, n = 0, bs;
909: IS is;
910: PetscBool issame = PETSC_FALSE;
911: PetscInt N = 0;
913: PetscFunctionBegin;
914: /* check if idxm < bx->nb */
915: 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);
917: PetscCall(PetscObjectReference((PetscObject)x));
918: PetscCall(VecDestroy(&bx->v[idxm]));
919: bx->v[idxm] = x;
921: /* check if we need to update the IS for the block */
922: offset = X->map->rstart;
923: for (i = 0; i < idxm; i++) {
924: n = 0;
925: PetscCall(VecGetLocalSize(bx->v[i], &n));
926: offset += n;
927: }
929: /* get the local size and block size */
930: PetscCall(VecGetLocalSize(x, &n));
931: PetscCall(VecGetBlockSize(x, &bs));
933: /* create the new IS */
934: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
935: PetscCall(ISSetBlockSize(is, bs));
937: /* check if they are equal */
938: PetscCall(ISEqual(is, bx->is[idxm], &issame));
940: if (!issame) {
941: /* The IS of given vector has a different layout compared to the existing block vector.
942: Destroy the existing reference and update the IS. */
943: PetscCall(ISDestroy(&bx->is[idxm]));
944: PetscCall(ISDuplicate(is, &bx->is[idxm]));
945: PetscCall(ISCopy(is, bx->is[idxm]));
947: offset += n;
948: /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
949: for (i = idxm + 1; i < bx->nb; i++) {
950: /* get the local size and block size */
951: PetscCall(VecGetLocalSize(bx->v[i], &n));
952: PetscCall(VecGetBlockSize(bx->v[i], &bs));
954: /* destroy the old and create the new IS */
955: PetscCall(ISDestroy(&bx->is[i]));
956: PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
957: PetscCall(ISSetBlockSize(bx->is[i], bs));
959: offset += n;
960: }
962: n = 0;
963: PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
964: PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
965: PetscCall(PetscLayoutSetSize(X->map, N));
966: PetscCall(PetscLayoutSetLocalSize(X->map, n));
967: }
969: PetscCall(ISDestroy(&is));
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
974: {
975: PetscFunctionBegin;
976: PetscCall(PetscObjectStateIncrease((PetscObject)X));
977: PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: /*@
982: VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
984: Not Collective
986: Input Parameters:
987: + X - nest vector
988: . idxm - index of the vector within the nest vector
989: - sx - vector at index `idxm` within the nest vector
991: Level: developer
993: Notes:
994: The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
996: The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.
998: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
999: @*/
1000: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
1001: {
1002: PetscFunctionBegin;
1003: PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1008: {
1009: PetscInt i;
1011: PetscFunctionBegin;
1012: PetscCall(PetscObjectStateIncrease((PetscObject)X));
1013: for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: /*@
1018: VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1020: Not Collective
1022: Input Parameters:
1023: + X - nest vector
1024: . N - number of component vecs in `sx`
1025: . idxm - indices of component vectors that are to be replaced
1026: - sx - array of vectors
1028: Level: developer
1030: Notes:
1031: The components in the vector array `sx` do not have to be of the same size as corresponding
1032: components in `X`. The user can also free the array `sx` after the call.
1034: The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.
1036: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1037: @*/
1038: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
1039: {
1040: PetscFunctionBegin;
1041: PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1046: {
1047: Vec_Nest *b = (Vec_Nest *)X->data;
1049: PetscFunctionBegin;
1050: *N = b->nb;
1051: PetscFunctionReturn(PETSC_SUCCESS);
1052: }
1054: /*@
1055: VecNestGetSize - Returns the size of the nest vector.
1057: Not Collective
1059: Input Parameter:
1060: . X - nest vector
1062: Output Parameter:
1063: . N - number of nested vecs
1065: Level: developer
1067: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1068: @*/
1069: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1070: {
1071: PetscFunctionBegin;
1073: PetscAssertPointer(N, 2);
1074: PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1075: PetscFunctionReturn(PETSC_SUCCESS);
1076: }
1078: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1079: {
1080: Vec_Nest *ctx = (Vec_Nest *)V->data;
1081: PetscInt i;
1083: PetscFunctionBegin;
1084: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1086: ctx->nb = nb;
1087: PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1089: /* Create space */
1090: PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1091: for (i = 0; i < ctx->nb; i++) {
1092: ctx->v[i] = x[i];
1093: PetscCall(PetscObjectReference((PetscObject)x[i]));
1094: /* Do not allocate memory for internal sub blocks */
1095: }
1097: PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1099: ctx->setup_called = PETSC_TRUE;
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1104: {
1105: Vec_Nest *ctx = (Vec_Nest *)V->data;
1106: PetscInt i, offset, m, n, M, N;
1108: PetscFunctionBegin;
1109: (void)nb;
1110: if (is) { /* Do some consistency checks and reference the is */
1111: offset = V->map->rstart;
1112: for (i = 0; i < ctx->nb; i++) {
1113: PetscCall(ISGetSize(is[i], &M));
1114: PetscCall(VecGetSize(ctx->v[i], &N));
1115: 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);
1116: PetscCall(ISGetLocalSize(is[i], &m));
1117: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1118: 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);
1119: PetscCall(PetscObjectReference((PetscObject)is[i]));
1120: ctx->is[i] = is[i];
1121: offset += n;
1122: }
1123: } else { /* Create a contiguous ISStride for each entry */
1124: offset = V->map->rstart;
1125: for (i = 0; i < ctx->nb; i++) {
1126: PetscInt bs;
1127: PetscCall(VecGetLocalSize(ctx->v[i], &n));
1128: PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1129: PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1130: PetscCall(ISSetBlockSize(ctx->is[i], bs));
1131: offset += n;
1132: }
1133: }
1134: PetscFunctionReturn(PETSC_SUCCESS);
1135: }
1137: /*MC
1138: VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1140: Level: intermediate
1142: Notes:
1143: This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1144: It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1146: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1147: M*/
1149: /*@C
1150: VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1152: Collective
1154: Input Parameters:
1155: + comm - Communicator for the new `Vec`
1156: . nb - number of nested blocks
1157: . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1158: - x - array of `nb` sub-vectors
1160: Output Parameter:
1161: . Y - new vector
1163: Level: advanced
1165: .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1166: @*/
1167: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1168: {
1169: Vec V;
1170: Vec_Nest *s;
1171: PetscInt n, N;
1173: PetscFunctionBegin;
1174: PetscCall(VecCreate(comm, &V));
1175: PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1177: /* allocate and set pointer for implementation data */
1178: PetscCall(PetscNew(&s));
1179: V->data = (void *)s;
1180: s->setup_called = PETSC_FALSE;
1181: s->nb = -1;
1182: s->v = NULL;
1184: PetscCall(VecSetUp_Nest_Private(V, nb, x));
1186: n = N = 0;
1187: PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1188: PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1189: PetscCall(PetscLayoutSetSize(V->map, N));
1190: PetscCall(PetscLayoutSetLocalSize(V->map, n));
1191: PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1192: PetscCall(PetscLayoutSetUp(V->map));
1194: PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1196: PetscCall(VecNestSetOps_Private(V->ops));
1197: V->petscnative = PETSC_FALSE;
1199: /* expose block api's */
1200: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1201: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1202: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1203: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1204: PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1206: *Y = V;
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }