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: VecGetOwnershipRange(vecs[0], &rstart[0], &rend[0]);
10: VecGetOwnershipRange(vecs[1], &rstart[1], &rend[1]);
11: if (!inv) {
12: ISCreateStride(PETSC_COMM_WORLD, rend[0] - rstart[0], rstart[0] + rstart[1], 1, &is[0]);
13: ISCreateStride(PETSC_COMM_WORLD, rend[1] - rstart[1], rend[0] + rstart[1], 1, &is[1]);
14: } else {
15: ISCreateStride(PETSC_COMM_WORLD, rend[0] - rstart[0], rend[0] + rend[1] - 1, -1, &is[0]);
16: ISCreateStride(PETSC_COMM_WORLD, rend[1] - rstart[1], rstart[0] + rend[1] - 1, -1, &is[1]);
17: }
18: return 0;
19: }
21: PetscErrorCode test_view(void)
22: {
23: Vec X, lX, a, b;
24: Vec c, d, e, f;
25: Vec tmp_buf[2];
26: IS tmp_is[2];
27: PetscInt index, n;
28: PetscReal val;
29: PetscInt list[] = {0, 1, 2};
30: PetscScalar vals[] = {0.5, 0.25, 0.125};
31: PetscScalar *x, *lx;
32: PetscBool explcit = PETSC_FALSE, inv = PETSC_FALSE;
34: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
36: VecCreate(PETSC_COMM_WORLD, &c);
37: VecSetSizes(c, PETSC_DECIDE, 3);
38: VecSetFromOptions(c);
39: VecDuplicate(c, &d);
40: VecDuplicate(c, &e);
41: VecDuplicate(c, &f);
43: VecSet(c, 1.0);
44: VecSet(d, 2.0);
45: VecSet(e, 3.0);
46: VecSetValues(f, 3, list, vals, INSERT_VALUES);
47: VecAssemblyBegin(f);
48: VecAssemblyEnd(f);
49: VecScale(f, 10.0);
51: tmp_buf[0] = e;
52: tmp_buf[1] = f;
53: PetscOptionsGetBool(NULL, NULL, "-explicit_is", &explcit, 0);
54: GetISs(tmp_buf, tmp_is, PETSC_FALSE);
55: VecCreateNest(PETSC_COMM_WORLD, 2, explcit ? tmp_is : NULL, tmp_buf, &b);
56: VecDestroy(&e);
57: VecDestroy(&f);
58: ISDestroy(&tmp_is[0]);
59: ISDestroy(&tmp_is[1]);
61: tmp_buf[0] = c;
62: tmp_buf[1] = d;
63: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &a);
64: VecDestroy(&c);
65: VecDestroy(&d);
67: PetscOptionsGetBool(NULL, NULL, "-inv", &inv, 0);
68: tmp_buf[0] = a;
69: tmp_buf[1] = b;
70: if (inv) {
71: GetISs(tmp_buf, tmp_is, inv);
72: VecCreateNest(PETSC_COMM_WORLD, 2, tmp_is, tmp_buf, &X);
73: ISDestroy(&tmp_is[0]);
74: ISDestroy(&tmp_is[1]);
75: } else {
76: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X);
77: }
78: VecDestroy(&a);
80: VecAssemblyBegin(X);
81: VecAssemblyEnd(X);
83: VecMax(b, &index, &val);
84: PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index);
86: VecMin(b, &index, &val);
87: PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index);
89: VecDestroy(&b);
91: VecMax(X, &index, &val);
92: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index);
93: VecMin(X, &index, &val);
94: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)val, index);
96: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
98: VecCreateLocalVector(X, &lX);
99: VecGetLocalVectorRead(X, lX);
100: VecGetLocalSize(lX, &n);
101: VecGetArrayRead(lX, (const PetscScalar **)&lx);
102: VecGetArrayRead(X, (const PetscScalar **)&x);
104: VecRestoreArrayRead(X, (const PetscScalar **)&x);
105: VecRestoreArrayRead(lX, (const PetscScalar **)&lx);
106: VecRestoreLocalVectorRead(X, lX);
108: VecDestroy(&lX);
109: VecDestroy(&X);
110: return 0;
111: }
113: #if 0
114: PetscErrorCode test_vec_ops(void)
115: {
116: Vec X, a,b;
117: Vec c,d,e,f;
118: PetscScalar val;
120: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);
122: VecCreate(PETSC_COMM_WORLD, &X);
123: VecSetSizes(X, 2, 2);
124: VecSetType(X, VECNEST);
126: VecCreate(PETSC_COMM_WORLD, &a);
127: VecSetSizes(a, 2, 2);
128: VecSetType(a, VECNEST);
130: VecCreate(PETSC_COMM_WORLD, &b);
131: VecSetSizes(b, 2, 2);
132: VecSetType(b, VECNEST);
134: /* assemble X */
135: VecNestSetSubVec(X, 0, a);
136: VecNestSetSubVec(X, 1, b);
137: VecAssemblyBegin(X);
138: VecAssemblyEnd(X);
140: VecCreate(PETSC_COMM_WORLD, &c);
141: VecSetSizes(c, 3, 3);
142: VecSetType(c, VECSEQ);
143: VecDuplicate(c, &d);
144: VecDuplicate(c, &e);
145: VecDuplicate(c, &f);
147: VecSet(c, 1.0);
148: VecSet(d, 2.0);
149: VecSet(e, 3.0);
150: VecSet(f, 4.0);
152: /* assemble a */
153: VecNestSetSubVec(a, 0, c);
154: VecNestSetSubVec(a, 1, d);
155: VecAssemblyBegin(a);
156: VecAssemblyEnd(a);
158: /* assemble b */
159: VecNestSetSubVec(b, 0, e);
160: VecNestSetSubVec(b, 1, f);
161: VecAssemblyBegin(b);
162: VecAssemblyEnd(b);
164: VecDot(X,X, &val);
165: PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
166: return 0;
167: }
168: #endif
170: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
171: {
172: PetscMPIInt size;
173: Vec v;
174: PetscInt i;
175: PetscScalar vx;
177: MPI_Comm_size(comm, &size);
178: VecCreate(comm, &v);
179: VecSetSizes(v, PETSC_DECIDE, length);
180: if (size == 1) VecSetType(v, VECSEQ);
181: else VecSetType(v, VECMPI);
183: for (i = 0; i < length; i++) {
184: vx = (PetscScalar)(start_value + i * stride);
185: VecSetValue(v, i, vx, INSERT_VALUES);
186: }
187: VecAssemblyBegin(v);
188: VecAssemblyEnd(v);
190: *_v = v;
191: return 0;
192: }
194: /*
195: X = ([0,1,2,3], [10,12,14,16,18])
196: Y = ([4,7,10,13], [5,6,7,8,9])
198: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
199: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])
201: */
202: PetscErrorCode test_axpy_dot_max(void)
203: {
204: Vec x1, y1, x2, y2;
205: Vec tmp_buf[2];
206: Vec X, Y;
207: PetscReal real, real2;
208: PetscScalar scalar;
209: PetscInt index;
211: PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);
213: gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
214: gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);
216: gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
217: gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);
219: tmp_buf[0] = x1;
220: tmp_buf[1] = x2;
221: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X);
222: VecAssemblyBegin(X);
223: VecAssemblyEnd(X);
224: VecDestroy(&x1);
225: VecDestroy(&x2);
227: tmp_buf[0] = y1;
228: tmp_buf[1] = y2;
229: VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &Y);
230: VecAssemblyBegin(Y);
231: VecAssemblyEnd(Y);
232: VecDestroy(&y1);
233: VecDestroy(&y2);
235: PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
236: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
237: VecNestGetSubVec(Y, 0, &y1);
238: VecNestGetSubVec(Y, 1, &y2);
239: PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
240: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
241: PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
242: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
243: VecDot(X, Y, &scalar);
245: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
247: VecDotNorm2(X, Y, &scalar, &real2);
248: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
250: VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
251: VecNestGetSubVec(Y, 0, &y1);
252: VecNestGetSubVec(Y, 1, &y2);
253: PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
254: VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
255: PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
256: VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
257: VecDot(X, Y, &scalar);
259: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
260: VecDotNorm2(X, Y, &scalar, &real2);
261: PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);
263: VecMax(X, &index, &real);
264: PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index);
265: VecMin(X, &index, &real);
266: PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n", (double)real, index);
268: VecDestroy(&X);
269: VecDestroy(&Y);
270: return 0;
271: }
273: int main(int argc, char **args)
274: {
276: PetscInitialize(&argc, &args, (char *)0, help);
277: test_view();
278: test_axpy_dot_max();
279: PetscFinalize();
280: return 0;
281: }
283: /*TEST
285: test:
286: args: -explicit_is 0
288: test:
289: suffix: 2
290: args: -explicit_is 1
291: output_file: output/ex37_1.out
293: test:
294: suffix: 3
295: nsize: 2
296: args: -explicit_is 0
298: testset:
299: nsize: 2
300: args: -explicit_is 1
301: output_file: output/ex37_4.out
302: filter: grep -v -e "type: mpi" -e "type=mpi"
304: test:
305: suffix: 4
307: test:
308: requires: cuda
309: suffix: 4_cuda
310: args: -vec_type cuda
312: test:
313: requires: kokkos_kernels
314: suffix: 4_kokkos
315: args: -vec_type kokkos
317: test:
318: requires: hip
319: suffix: 4_hip
320: args: -vec_type hip
322: testset:
323: nsize: 2
324: args: -explicit_is 1 -inv
325: output_file: output/ex37_5.out
326: filter: grep -v -e "type: mpi" -e "type=mpi"
328: test:
329: suffix: 5
331: test:
332: requires: cuda
333: suffix: 5_cuda
334: args: -vec_type cuda
336: test:
337: requires: kokkos_kernels
338: suffix: 5_kokkos
339: args: -vec_type kokkos
341: test:
342: requires: hip
343: suffix: 5_hip
344: args: -vec_type hip
346: TEST*/