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