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, NULL));
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, NULL));
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: {
133: PetscReal norm[2], ref[2];
135: PetscCall(VecNorm(A, NORM_FROBENIUS, norm));
136: PetscCall(VecNorm(vA, NORM_FROBENIUS, ref));
137: PetscCheck(PetscIsCloseAtTol(norm[0], ref[0], PETSC_SMALL, PETSC_SMALL), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Nested Frobenius norm differs from the CPU reference");
138: PetscCall(VecNorm(A, NORM_1_AND_2, norm));
139: PetscCall(VecNorm(vA, NORM_1_AND_2, ref));
140: PetscCheck(PetscIsCloseAtTol(norm[0], ref[0], PETSC_SMALL, PETSC_SMALL), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Nested 1-norm differs from the CPU reference");
141: PetscCheck(PetscIsCloseAtTol(norm[1], ref[1], PETSC_SMALL, PETSC_SMALL), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Nested 2-norm differs from the CPU reference");
142: }
143: PetscCall(VecCopy(X, X2));
144: PetscCall(VecDuplicate(vX, &vX2));
145: PetscCall(VecCopy(vX, vX2));
146: PetscCall(VecScale(X2, 2.0));
147: PetscCall(VecScale(vX2, 2.0));
148: for (int nt = 0; nt < 2; nt++) {
149: NormType norm = nt ? NORM_INFINITY : NORM_2;
150: for (int e = 0; e < 2; e++) {
151: for (int a = 0; a < 2; a++) {
152: for (int r = 0; r < 2; r++) {
153: PetscReal vn, vna, vnr, nn, nna, nnr;
154: PetscInt vn_loc, vna_loc, vnr_loc, nn_loc, nna_loc, nnr_loc;
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing Wnorms %s: E? %d A? %d R? %d\n", norm == NORM_2 ? "2" : "inf", e, a, r));
157: 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));
158: 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));
159: if (vn_loc != nn_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with total norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vn_loc, nn_loc));
160: if (vna_loc != nna_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with absolute norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vna_loc, nna_loc));
161: if (vnr_loc != nnr_loc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " error with relative norm loc: %" PetscInt_FMT " %" PetscInt_FMT "\n", vnr_loc, nnr_loc));
162: 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)));
163: 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)));
164: }
165: }
166: }
167: }
168: PetscCall(VecDestroy(&X2));
169: PetscCall(VecDestroy(&A));
170: PetscCall(VecDestroy(&R));
171: PetscCall(VecDestroy(&E));
172: PetscCall(VecDestroy(&vX));
173: PetscCall(VecDestroy(&vX2));
174: PetscCall(VecDestroy(&vA));
175: PetscCall(VecDestroy(&vR));
176: PetscCall(VecDestroy(&vE));
178: PetscCall(VecCreateLocalVector(X, &lX));
179: PetscCall(VecGetLocalVectorRead(X, lX));
180: PetscCall(VecGetLocalSize(lX, &n));
181: PetscCall(VecGetArrayRead(lX, (const PetscScalar **)&lx));
182: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
183: for (PetscInt i = 0; i < n; i++) PetscCheck(lx[i] == x[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid data");
184: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
185: PetscCall(VecRestoreArrayRead(lX, (const PetscScalar **)&lx));
186: PetscCall(VecRestoreLocalVectorRead(X, lX));
188: PetscCall(VecDestroy(&lX));
189: PetscCall(VecDestroy(&X));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: #if 0
194: PetscErrorCode test_vec_ops(void)
195: {
196: Vec X, a,b;
197: Vec c,d,e,f;
198: PetscScalar val;
200: PetscFunctionBegin;
201: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME));
203: PetscCall(VecCreate(PETSC_COMM_WORLD, &X));
204: PetscCall(VecSetSizes(X, 2, 2));
205: PetscCall(VecSetType(X, VECNEST));
207: PetscCall(VecCreate(PETSC_COMM_WORLD, &a));
208: PetscCall(VecSetSizes(a, 2, 2));
209: PetscCall(VecSetType(a, VECNEST));
211: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
212: PetscCall(VecSetSizes(b, 2, 2));
213: PetscCall(VecSetType(b, VECNEST));
215: /* assemble X */
216: PetscCall(VecNestSetSubVec(X, 0, a));
217: PetscCall(VecNestSetSubVec(X, 1, b));
218: PetscCall(VecAssemblyBegin(X));
219: PetscCall(VecAssemblyEnd(X));
221: PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
222: PetscCall(VecSetSizes(c, 3, 3));
223: PetscCall(VecSetType(c, VECSEQ));
224: PetscCall(VecDuplicate(c, &d));
225: PetscCall(VecDuplicate(c, &e));
226: PetscCall(VecDuplicate(c, &f));
228: PetscCall(VecSet(c, 1.0));
229: PetscCall(VecSet(d, 2.0));
230: PetscCall(VecSet(e, 3.0));
231: PetscCall(VecSet(f, 4.0));
233: /* assemble a */
234: PetscCall(VecNestSetSubVec(a, 0, c));
235: PetscCall(VecNestSetSubVec(a, 1, d));
236: PetscCall(VecAssemblyBegin(a));
237: PetscCall(VecAssemblyEnd(a));
239: /* assemble b */
240: PetscCall(VecNestSetSubVec(b, 0, e));
241: PetscCall(VecNestSetSubVec(b, 1, f));
242: PetscCall(VecAssemblyBegin(b));
243: PetscCall(VecAssemblyEnd(b));
245: PetscCall(VecDot(X,X, &val));
246: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
249: #endif
251: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
252: {
253: PetscMPIInt size;
254: Vec v;
255: PetscScalar vx;
257: PetscFunctionBegin;
258: PetscCallMPI(MPI_Comm_size(comm, &size));
259: PetscCall(VecCreate(comm, &v));
260: PetscCall(VecSetSizes(v, PETSC_DECIDE, length));
261: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
262: else PetscCall(VecSetType(v, VECMPI));
264: for (PetscInt i = 0; i < length; i++) {
265: vx = (PetscScalar)(start_value + i * stride);
266: PetscCall(VecSetValue(v, i, vx, INSERT_VALUES));
267: }
268: PetscCall(VecAssemblyBegin(v));
269: PetscCall(VecAssemblyEnd(v));
271: *_v = v;
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*
276: X = ([0,1,2,3], [10,12,14,16,18])
277: Y = ([4,7,10,13], [5,6,7,8,9])
279: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
280: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
282: */
283: PetscErrorCode test_axpy_dot_max(void)
284: {
285: Vec x1, y1, x2, y2;
286: Vec tmp_buf[2];
287: Vec X, Y;
288: PetscReal real, real2;
289: PetscScalar scalar;
290: PetscInt index;
292: PetscFunctionBegin;
293: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
295: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1));
296: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2));
298: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1));
299: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2));
301: tmp_buf[0] = x1;
302: tmp_buf[1] = x2;
303: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
304: PetscCall(VecAssemblyBegin(X));
305: PetscCall(VecAssemblyEnd(X));
306: PetscCall(VecDestroy(&x1));
307: PetscCall(VecDestroy(&x2));
309: tmp_buf[0] = y1;
310: tmp_buf[1] = y2;
311: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &Y));
312: PetscCall(VecAssemblyBegin(Y));
313: PetscCall(VecAssemblyEnd(Y));
314: PetscCall(VecDestroy(&y1));
315: PetscCall(VecDestroy(&y2));
317: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n"));
318: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
319: PetscCall(VecNestGetSubVec(Y, 0, &y1));
320: PetscCall(VecNestGetSubVec(Y, 1, &y2));
321: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n"));
322: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
323: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n"));
324: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
325: PetscCall(VecDot(X, Y, &scalar));
327: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
329: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
330: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
332: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
333: PetscCall(VecNestGetSubVec(Y, 0, &y1));
334: PetscCall(VecNestGetSubVec(Y, 1, &y2));
335: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n"));
336: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
337: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n"));
338: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
339: PetscCall(VecDot(X, Y, &scalar));
341: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
342: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
343: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
345: PetscCall(VecMax(X, &index, &real));
346: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
347: PetscCall(VecMin(X, &index, &real));
348: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
350: PetscCall(VecDestroy(&X));
351: PetscCall(VecDestroy(&Y));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: int main(int argc, char **args)
356: {
357: PetscFunctionBeginUser;
358: PetscCall(PetscInitialize(&argc, &args, NULL, help));
359: PetscCall(test_view());
360: PetscCall(test_axpy_dot_max());
361: PetscCall(PetscFinalize());
362: return 0;
363: }
365: /*TEST
367: test:
368: args: -explicit_is 0
370: test:
371: suffix: 2
372: args: -explicit_is 1
373: output_file: output/ex37_1.out
375: test:
376: suffix: 3
377: nsize: 2
378: args: -explicit_is 0
380: testset:
381: nsize: 2
382: args: -explicit_is 1
383: output_file: output/ex37_4.out
384: filter: grep -v -e "type: mpi" -e "type=mpi"
386: test:
387: suffix: 4
389: test:
390: requires: cuda
391: suffix: 4_cuda
392: args: -vec_type cuda
394: test:
395: requires: kokkos_kernels
396: suffix: 4_kokkos
397: args: -vec_type kokkos
399: test:
400: requires: hip
401: suffix: 4_hip
402: args: -vec_type hip
404: testset:
405: nsize: 2
406: args: -explicit_is 1 -inv
407: output_file: output/ex37_5.out
408: filter: grep -v -e "type: mpi" -e "type=mpi"
410: test:
411: suffix: 5
413: test:
414: requires: cuda
415: suffix: 5_cuda
416: args: -vec_type cuda
418: test:
419: requires: kokkos_kernels
420: suffix: 5_kokkos
421: args: -vec_type kokkos
423: test:
424: requires: hip
425: suffix: 5_hip
426: args: -vec_type hip
428: TEST*/