Actual source code: ex37.c
1: static char help[] = "Nest vector functionality.\n\n";
3: #include <petscvec.h>
5: static PetscErrorCode GetISs(Vec vecs[], IS is[], PetscBool inv)
6: {
7: PetscInt rstart[2], rend[2];
9: PetscFunctionBegin;
10: PetscCall(VecGetOwnershipRange(vecs[0], &rstart[0], &rend[0]));
11: PetscCall(VecGetOwnershipRange(vecs[1], &rstart[1], &rend[1]));
12: if (!inv) {
13: PetscCall(ISCreateStride(PETSC_COMM_WORLD, rend[0] - rstart[0], rstart[0] + rstart[1], 1, &is[0]));
14: PetscCall(ISCreateStride(PETSC_COMM_WORLD, rend[1] - rstart[1], rend[0] + rstart[1], 1, &is[1]));
15: } else {
16: PetscCall(ISCreateStride(PETSC_COMM_WORLD, rend[0] - rstart[0], rend[0] + rend[1] - 1, -1, &is[0]));
17: PetscCall(ISCreateStride(PETSC_COMM_WORLD, rend[1] - rstart[1], rstart[0] + rend[1] - 1, -1, &is[1]));
18: }
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: PetscErrorCode convert_from_nest(Vec X, Vec *Y)
23: {
24: const PetscScalar *v;
25: PetscInt rstart, n, N;
27: PetscFunctionBegin;
28: PetscCall(VecGetLocalSize(X, &n));
29: PetscCall(VecGetSize(X, &N));
30: PetscCall(VecGetOwnershipRange(X, &rstart, NULL));
31: PetscCall(VecCreate(PetscObjectComm((PetscObject)X), Y));
32: PetscCall(VecSetSizes(*Y, n, N));
33: PetscCall(VecSetType(*Y, VECSTANDARD)); // We always use a CPU only version
34: PetscCall(VecGetArrayRead(X, &v));
35: for (PetscInt i = 0; i < n; i++) PetscCall(VecSetValue(*Y, rstart + i, v[i], INSERT_VALUES));
36: PetscCall(VecRestoreArrayRead(X, &v));
37: PetscCall(VecAssemblyBegin(*Y));
38: PetscCall(VecAssemblyEnd(*Y));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: PetscErrorCode test_view(void)
43: {
44: Vec X, lX, a, b;
45: Vec c, d, e, f;
46: Vec tmp_buf[2];
47: IS tmp_is[2];
48: PetscInt index, n;
49: PetscReal val;
50: PetscInt list[] = {0, 1, 2};
51: PetscScalar vals[] = {0.5, 0.25, 0.125};
52: PetscScalar *x, *lx;
53: PetscBool explcit = PETSC_FALSE, inv = PETSC_FALSE;
55: PetscFunctionBegin;
56: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
58: PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
59: PetscCall(VecSetSizes(c, PETSC_DECIDE, 3));
60: PetscCall(VecSetFromOptions(c));
61: PetscCall(VecDuplicate(c, &d));
62: PetscCall(VecDuplicate(c, &e));
63: PetscCall(VecDuplicate(c, &f));
65: PetscCall(VecSet(c, 1.0));
66: PetscCall(VecSet(d, 2.0));
67: PetscCall(VecSet(e, 3.0));
68: PetscCall(VecSetValues(f, 3, list, vals, INSERT_VALUES));
69: PetscCall(VecAssemblyBegin(f));
70: PetscCall(VecAssemblyEnd(f));
71: PetscCall(VecScale(f, 10.0));
73: tmp_buf[0] = e;
74: tmp_buf[1] = f;
75: PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicit_is", &explcit, 0));
76: PetscCall(GetISs(tmp_buf, tmp_is, PETSC_FALSE));
77: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, explcit ? tmp_is : NULL, tmp_buf, &b));
78: PetscCall(VecDestroy(&e));
79: PetscCall(VecDestroy(&f));
80: PetscCall(ISDestroy(&tmp_is[0]));
81: PetscCall(ISDestroy(&tmp_is[1]));
83: tmp_buf[0] = c;
84: tmp_buf[1] = d;
85: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &a));
86: PetscCall(VecDestroy(&c));
87: PetscCall(VecDestroy(&d));
89: PetscCall(PetscOptionsGetBool(NULL, NULL, "-inv", &inv, 0));
90: tmp_buf[0] = a;
91: tmp_buf[1] = b;
92: if (inv) {
93: PetscCall(GetISs(tmp_buf, tmp_is, inv));
94: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, tmp_is, tmp_buf, &X));
95: PetscCall(ISDestroy(&tmp_is[0]));
96: PetscCall(ISDestroy(&tmp_is[1]));
97: } else {
98: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
99: }
100: PetscCall(VecDestroy(&a));
102: PetscCall(VecAssemblyBegin(X));
103: PetscCall(VecAssemblyEnd(X));
105: PetscCall(VecMax(b, &index, &val));
106: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
108: PetscCall(VecMin(b, &index, &val));
109: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
111: PetscCall(VecDestroy(&b));
113: PetscCall(VecMax(X, &index, &val));
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
115: PetscCall(VecMin(X, &index, &val));
116: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
118: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
120: Vec X2, A, R, E, vX, vX2, vA, vR, vE;
121: PetscCall(convert_from_nest(X, &vX));
122: PetscCall(VecDuplicate(X, &X2));
123: PetscCall(VecDuplicate(X, &A));
124: PetscCall(VecDuplicate(X, &R));
125: PetscCall(VecDuplicate(X, &E));
126: PetscCall(VecSetRandom(A, NULL));
127: PetscCall(VecSetRandom(R, NULL));
128: PetscCall(VecSetRandom(E, NULL));
129: PetscCall(convert_from_nest(A, &vA));
130: PetscCall(convert_from_nest(R, &vR));
131: PetscCall(convert_from_nest(E, &vE));
132: PetscCall(VecCopy(X, X2));
133: PetscCall(VecDuplicate(vX, &vX2));
134: PetscCall(VecCopy(vX, vX2));
135: PetscCall(VecScale(X2, 2.0));
136: PetscCall(VecScale(vX2, 2.0));
137: for (int nt = 0; nt < 2; nt++) {
138: NormType norm = nt ? NORM_INFINITY : NORM_2;
139: for (int e = 0; e < 2; e++) {
140: for (int a = 0; a < 2; a++) {
141: for (int r = 0; r < 2; r++) {
142: PetscReal vn, vna, vnr, nn, nna, nnr;
143: PetscInt vn_loc, vna_loc, vnr_loc, nn_loc, nna_loc, nnr_loc;
145: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing Wnorms %s: E? %d A? %d R? %d\n", norm == NORM_2 ? "2" : "inf", e, a, r));
146: PetscCall(VecErrorWeightedNorms(vX, vX2, e ? vE : NULL, norm, 0.5, a ? vA : NULL, 0.5, r ? vR : NULL, 0.0, &vn, &vn_loc, &vna, &vna_loc, &vnr, &vnr_loc));
147: PetscCall(VecErrorWeightedNorms(X, X2, e ? E : NULL, norm, 0.5, a ? A : NULL, 0.5, r ? R : NULL, 0.0, &nn, &nn_loc, &nna, &nna_loc, &nnr, &nnr_loc));
148: if (vn_loc != nn_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with total norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vn_loc, nn_loc));
149: if (vna_loc != nna_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with absolute norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vna_loc, nna_loc));
150: if (vnr_loc != nnr_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with relative norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vnr_loc, nnr_loc));
151: if (!PetscIsCloseAtTol(vna, nna, 0, PETSC_SQRT_MACHINE_EPSILON)) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with absolute norm: %1.16e %1.16e diff %1.16e\n", (double)vna, (double)nna, (double)(vna - nna)));
152: if (!PetscIsCloseAtTol(vnr, nnr, 0, PETSC_SQRT_MACHINE_EPSILON)) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with relative norm: %1.16e %1.16e diff %1.16e\n", (double)vnr, (double)nnr, (double)(vnr - nnr)));
153: }
154: }
155: }
156: }
157: PetscCall(VecDestroy(&X2));
158: PetscCall(VecDestroy(&A));
159: PetscCall(VecDestroy(&R));
160: PetscCall(VecDestroy(&E));
161: PetscCall(VecDestroy(&vX));
162: PetscCall(VecDestroy(&vX2));
163: PetscCall(VecDestroy(&vA));
164: PetscCall(VecDestroy(&vR));
165: PetscCall(VecDestroy(&vE));
167: PetscCall(VecCreateLocalVector(X, &lX));
168: PetscCall(VecGetLocalVectorRead(X, lX));
169: PetscCall(VecGetLocalSize(lX, &n));
170: PetscCall(VecGetArrayRead(lX, (const PetscScalar **)&lx));
171: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
172: for (PetscInt i = 0; i < n; i++) PetscCheck(lx[i] == x[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid data");
173: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
174: PetscCall(VecRestoreArrayRead(lX, (const PetscScalar **)&lx));
175: PetscCall(VecRestoreLocalVectorRead(X, lX));
177: PetscCall(VecDestroy(&lX));
178: PetscCall(VecDestroy(&X));
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: #if 0
183: PetscErrorCode test_vec_ops(void)
184: {
185: Vec X, a,b;
186: Vec c,d,e,f;
187: PetscScalar val;
189: PetscFunctionBegin;
190: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME));
192: PetscCall(VecCreate(PETSC_COMM_WORLD, &X));
193: PetscCall(VecSetSizes(X, 2, 2));
194: PetscCall(VecSetType(X, VECNEST));
196: PetscCall(VecCreate(PETSC_COMM_WORLD, &a));
197: PetscCall(VecSetSizes(a, 2, 2));
198: PetscCall(VecSetType(a, VECNEST));
200: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
201: PetscCall(VecSetSizes(b, 2, 2));
202: PetscCall(VecSetType(b, VECNEST));
204: /* assemble X */
205: PetscCall(VecNestSetSubVec(X, 0, a));
206: PetscCall(VecNestSetSubVec(X, 1, b));
207: PetscCall(VecAssemblyBegin(X));
208: PetscCall(VecAssemblyEnd(X));
210: PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
211: PetscCall(VecSetSizes(c, 3, 3));
212: PetscCall(VecSetType(c, VECSEQ));
213: PetscCall(VecDuplicate(c, &d));
214: PetscCall(VecDuplicate(c, &e));
215: PetscCall(VecDuplicate(c, &f));
217: PetscCall(VecSet(c, 1.0));
218: PetscCall(VecSet(d, 2.0));
219: PetscCall(VecSet(e, 3.0));
220: PetscCall(VecSet(f, 4.0));
222: /* assemble a */
223: PetscCall(VecNestSetSubVec(a, 0, c));
224: PetscCall(VecNestSetSubVec(a, 1, d));
225: PetscCall(VecAssemblyBegin(a));
226: PetscCall(VecAssemblyEnd(a));
228: /* assemble b */
229: PetscCall(VecNestSetSubVec(b, 0, e));
230: PetscCall(VecNestSetSubVec(b, 1, f));
231: PetscCall(VecAssemblyBegin(b));
232: PetscCall(VecAssemblyEnd(b));
234: PetscCall(VecDot(X,X, &val));
235: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
238: #endif
240: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
241: {
242: PetscMPIInt size;
243: Vec v;
244: PetscInt i;
245: PetscScalar vx;
247: PetscFunctionBegin;
248: PetscCallMPI(MPI_Comm_size(comm, &size));
249: PetscCall(VecCreate(comm, &v));
250: PetscCall(VecSetSizes(v, PETSC_DECIDE, length));
251: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
252: else PetscCall(VecSetType(v, VECMPI));
254: for (i = 0; i < length; i++) {
255: vx = (PetscScalar)(start_value + i * stride);
256: PetscCall(VecSetValue(v, i, vx, INSERT_VALUES));
257: }
258: PetscCall(VecAssemblyBegin(v));
259: PetscCall(VecAssemblyEnd(v));
261: *_v = v;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*
266: X = ([0,1,2,3], [10,12,14,16,18])
267: Y = ([4,7,10,13], [5,6,7,8,9])
269: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
270: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
272: */
273: PetscErrorCode test_axpy_dot_max(void)
274: {
275: Vec x1, y1, x2, y2;
276: Vec tmp_buf[2];
277: Vec X, Y;
278: PetscReal real, real2;
279: PetscScalar scalar;
280: PetscInt index;
282: PetscFunctionBegin;
283: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
285: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1));
286: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2));
288: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1));
289: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2));
291: tmp_buf[0] = x1;
292: tmp_buf[1] = x2;
293: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
294: PetscCall(VecAssemblyBegin(X));
295: PetscCall(VecAssemblyEnd(X));
296: PetscCall(VecDestroy(&x1));
297: PetscCall(VecDestroy(&x2));
299: tmp_buf[0] = y1;
300: tmp_buf[1] = y2;
301: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &Y));
302: PetscCall(VecAssemblyBegin(Y));
303: PetscCall(VecAssemblyEnd(Y));
304: PetscCall(VecDestroy(&y1));
305: PetscCall(VecDestroy(&y2));
307: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n"));
308: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
309: PetscCall(VecNestGetSubVec(Y, 0, &y1));
310: PetscCall(VecNestGetSubVec(Y, 1, &y2));
311: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n"));
312: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
313: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n"));
314: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
315: PetscCall(VecDot(X, Y, &scalar));
317: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
319: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
320: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
322: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
323: PetscCall(VecNestGetSubVec(Y, 0, &y1));
324: PetscCall(VecNestGetSubVec(Y, 1, &y2));
325: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n"));
326: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
327: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n"));
328: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
329: PetscCall(VecDot(X, Y, &scalar));
331: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
332: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
333: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
335: PetscCall(VecMax(X, &index, &real));
336: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
337: PetscCall(VecMin(X, &index, &real));
338: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
340: PetscCall(VecDestroy(&X));
341: PetscCall(VecDestroy(&Y));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: int main(int argc, char **args)
346: {
347: PetscFunctionBeginUser;
348: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
349: PetscCall(test_view());
350: PetscCall(test_axpy_dot_max());
351: PetscCall(PetscFinalize());
352: return 0;
353: }
355: /*TEST
357: test:
358: args: -explicit_is 0
360: test:
361: suffix: 2
362: args: -explicit_is 1
363: output_file: output/ex37_1.out
365: test:
366: suffix: 3
367: nsize: 2
368: args: -explicit_is 0
370: testset:
371: nsize: 2
372: args: -explicit_is 1
373: output_file: output/ex37_4.out
374: filter: grep -v -e "type: mpi" -e "type=mpi"
376: test:
377: suffix: 4
379: test:
380: requires: cuda
381: suffix: 4_cuda
382: args: -vec_type cuda
384: test:
385: requires: kokkos_kernels
386: suffix: 4_kokkos
387: args: -vec_type kokkos
389: test:
390: requires: hip
391: suffix: 4_hip
392: args: -vec_type hip
394: testset:
395: nsize: 2
396: args: -explicit_is 1 -inv
397: output_file: output/ex37_5.out
398: filter: grep -v -e "type: mpi" -e "type=mpi"
400: test:
401: suffix: 5
403: test:
404: requires: cuda
405: suffix: 5_cuda
406: args: -vec_type cuda
408: test:
409: requires: kokkos_kernels
410: suffix: 5_kokkos
411: args: -vec_type kokkos
413: test:
414: requires: hip
415: suffix: 5_hip
416: args: -vec_type hip
418: TEST*/