1 static char help[] = "Nest vector functionality.\n\n";
2
3 #include <petscvec.h>
4
GetISs(Vec vecs[],IS is[],PetscBool inv)5 static PetscErrorCode GetISs(Vec vecs[], IS is[], PetscBool inv)
6 {
7 PetscInt rstart[2], rend[2];
8
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 }
21
convert_from_nest(Vec X,Vec * Y)22 PetscErrorCode convert_from_nest(Vec X, Vec *Y)
23 {
24 const PetscScalar *v;
25 PetscInt rstart, n, N;
26
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 }
41
test_view(void)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;
54
55 PetscFunctionBegin;
56 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
57
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));
64
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));
72
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]));
82
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));
88
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));
101
102 PetscCall(VecAssemblyBegin(X));
103 PetscCall(VecAssemblyEnd(X));
104
105 PetscCall(VecMax(b, &index, &val));
106 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
107
108 PetscCall(VecMin(b, &index, &val));
109 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n", (double)val, index));
110
111 PetscCall(VecDestroy(&b));
112
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));
117
118 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
119
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;
144
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));
166
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));
176
177 PetscCall(VecDestroy(&lX));
178 PetscCall(VecDestroy(&X));
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
182 #if 0
183 PetscErrorCode test_vec_ops(void)
184 {
185 Vec X, a,b;
186 Vec c,d,e,f;
187 PetscScalar val;
188
189 PetscFunctionBegin;
190 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME));
191
192 PetscCall(VecCreate(PETSC_COMM_WORLD, &X));
193 PetscCall(VecSetSizes(X, 2, 2));
194 PetscCall(VecSetType(X, VECNEST));
195
196 PetscCall(VecCreate(PETSC_COMM_WORLD, &a));
197 PetscCall(VecSetSizes(a, 2, 2));
198 PetscCall(VecSetType(a, VECNEST));
199
200 PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
201 PetscCall(VecSetSizes(b, 2, 2));
202 PetscCall(VecSetType(b, VECNEST));
203
204 /* assemble X */
205 PetscCall(VecNestSetSubVec(X, 0, a));
206 PetscCall(VecNestSetSubVec(X, 1, b));
207 PetscCall(VecAssemblyBegin(X));
208 PetscCall(VecAssemblyEnd(X));
209
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));
216
217 PetscCall(VecSet(c, 1.0));
218 PetscCall(VecSet(d, 2.0));
219 PetscCall(VecSet(e, 3.0));
220 PetscCall(VecSet(f, 4.0));
221
222 /* assemble a */
223 PetscCall(VecNestSetSubVec(a, 0, c));
224 PetscCall(VecNestSetSubVec(a, 1, d));
225 PetscCall(VecAssemblyBegin(a));
226 PetscCall(VecAssemblyEnd(a));
227
228 /* assemble b */
229 PetscCall(VecNestSetSubVec(b, 0, e));
230 PetscCall(VecNestSetSubVec(b, 1, f));
231 PetscCall(VecAssemblyBegin(b));
232 PetscCall(VecAssemblyEnd(b));
233
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
239
gen_test_vector(MPI_Comm comm,PetscInt length,PetscInt start_value,PetscInt stride,Vec * _v)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;
246
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));
253
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));
260
261 *_v = v;
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
265 /*
266 X = ([0,1,2,3], [10,12,14,16,18])
267 Y = ([4,7,10,13], [5,6,7,8,9])
268
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])
271
272 */
test_axpy_dot_max(void)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;
281
282 PetscFunctionBegin;
283 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME));
284
285 PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1));
286 PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2));
287
288 PetscCall(gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1));
289 PetscCall(gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2));
290
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));
298
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));
306
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));
316
317 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar)));
318
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));
321
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));
330
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));
334
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));
339
340 PetscCall(VecDestroy(&X));
341 PetscCall(VecDestroy(&Y));
342 PetscFunctionReturn(PETSC_SUCCESS);
343 }
344
main(int argc,char ** args)345 int main(int argc, char **args)
346 {
347 PetscFunctionBeginUser;
348 PetscCall(PetscInitialize(&argc, &args, NULL, help));
349 PetscCall(test_view());
350 PetscCall(test_axpy_dot_max());
351 PetscCall(PetscFinalize());
352 return 0;
353 }
354
355 /*TEST
356
357 test:
358 args: -explicit_is 0
359
360 test:
361 suffix: 2
362 args: -explicit_is 1
363 output_file: output/ex37_1.out
364
365 test:
366 suffix: 3
367 nsize: 2
368 args: -explicit_is 0
369
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"
375
376 test:
377 suffix: 4
378
379 test:
380 requires: cuda
381 suffix: 4_cuda
382 args: -vec_type cuda
383
384 test:
385 requires: kokkos_kernels
386 suffix: 4_kokkos
387 args: -vec_type kokkos
388
389 test:
390 requires: hip
391 suffix: 4_hip
392 args: -vec_type hip
393
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"
399
400 test:
401 suffix: 5
402
403 test:
404 requires: cuda
405 suffix: 5_cuda
406 args: -vec_type cuda
407
408 test:
409 requires: kokkos_kernels
410 suffix: 5_kokkos
411 args: -vec_type kokkos
412
413 test:
414 requires: hip
415 suffix: 5_hip
416 args: -vec_type hip
417
418 TEST*/
419