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