1 #include <../src/vec/vec/impls/nest/vecnestimpl.h> /*I "petscvec.h" I*/
2
3 /* check all blocks are filled */
VecAssemblyBegin_Nest(Vec v)4 static PetscErrorCode VecAssemblyBegin_Nest(Vec v)
5 {
6 Vec_Nest *vs = (Vec_Nest *)v->data;
7 PetscInt i;
8
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 }
16
VecAssemblyEnd_Nest(Vec v)17 static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
18 {
19 Vec_Nest *vs = (Vec_Nest *)v->data;
20 PetscInt i;
21
22 PetscFunctionBegin;
23 for (i = 0; i < vs->nb; i++) PetscCall(VecAssemblyEnd(vs->v[i]));
24 PetscFunctionReturn(PETSC_SUCCESS);
25 }
26
VecDestroy_Nest(Vec v)27 static PetscErrorCode VecDestroy_Nest(Vec v)
28 {
29 Vec_Nest *vs = (Vec_Nest *)v->data;
30 PetscInt i;
31
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));
44
45 PetscCall(PetscFree(vs));
46 PetscFunctionReturn(PETSC_SUCCESS);
47 }
48
VecCopy_Nest(Vec x,Vec y)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;
54
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 }
60
VecDuplicate_Nest(Vec x,Vec * y)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;
67
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 }
77
VecDot_Nest(Vec x,Vec y,PetscScalar * val)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;
84
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 }
96
VecTDot_Nest(Vec x,Vec y,PetscScalar * val)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;
103
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 }
115
VecDotNorm2_Nest(Vec x,Vec y,PetscScalar * dp,PetscScalar * nm)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;
123
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 }
138
VecAXPY_Nest(Vec y,PetscScalar alpha,Vec x)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;
144
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 }
151
VecAYPX_Nest(Vec y,PetscScalar alpha,Vec x)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;
157
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 }
164
VecAXPBY_Nest(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)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;
170
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 }
177
VecAXPBYPCZ_Nest(Vec z,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec x,Vec y)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;
184
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 }
191
VecScale_Nest(Vec x,PetscScalar alpha)192 static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
193 {
194 Vec_Nest *bx = (Vec_Nest *)x->data;
195 PetscInt i, nr;
196
197 PetscFunctionBegin;
198 nr = bx->nb;
199 for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
VecPointwiseMult_Nest(Vec w,Vec x,Vec y)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;
209
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 }
216
VecPointwiseDivide_Nest(Vec w,Vec x,Vec y)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;
223
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 }
230
VecReciprocal_Nest(Vec x)231 static PetscErrorCode VecReciprocal_Nest(Vec x)
232 {
233 Vec_Nest *bx = (Vec_Nest *)x->data;
234 PetscInt i, nr;
235
236 PetscFunctionBegin;
237 nr = bx->nb;
238 for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
VecNorm_Nest(Vec xin,NormType type,PetscReal * z)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;
248
249 PetscFunctionBegin;
250 nr = bx->nb;
251 _z = 0.0;
252
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 }
268
269 *z = _z;
270 PetscFunctionReturn(PETSC_SUCCESS);
271 }
272
VecMAXPY_Nest(Vec y,PetscInt nv,const PetscScalar alpha[],Vec * x)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 }
280
VecMDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar * val)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 }
288
VecMTDot_Nest(Vec x,PetscInt nv,const Vec y[],PetscScalar * val)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 }
296
VecSet_Nest(Vec x,PetscScalar alpha)297 static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
298 {
299 Vec_Nest *bx = (Vec_Nest *)x->data;
300 PetscInt i, nr;
301
302 PetscFunctionBegin;
303 nr = bx->nb;
304 for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
305 PetscFunctionReturn(PETSC_SUCCESS);
306 }
307
VecConjugate_Nest(Vec x)308 static PetscErrorCode VecConjugate_Nest(Vec x)
309 {
310 Vec_Nest *bx = (Vec_Nest *)x->data;
311 PetscInt j, nr;
312
313 PetscFunctionBegin;
314 nr = bx->nb;
315 for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
316 PetscFunctionReturn(PETSC_SUCCESS);
317 }
318
VecSwap_Nest(Vec x,Vec y)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;
324
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 }
331
VecWAXPY_Nest(Vec w,PetscScalar alpha,Vec x,Vec y)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;
338
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 }
345
VecMin_Nest(Vec x,PetscInt * p,PetscReal * v)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;
351
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;
366
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 }
378
VecMax_Nest(Vec x,PetscInt * p,PetscReal * v)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;
384
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;
399
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 }
411
VecView_Nest(Vec x,PetscViewer viewer)412 static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
413 {
414 Vec_Nest *bx = (Vec_Nest *)x->data;
415 PetscBool isascii;
416 PetscInt i;
417
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;
427
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));
432
433 PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT "\n", i, name, prefix, type, NR));
434
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 }
443
VecSize_Nest_Recursive(Vec x,PetscBool globalsize,PetscInt * L)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;
449
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 }
459
460 /* Otherwise we have a nest */
461 bx = (Vec_Nest *)x->data;
462 nr = bx->nb;
463
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 }
468
469 /* Returns the sum of the global size of all the constituent vectors in the nest */
VecGetSize_Nest(Vec x,PetscInt * N)470 static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
471 {
472 PetscFunctionBegin;
473 *N = x->map->N;
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
477 /* Returns the sum of the local size of all the constituent vectors in the nest */
VecGetLocalSize_Nest(Vec x,PetscInt * n)478 static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
479 {
480 PetscFunctionBegin;
481 *n = x->map->n;
482 PetscFunctionReturn(PETSC_SUCCESS);
483 }
484
VecMaxPointwiseDivide_Nest(Vec x,Vec y,PetscReal * max)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;
491
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 }
503
VecGetSubVector_Nest(Vec X,IS is,Vec * x)504 static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
505 {
506 Vec_Nest *bx = (Vec_Nest *)X->data;
507 PetscInt i;
508
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 }
523
VecRestoreSubVector_Nest(Vec X,IS is,Vec * x)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 }
531
VecGetArray_Nest(Vec X,PetscScalar ** x)532 static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
533 {
534 Vec_Nest *bx = (Vec_Nest *)X->data;
535 PetscInt i, m, rstart, rend;
536
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 }
560
VecRestoreArray_Nest(Vec X,PetscScalar ** x)561 static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
562 {
563 Vec_Nest *bx = (Vec_Nest *)X->data;
564 PetscInt i, m, rstart, rend;
565
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 }
590
VecRestoreArrayRead_Nest(Vec X,const PetscScalar ** x)591 static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
592 {
593 PetscFunctionBegin;
594 PetscCall(PetscFree(*x));
595 PetscFunctionReturn(PETSC_SUCCESS);
596 }
597
VecConcatenate_Nest(PetscInt nx,const Vec X[],Vec * Y,IS * x_is[])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 }
604
VecCreateLocalVector_Nest(Vec v,Vec * w)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;
612
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;
623
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 }
636
VecGetLocalVector_Nest(Vec v,Vec w)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;
642
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 }
649
VecRestoreLocalVector_Nest(Vec v,Vec w)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;
655
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 }
663
VecGetLocalVectorRead_Nest(Vec v,Vec w)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;
669
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 }
676
VecRestoreLocalVectorRead_Nest(Vec v,Vec w)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;
682
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 }
689
VecSetRandom_Nest(Vec v,PetscRandom r)690 static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
691 {
692 Vec_Nest *bv = (Vec_Nest *)v->data;
693
694 PetscFunctionBegin;
695 for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
696 PetscFunctionReturn(PETSC_SUCCESS);
697 }
698
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)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;
706
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;
721
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 }
743
VecNestSetOps_Private(struct _VecOps * ops)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->resetarray = NULL;
794 ops->setfromoptions = NULL;
795 ops->maxpointwisedivide = VecMaxPointwiseDivide_Nest;
796 ops->load = NULL;
797 ops->pointwisemax = NULL;
798 ops->pointwisemaxabs = NULL;
799 ops->pointwisemin = NULL;
800 ops->getvalues = NULL;
801 ops->sqrt = NULL;
802 ops->abs = NULL;
803 ops->exp = NULL;
804 ops->shift = NULL;
805 ops->create = NULL;
806 ops->stridegather = NULL;
807 ops->stridescatter = NULL;
808 ops->dotnorm2 = VecDotNorm2_Nest;
809 ops->getsubvector = VecGetSubVector_Nest;
810 ops->restoresubvector = VecRestoreSubVector_Nest;
811 ops->axpbypcz = VecAXPBYPCZ_Nest;
812 ops->concatenate = VecConcatenate_Nest;
813 ops->createlocalvector = VecCreateLocalVector_Nest;
814 ops->getlocalvector = VecGetLocalVector_Nest;
815 ops->getlocalvectorread = VecGetLocalVectorRead_Nest;
816 ops->restorelocalvector = VecRestoreLocalVector_Nest;
817 ops->restorelocalvectorread = VecRestoreLocalVectorRead_Nest;
818 ops->setrandom = VecSetRandom_Nest;
819 ops->errorwnorm = VecErrorWeightedNorms_Nest;
820 PetscFunctionReturn(PETSC_SUCCESS);
821 }
822
VecNestGetSubVecs_Private(Vec x,PetscInt m,const PetscInt idxm[],Vec vec[])823 static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
824 {
825 Vec_Nest *b = (Vec_Nest *)x->data;
826 PetscInt i;
827 PetscInt row;
828
829 PetscFunctionBegin;
830 if (!m) PetscFunctionReturn(PETSC_SUCCESS);
831 for (i = 0; i < m; i++) {
832 row = idxm[i];
833 PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
834 vec[i] = b->v[row];
835 }
836 PetscFunctionReturn(PETSC_SUCCESS);
837 }
838
VecNestGetSubVec_Nest(Vec X,PetscInt idxm,Vec * sx)839 static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
840 {
841 PetscFunctionBegin;
842 PetscCall(PetscObjectStateIncrease((PetscObject)X));
843 PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
844 PetscFunctionReturn(PETSC_SUCCESS);
845 }
846
847 /*@
848 VecNestGetSubVec - Returns a single, sub-vector from a nest vector.
849
850 Not Collective
851
852 Input Parameters:
853 + X - nest vector
854 - idxm - index of the vector within the nest
855
856 Output Parameter:
857 . sx - vector at index `idxm` within the nest
858
859 Level: developer
860
861 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
862 @*/
VecNestGetSubVec(Vec X,PetscInt idxm,Vec * sx)863 PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
864 {
865 PetscFunctionBegin;
866 PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
867 PetscFunctionReturn(PETSC_SUCCESS);
868 }
869
VecNestGetSubVecs_Nest(Vec X,PetscInt * N,Vec ** sx)870 static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
871 {
872 Vec_Nest *b = (Vec_Nest *)X->data;
873
874 PetscFunctionBegin;
875 PetscCall(PetscObjectStateIncrease((PetscObject)X));
876 if (N) *N = b->nb;
877 if (sx) *sx = b->v;
878 PetscFunctionReturn(PETSC_SUCCESS);
879 }
880
881 /*@C
882 VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.
883
884 Not Collective
885
886 Input Parameter:
887 . X - nest vector
888
889 Output Parameters:
890 + N - number of nested vecs
891 - sx - array of vectors, can pass in `NULL`
892
893 Level: developer
894
895 Note:
896 The user should not free the array `sx`.
897
898 Fortran Notes:
899 The caller must allocate the array to hold the subvectors and pass it in.
900
901 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
902 @*/
VecNestGetSubVecs(Vec X,PetscInt * N,Vec * sx[])903 PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
904 {
905 PetscFunctionBegin;
906 PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
907 PetscFunctionReturn(PETSC_SUCCESS);
908 }
909
VecNestSetSubVec_Private(Vec X,PetscInt idxm,Vec x)910 static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
911 {
912 Vec_Nest *bx = (Vec_Nest *)X->data;
913 PetscInt i, offset = 0, n = 0, bs;
914 IS is;
915 PetscBool issame = PETSC_FALSE;
916 PetscInt N = 0;
917
918 PetscFunctionBegin;
919 /* check if idxm < bx->nb */
920 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);
921
922 PetscCall(PetscObjectReference((PetscObject)x));
923 PetscCall(VecDestroy(&bx->v[idxm]));
924 bx->v[idxm] = x;
925
926 /* check if we need to update the IS for the block */
927 offset = X->map->rstart;
928 for (i = 0; i < idxm; i++) {
929 n = 0;
930 PetscCall(VecGetLocalSize(bx->v[i], &n));
931 offset += n;
932 }
933
934 /* get the local size and block size */
935 PetscCall(VecGetLocalSize(x, &n));
936 PetscCall(VecGetBlockSize(x, &bs));
937
938 /* create the new IS */
939 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
940 PetscCall(ISSetBlockSize(is, bs));
941
942 /* check if they are equal */
943 PetscCall(ISEqual(is, bx->is[idxm], &issame));
944
945 if (!issame) {
946 /* The IS of given vector has a different layout compared to the existing block vector.
947 Destroy the existing reference and update the IS. */
948 PetscCall(ISDestroy(&bx->is[idxm]));
949 PetscCall(ISDuplicate(is, &bx->is[idxm]));
950 PetscCall(ISCopy(is, bx->is[idxm]));
951
952 offset += n;
953 /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
954 for (i = idxm + 1; i < bx->nb; i++) {
955 /* get the local size and block size */
956 PetscCall(VecGetLocalSize(bx->v[i], &n));
957 PetscCall(VecGetBlockSize(bx->v[i], &bs));
958
959 /* destroy the old and create the new IS */
960 PetscCall(ISDestroy(&bx->is[i]));
961 PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
962 PetscCall(ISSetBlockSize(bx->is[i], bs));
963
964 offset += n;
965 }
966
967 n = 0;
968 PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
969 PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
970 PetscCall(PetscLayoutSetSize(X->map, N));
971 PetscCall(PetscLayoutSetLocalSize(X->map, n));
972 }
973
974 PetscCall(ISDestroy(&is));
975 PetscFunctionReturn(PETSC_SUCCESS);
976 }
977
VecNestSetSubVec_Nest(Vec X,PetscInt idxm,Vec sx)978 static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
979 {
980 PetscFunctionBegin;
981 PetscCall(PetscObjectStateIncrease((PetscObject)X));
982 PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
983 PetscFunctionReturn(PETSC_SUCCESS);
984 }
985
986 /*@
987 VecNestSetSubVec - Set a single component vector in a nest vector at specified index.
988
989 Not Collective
990
991 Input Parameters:
992 + X - nest vector
993 . idxm - index of the vector within the nest vector
994 - sx - vector at index `idxm` within the nest vector
995
996 Level: developer
997
998 Notes:
999 The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.
1000
1001 The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.
1002
1003 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
1004 @*/
VecNestSetSubVec(Vec X,PetscInt idxm,Vec sx)1005 PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
1006 {
1007 PetscFunctionBegin;
1008 PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
1009 PetscFunctionReturn(PETSC_SUCCESS);
1010 }
1011
VecNestSetSubVecs_Nest(Vec X,PetscInt N,PetscInt * idxm,Vec * sx)1012 static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
1013 {
1014 PetscInt i;
1015
1016 PetscFunctionBegin;
1017 PetscCall(PetscObjectStateIncrease((PetscObject)X));
1018 for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
1019 PetscFunctionReturn(PETSC_SUCCESS);
1020 }
1021
1022 /*@
1023 VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.
1024
1025 Not Collective
1026
1027 Input Parameters:
1028 + X - nest vector
1029 . N - number of component vecs in `sx`
1030 . idxm - indices of component vectors that are to be replaced
1031 - sx - array of vectors
1032
1033 Level: developer
1034
1035 Notes:
1036 The components in the vector array `sx` do not have to be of the same size as corresponding
1037 components in `X`. The user can also free the array `sx` after the call.
1038
1039 The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.
1040
1041 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
1042 @*/
VecNestSetSubVecs(Vec X,PetscInt N,PetscInt idxm[],Vec sx[])1043 PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
1044 {
1045 PetscFunctionBegin;
1046 PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
1047 PetscFunctionReturn(PETSC_SUCCESS);
1048 }
1049
VecNestGetSize_Nest(Vec X,PetscInt * N)1050 static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
1051 {
1052 Vec_Nest *b = (Vec_Nest *)X->data;
1053
1054 PetscFunctionBegin;
1055 *N = b->nb;
1056 PetscFunctionReturn(PETSC_SUCCESS);
1057 }
1058
1059 /*@
1060 VecNestGetSize - Returns the size of the nest vector.
1061
1062 Not Collective
1063
1064 Input Parameter:
1065 . X - nest vector
1066
1067 Output Parameter:
1068 . N - number of nested vecs
1069
1070 Level: developer
1071
1072 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
1073 @*/
VecNestGetSize(Vec X,PetscInt * N)1074 PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
1075 {
1076 PetscFunctionBegin;
1077 PetscValidHeaderSpecific(X, VEC_CLASSID, 1);
1078 PetscAssertPointer(N, 2);
1079 PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
1080 PetscFunctionReturn(PETSC_SUCCESS);
1081 }
1082
VecSetUp_Nest_Private(Vec V,PetscInt nb,Vec x[])1083 static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
1084 {
1085 Vec_Nest *ctx = (Vec_Nest *)V->data;
1086 PetscInt i;
1087
1088 PetscFunctionBegin;
1089 if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
1090
1091 ctx->nb = nb;
1092 PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");
1093
1094 /* Create space */
1095 PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
1096 for (i = 0; i < ctx->nb; i++) {
1097 ctx->v[i] = x[i];
1098 PetscCall(PetscObjectReference((PetscObject)x[i]));
1099 /* Do not allocate memory for internal sub blocks */
1100 }
1101
1102 PetscCall(PetscMalloc1(ctx->nb, &ctx->is));
1103
1104 ctx->setup_called = PETSC_TRUE;
1105 PetscFunctionReturn(PETSC_SUCCESS);
1106 }
1107
VecSetUp_NestIS_Private(Vec V,PetscInt nb,IS is[])1108 static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
1109 {
1110 Vec_Nest *ctx = (Vec_Nest *)V->data;
1111 PetscInt i, offset, m, n, M, N;
1112
1113 PetscFunctionBegin;
1114 (void)nb;
1115 if (is) { /* Do some consistency checks and reference the is */
1116 offset = V->map->rstart;
1117 for (i = 0; i < ctx->nb; i++) {
1118 PetscCall(ISGetSize(is[i], &M));
1119 PetscCall(VecGetSize(ctx->v[i], &N));
1120 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);
1121 PetscCall(ISGetLocalSize(is[i], &m));
1122 PetscCall(VecGetLocalSize(ctx->v[i], &n));
1123 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);
1124 PetscCall(PetscObjectReference((PetscObject)is[i]));
1125 ctx->is[i] = is[i];
1126 offset += n;
1127 }
1128 } else { /* Create a contiguous ISStride for each entry */
1129 offset = V->map->rstart;
1130 for (i = 0; i < ctx->nb; i++) {
1131 PetscInt bs;
1132 PetscCall(VecGetLocalSize(ctx->v[i], &n));
1133 PetscCall(VecGetBlockSize(ctx->v[i], &bs));
1134 PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
1135 PetscCall(ISSetBlockSize(ctx->is[i], bs));
1136 offset += n;
1137 }
1138 }
1139 PetscFunctionReturn(PETSC_SUCCESS);
1140 }
1141
1142 /*MC
1143 VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.
1144
1145 Level: intermediate
1146
1147 Notes:
1148 This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
1149 It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.
1150
1151 .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
1152 M*/
1153
1154 /*@C
1155 VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately
1156
1157 Collective
1158
1159 Input Parameters:
1160 + comm - Communicator for the new `Vec`
1161 . nb - number of nested blocks
1162 . is - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
1163 - x - array of `nb` sub-vectors
1164
1165 Output Parameter:
1166 . Y - new vector
1167
1168 Level: advanced
1169
1170 .seealso: `VECNEST`, [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
1171 @*/
VecCreateNest(MPI_Comm comm,PetscInt nb,IS is[],Vec x[],Vec * Y)1172 PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
1173 {
1174 Vec V;
1175 Vec_Nest *s;
1176 PetscInt n, N;
1177
1178 PetscFunctionBegin;
1179 PetscCall(VecCreate(comm, &V));
1180 PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));
1181
1182 /* allocate and set pointer for implementation data */
1183 PetscCall(PetscNew(&s));
1184 V->data = (void *)s;
1185 s->setup_called = PETSC_FALSE;
1186 s->nb = -1;
1187 s->v = NULL;
1188
1189 PetscCall(VecSetUp_Nest_Private(V, nb, x));
1190
1191 n = N = 0;
1192 PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
1193 PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
1194 PetscCall(PetscLayoutSetSize(V->map, N));
1195 PetscCall(PetscLayoutSetLocalSize(V->map, n));
1196 PetscCall(PetscLayoutSetBlockSize(V->map, 1));
1197 PetscCall(PetscLayoutSetUp(V->map));
1198
1199 PetscCall(VecSetUp_NestIS_Private(V, nb, is));
1200
1201 PetscCall(VecNestSetOps_Private(V->ops));
1202 V->petscnative = PETSC_FALSE;
1203
1204 /* expose block api's */
1205 PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
1206 PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
1207 PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
1208 PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
1209 PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));
1210
1211 *Y = V;
1212 PetscFunctionReturn(PETSC_SUCCESS);
1213 }
1214