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 test_view(void)
23: {
24: Vec X, lX, a, b;
25: Vec c, d, e, f;
26: Vec tmp_buf[2];
27: IS tmp_is[2];
28: PetscInt index, n;
29: PetscReal val;
30: PetscInt list[] = {0, 1, 2};
31: PetscScalar vals[] = {0.5, 0.25, 0.125};
32: PetscScalar *x, *lx;
33: PetscBool explcit = PETSC_FALSE, inv = PETSC_FALSE;
35: PetscFunctionBegin;
36: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
38: PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
39: PetscCall(VecSetSizes(c, PETSC_DECIDE, 3));
40: PetscCall(VecSetFromOptions(c));
41: PetscCall(VecDuplicate(c, &d));
42: PetscCall(VecDuplicate(c, &e));
43: PetscCall(VecDuplicate(c, &f));
45: PetscCall(VecSet(c, 1.0));
46: PetscCall(VecSet(d, 2.0));
47: PetscCall(VecSet(e, 3.0));
48: PetscCall(VecSetValues(f, 3, list, vals, INSERT_VALUES));
49: PetscCall(VecAssemblyBegin(f));
50: PetscCall(VecAssemblyEnd(f));
51: PetscCall(VecScale(f, 10.0));
53: tmp_buf[0] = e;
54: tmp_buf[1] = f;
55: PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicit_is", &explcit, 0));
56: PetscCall(GetISs(tmp_buf, tmp_is, PETSC_FALSE));
57: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, explcit ? tmp_is : NULL, tmp_buf, &b));
58: PetscCall(VecDestroy(&e));
59: PetscCall(VecDestroy(&f));
60: PetscCall(ISDestroy(&tmp_is[0]));
61: PetscCall(ISDestroy(&tmp_is[1]));
63: tmp_buf[0] = c;
64: tmp_buf[1] = d;
65: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &a));
66: PetscCall(VecDestroy(&c));
67: PetscCall(VecDestroy(&d));
69: PetscCall(PetscOptionsGetBool(NULL, NULL, "-inv", &inv, 0));
70: tmp_buf[0] = a;
71: tmp_buf[1] = b;
72: if (inv) {
73: PetscCall(GetISs(tmp_buf, tmp_is, inv));
74: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, tmp_is, tmp_buf, &X));
75: PetscCall(ISDestroy(&tmp_is[0]));
76: PetscCall(ISDestroy(&tmp_is[1]));
77: } else {
78: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
79: }
80: PetscCall(VecDestroy(&a));
82: PetscCall(VecAssemblyBegin(X));
83: PetscCall(VecAssemblyEnd(X));
85: PetscCall(VecMax(b, &index, &val));
86: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
88: PetscCall(VecMin(b, &index, &val));
89: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
91: PetscCall(VecDestroy(&b));
93: PetscCall(VecMax(X, &index, &val));
94: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
95: PetscCall(VecMin(X, &index, &val));
96: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
98: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
100: PetscCall(VecCreateLocalVector(X, &lX));
101: PetscCall(VecGetLocalVectorRead(X, lX));
102: PetscCall(VecGetLocalSize(lX, &n));
103: PetscCall(VecGetArrayRead(lX, (const PetscScalar **)&lx));
104: PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
105: for (PetscInt i = 0; i < n; i++) PetscCheck(lx[i] == x[i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid data");
106: PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&x));
107: PetscCall(VecRestoreArrayRead(lX, (const PetscScalar **)&lx));
108: PetscCall(VecRestoreLocalVectorRead(X, lX));
110: PetscCall(VecDestroy(&lX));
111: PetscCall(VecDestroy(&X));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: #if 0
116: PetscErrorCode test_vec_ops(void)
117: {
118: Vec X, a,b;
119: Vec c,d,e,f;
120: PetscScalar val;
122: PetscFunctionBegin;
123: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME));
125: PetscCall(VecCreate(PETSC_COMM_WORLD, &X));
126: PetscCall(VecSetSizes(X, 2, 2));
127: PetscCall(VecSetType(X, VECNEST));
129: PetscCall(VecCreate(PETSC_COMM_WORLD, &a));
130: PetscCall(VecSetSizes(a, 2, 2));
131: PetscCall(VecSetType(a, VECNEST));
133: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
134: PetscCall(VecSetSizes(b, 2, 2));
135: PetscCall(VecSetType(b, VECNEST));
137: /* assemble X */
138: PetscCall(VecNestSetSubVec(X, 0, a));
139: PetscCall(VecNestSetSubVec(X, 1, b));
140: PetscCall(VecAssemblyBegin(X));
141: PetscCall(VecAssemblyEnd(X));
143: PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
144: PetscCall(VecSetSizes(c, 3, 3));
145: PetscCall(VecSetType(c, VECSEQ));
146: PetscCall(VecDuplicate(c, &d));
147: PetscCall(VecDuplicate(c, &e));
148: PetscCall(VecDuplicate(c, &f));
150: PetscCall(VecSet(c, 1.0));
151: PetscCall(VecSet(d, 2.0));
152: PetscCall(VecSet(e, 3.0));
153: PetscCall(VecSet(f, 4.0));
155: /* assemble a */
156: PetscCall(VecNestSetSubVec(a, 0, c));
157: PetscCall(VecNestSetSubVec(a, 1, d));
158: PetscCall(VecAssemblyBegin(a));
159: PetscCall(VecAssemblyEnd(a));
161: /* assemble b */
162: PetscCall(VecNestSetSubVec(b, 0, e));
163: PetscCall(VecNestSetSubVec(b, 1, f));
164: PetscCall(VecAssemblyBegin(b));
165: PetscCall(VecAssemblyEnd(b));
167: PetscCall(VecDot(X,X, &val));
168: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
171: #endif
173: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
174: {
175: PetscMPIInt size;
176: Vec v;
177: PetscInt i;
178: PetscScalar vx;
180: PetscFunctionBegin;
181: PetscCallMPI(MPI_Comm_size(comm, &size));
182: PetscCall(VecCreate(comm, &v));
183: PetscCall(VecSetSizes(v, PETSC_DECIDE, length));
184: if (size == 1) PetscCall(VecSetType(v, VECSEQ));
185: else PetscCall(VecSetType(v, VECMPI));
187: for (i = 0; i < length; i++) {
188: vx = (PetscScalar)(start_value + i * stride);
189: PetscCall(VecSetValue(v, i, vx, INSERT_VALUES));
190: }
191: PetscCall(VecAssemblyBegin(v));
192: PetscCall(VecAssemblyEnd(v));
194: *_v = v;
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*
199: X = ([0,1,2,3], [10,12,14,16,18])
200: Y = ([4,7,10,13], [5,6,7,8,9])
202: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
203: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
205: */
206: PetscErrorCode test_axpy_dot_max(void)
207: {
208: Vec x1, y1, x2, y2;
209: Vec tmp_buf[2];
210: Vec X, Y;
211: PetscReal real, real2;
212: PetscScalar scalar;
213: PetscInt index;
215: PetscFunctionBegin;
216: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
218: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1));
219: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2));
221: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1));
222: PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2));
224: tmp_buf[0] = x1;
225: tmp_buf[1] = x2;
226: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
227: PetscCall(VecAssemblyBegin(X));
228: PetscCall(VecAssemblyEnd(X));
229: PetscCall(VecDestroy(&x1));
230: PetscCall(VecDestroy(&x2));
232: tmp_buf[0] = y1;
233: tmp_buf[1] = y2;
234: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &Y));
235: PetscCall(VecAssemblyBegin(Y));
236: PetscCall(VecAssemblyEnd(Y));
237: PetscCall(VecDestroy(&y1));
238: PetscCall(VecDestroy(&y2));
240: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n"));
241: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
242: PetscCall(VecNestGetSubVec(Y, 0, &y1));
243: PetscCall(VecNestGetSubVec(Y, 1, &y2));
244: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n"));
245: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
246: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n"));
247: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
248: PetscCall(VecDot(X, Y, &scalar));
250: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
252: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
253: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
255: PetscCall(VecAXPY(Y, 1.0, X)); /* Y <- a X + Y */
256: PetscCall(VecNestGetSubVec(Y, 0, &y1));
257: PetscCall(VecNestGetSubVec(Y, 1, &y2));
258: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n"));
259: PetscCall(VecView(y1, PETSC_VIEWER_STDOUT_WORLD));
260: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n"));
261: PetscCall(VecView(y2, PETSC_VIEWER_STDOUT_WORLD));
262: PetscCall(VecDot(X, Y, &scalar));
264: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
265: PetscCall(VecDotNorm2(X, Y, &scalar, &real2));
266: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2));
268: PetscCall(VecMax(X, &index, &real));
269: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
270: PetscCall(VecMin(X, &index, &real));
271: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index));
273: PetscCall(VecDestroy(&X));
274: PetscCall(VecDestroy(&Y));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: int main(int argc, char **args)
279: {
280: PetscFunctionBeginUser;
281: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
282: PetscCall(test_view());
283: PetscCall(test_axpy_dot_max());
284: PetscCall(PetscFinalize());
285: return 0;
286: }
288: /*TEST
290: test:
291: args: -explicit_is 0
293: test:
294: suffix: 2
295: args: -explicit_is 1
296: output_file: output/ex37_1.out
298: test:
299: suffix: 3
300: nsize: 2
301: args: -explicit_is 0
303: testset:
304: nsize: 2
305: args: -explicit_is 1
306: output_file: output/ex37_4.out
307: filter: grep -v -e "type: mpi" -e "type=mpi"
309: test:
310: suffix: 4
312: test:
313: requires: cuda
314: suffix: 4_cuda
315: args: -vec_type cuda
317: test:
318: requires: kokkos_kernels
319: suffix: 4_kokkos
320: args: -vec_type kokkos
322: test:
323: requires: hip
324: suffix: 4_hip
325: args: -vec_type hip
327: testset:
328: nsize: 2
329: args: -explicit_is 1 -inv
330: output_file: output/ex37_5.out
331: filter: grep -v -e "type: mpi" -e "type=mpi"
333: test:
334: suffix: 5
336: test:
337: requires: cuda
338: suffix: 5_cuda
339: args: -vec_type cuda
341: test:
342: requires: kokkos_kernels
343: suffix: 5_kokkos
344: args: -vec_type kokkos
346: test:
347: requires: hip
348: suffix: 5_hip
349: args: -vec_type hip
351: TEST*/