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, NULL, 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*/