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[])
  6: {
  7:   PetscInt       rstart[2],rend[2];

  9:   VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
 10:   VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
 11:   ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
 12:   ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
 13:   return 0;
 14: }

 16: PetscErrorCode test_view(void)
 17: {
 18:   Vec            X, a,b;
 19:   Vec            c,d,e,f;
 20:   Vec            tmp_buf[2];
 21:   IS             tmp_is[2];
 22:   PetscInt       index;
 23:   PetscReal      val;
 24:   PetscInt       list[]={0,1,2};
 25:   PetscScalar    vals[]={0.720032,0.061794,0.0100223};
 26:   PetscBool      explcit = PETSC_FALSE;

 28:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);

 30:   VecCreate(PETSC_COMM_WORLD, &c);
 31:   VecSetSizes(c, PETSC_DECIDE, 3);
 32:   VecSetFromOptions(c);
 33:   VecDuplicate(c, &d);
 34:   VecDuplicate(c, &e);
 35:   VecDuplicate(c, &f);

 37:   VecSet(c, 1.0);
 38:   VecSet(d, 2.0);
 39:   VecSet(e, 3.0);
 40:   VecSetValues(f,3,list,vals,INSERT_VALUES);
 41:   VecAssemblyBegin(f);
 42:   VecAssemblyEnd(f);
 43:   VecScale(f, 10.0);

 45:   tmp_buf[0] = e;
 46:   tmp_buf[1] = f;
 47:   PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
 48:   GetISs(tmp_buf,tmp_is);
 49:   VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
 50:   VecDestroy(&e);
 51:   VecDestroy(&f);
 52:   ISDestroy(&tmp_is[0]);
 53:   ISDestroy(&tmp_is[1]);

 55:   tmp_buf[0] = c;
 56:   tmp_buf[1] = d;
 57:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
 58:   VecDestroy(&c));   PetscCall(VecDestroy(&d);

 60:   tmp_buf[0] = a;
 61:   tmp_buf[1] = b;
 62:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
 63:   VecDestroy(&a);

 65:   VecAssemblyBegin(X);
 66:   VecAssemblyEnd(X);

 68:   VecMax(b, &index, &val);
 69:   PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);

 71:   VecMin(b, &index, &val);
 72:   PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %" PetscInt_FMT " \n",(double) val, index);

 74:   VecDestroy(&b);

 76:   VecMax(X, &index, &val);
 77:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
 78:   VecMin(X, &index, &val);
 79:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);

 81:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

 83:   VecDestroy(&X);
 84:   return 0;
 85: }

 87: #if 0
 88: PetscErrorCode test_vec_ops(void)
 89: {
 90:   Vec            X, a,b;
 91:   Vec            c,d,e,f;
 92:   PetscScalar    val;

 94:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n",PETSC_FUNCTION_NAME);

 96:   VecCreate(PETSC_COMM_WORLD, &X);
 97:   VecSetSizes(X, 2, 2);
 98:   VecSetType(X, VECNEST);

100:   VecCreate(PETSC_COMM_WORLD, &a);
101:   VecSetSizes(a, 2, 2);
102:   VecSetType(a, VECNEST);

104:   VecCreate(PETSC_COMM_WORLD, &b);
105:   VecSetSizes(b, 2, 2);
106:   VecSetType(b, VECNEST);

108:   /* assemble X */
109:   VecNestSetSubVec(X, 0, a);
110:   VecNestSetSubVec(X, 1, b);
111:   VecAssemblyBegin(X);
112:   VecAssemblyEnd(X);

114:   VecCreate(PETSC_COMM_WORLD, &c);
115:   VecSetSizes(c, 3, 3);
116:   VecSetType(c, VECSEQ);
117:   VecDuplicate(c, &d);
118:   VecDuplicate(c, &e);
119:   VecDuplicate(c, &f);

121:   VecSet(c, 1.0);
122:   VecSet(d, 2.0);
123:   VecSet(e, 3.0);
124:   VecSet(f, 4.0);

126:   /* assemble a */
127:   VecNestSetSubVec(a, 0, c);
128:   VecNestSetSubVec(a, 1, d);
129:   VecAssemblyBegin(a);
130:   VecAssemblyEnd(a);

132:   /* assemble b */
133:   VecNestSetSubVec(b, 0, e);
134:   VecNestSetSubVec(b, 1, f);
135:   VecAssemblyBegin(b);
136:   VecAssemblyEnd(b);

138:   VecDot(X,X, &val);
139:   PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
140:   return 0;
141: }
142: #endif

144: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
145: {
146:   int            size;
147:   Vec            v;
148:   PetscInt       i;
149:   PetscScalar    vx;

151:   MPI_Comm_size(comm, &size);
152:   VecCreate(comm, &v);
153:   VecSetSizes(v, PETSC_DECIDE, length);
154:   if (size == 1) VecSetType(v, VECSEQ);
155:   else VecSetType(v, VECMPI);

157:   for (i=0; i<length; i++) {
158:     vx   = (PetscScalar)(start_value + i * stride);
159:     VecSetValue(v, i, vx, INSERT_VALUES);
160:   }
161:   VecAssemblyBegin(v);
162:   VecAssemblyEnd(v);

164:   *_v = v;
165:   return 0;
166: }

168: /*
169: X = ([0,1,2,3], [10,12,14,16,18])
170: Y = ([4,7,10,13], [5,6,7,8,9])

172: Y = aX + y = ([4,8,12,16], (15,18,21,24,27])
173: Y = aX + y = ([4,9,14,19], (25,30,35,40,45])

175: */
176: PetscErrorCode test_axpy_dot_max(void)
177: {
178:   Vec            x1,y1, x2,y2;
179:   Vec            tmp_buf[2];
180:   Vec            X, Y;
181:   PetscReal      real,real2;
182:   PetscScalar    scalar;
183:   PetscInt       index;

185:   PetscPrintf(PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME);

187:   gen_test_vector(PETSC_COMM_WORLD, 4, 0, 1, &x1);
188:   gen_test_vector(PETSC_COMM_WORLD, 5, 10, 2, &x2);

190:   gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
191:   gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);

193:   tmp_buf[0] = x1;
194:   tmp_buf[1] = x2;
195:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
196:   VecAssemblyBegin(X);
197:   VecAssemblyEnd(X);
198:   VecDestroy(&x1);
199:   VecDestroy(&x2);

201:   tmp_buf[0] = y1;
202:   tmp_buf[1] = y2;
203:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
204:   VecAssemblyBegin(Y);
205:   VecAssemblyEnd(Y);
206:   VecDestroy(&y1);
207:   VecDestroy(&y2);

209:   PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
210:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
211:   VecNestGetSubVec(Y, 0, &y1);
212:   VecNestGetSubVec(Y, 1, &y2);
213:   PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
214:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
215:   PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
216:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
217:   VecDot(X,Y, &scalar);

219:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));

221:   VecDotNorm2(X,Y, &scalar, &real2);
222:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi     norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);

224:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
225:   VecNestGetSubVec(Y, 0, &y1);
226:   VecNestGetSubVec(Y, 1, &y2);
227:   PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
228:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
229:   PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
230:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
231:   VecDot(X,Y, &scalar);

233:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar));
234:   VecDotNorm2(X,Y, &scalar, &real2);
235:   PetscPrintf(PETSC_COMM_WORLD, "X.Y = %lf + %lfi     norm2(Y) = %lf\n", (double)PetscRealPart(scalar), (double)PetscImaginaryPart(scalar), (double)real2);

237:   VecMax(X, &index, &real);
238:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
239:   VecMin(X, &index, &real);
240:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);

242:   VecDestroy(&X);
243:   VecDestroy(&Y);
244:   return 0;
245: }

247: int main(int argc, char **args)
248: {

250:   PetscInitialize(&argc, &args,(char*)0, help);
251:   test_view();
252:   test_axpy_dot_max();
253:   PetscFinalize();
254:   return 0;
255: }

257: /*TEST

259:    test:
260:       args: -explicit_is 0

262:    test:
263:       suffix: 2
264:       args: -explicit_is 1
265:       output_file: output/ex37_1.out

267:    test:
268:       suffix: 3
269:       nsize: 2
270:       args: -explicit_is 0

272:    testset:
273:       nsize: 2
274:       args: -explicit_is 1
275:       output_file: output/ex37_4.out
276:       filter: grep -v -e "type: mpi" -e "type=mpi"

278:       test:
279:         suffix: 4

281:       test:
282:         requires: kokkos_kernels
283:         suffix: kokkos
284:         args: -vec_type kokkos

286:       test:
287:         requires: hip
288:         suffix: hip
289:         args: -vec_type hip

291: TEST*/