Actual source code: ex37.c

  1: static char help[] = "Nest vector functionality.\n\n";

  3: /*T
  4:    Concepts: vectors^block operators
  5:    Concepts: vectors^setting values
  6:    Concepts: vectors^local access to
  7:    Processors: n
  8: T*/

 10: #include <petscvec.h>

 12: static PetscErrorCode GetISs(Vec vecs[],IS is[])
 13: {
 14:   PetscInt       rstart[2],rend[2];

 16:   VecGetOwnershipRange(vecs[0],&rstart[0],&rend[0]);
 17:   VecGetOwnershipRange(vecs[1],&rstart[1],&rend[1]);
 18:   ISCreateStride(PETSC_COMM_WORLD,rend[0]-rstart[0],rstart[0]+rstart[1],1,&is[0]);
 19:   ISCreateStride(PETSC_COMM_WORLD,rend[1]-rstart[1],rend[0]+rstart[1],1,&is[1]);
 20:   return 0;
 21: }

 23: PetscErrorCode test_view(void)
 24: {
 25:   Vec            X, a,b;
 26:   Vec            c,d,e,f;
 27:   Vec            tmp_buf[2];
 28:   IS             tmp_is[2];
 29:   PetscInt       index;
 30:   PetscReal      val;
 31:   PetscInt       list[]={0,1,2};
 32:   PetscScalar    vals[]={0.720032,0.061794,0.0100223};
 33:   PetscBool      explcit = PETSC_FALSE;

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

 37:   VecCreate(PETSC_COMM_WORLD, &c);
 38:   VecSetSizes(c, PETSC_DECIDE, 3);
 39:   VecSetFromOptions(c);
 40:   VecDuplicate(c, &d);
 41:   VecDuplicate(c, &e);
 42:   VecDuplicate(c, &f);

 44:   VecSet(c, 1.0);
 45:   VecSet(d, 2.0);
 46:   VecSet(e, 3.0);
 47:   VecSetValues(f,3,list,vals,INSERT_VALUES);
 48:   VecAssemblyBegin(f);
 49:   VecAssemblyEnd(f);
 50:   VecScale(f, 10.0);

 52:   tmp_buf[0] = e;
 53:   tmp_buf[1] = f;
 54:   PetscOptionsGetBool(NULL,0,"-explicit_is",&explcit,0);
 55:   GetISs(tmp_buf,tmp_is);
 56:   VecCreateNest(PETSC_COMM_WORLD,2,explcit ? tmp_is : NULL,tmp_buf,&b);
 57:   VecDestroy(&e);
 58:   VecDestroy(&f);
 59:   ISDestroy(&tmp_is[0]);
 60:   ISDestroy(&tmp_is[1]);

 62:   tmp_buf[0] = c;
 63:   tmp_buf[1] = d;
 64:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
 65:   VecDestroy(&c));   PetscCall(VecDestroy(&d);

 67:   tmp_buf[0] = a;
 68:   tmp_buf[1] = b;
 69:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
 70:   VecDestroy(&a);

 72:   VecAssemblyBegin(X);
 73:   VecAssemblyEnd(X);

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

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

 81:   VecDestroy(&b);

 83:   VecMax(X, &index, &val);
 84:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);
 85:   VecMin(X, &index, &val);
 86:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) val, index);

 88:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

 90:   VecDestroy(&X);
 91:   return 0;
 92: }

 94: #if 0
 95: PetscErrorCode test_vec_ops(void)
 96: {
 97:   Vec            X, a,b;
 98:   Vec            c,d,e,f;
 99:   PetscScalar    val;

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

103:   VecCreate(PETSC_COMM_WORLD, &X);
104:   VecSetSizes(X, 2, 2);
105:   VecSetType(X, VECNEST);

107:   VecCreate(PETSC_COMM_WORLD, &a);
108:   VecSetSizes(a, 2, 2);
109:   VecSetType(a, VECNEST);

111:   VecCreate(PETSC_COMM_WORLD, &b);
112:   VecSetSizes(b, 2, 2);
113:   VecSetType(b, VECNEST);

115:   /* assemble X */
116:   VecNestSetSubVec(X, 0, a);
117:   VecNestSetSubVec(X, 1, b);
118:   VecAssemblyBegin(X);
119:   VecAssemblyEnd(X);

121:   VecCreate(PETSC_COMM_WORLD, &c);
122:   VecSetSizes(c, 3, 3);
123:   VecSetType(c, VECSEQ);
124:   VecDuplicate(c, &d);
125:   VecDuplicate(c, &e);
126:   VecDuplicate(c, &f);

128:   VecSet(c, 1.0);
129:   VecSet(d, 2.0);
130:   VecSet(e, 3.0);
131:   VecSet(f, 4.0);

133:   /* assemble a */
134:   VecNestSetSubVec(a, 0, c);
135:   VecNestSetSubVec(a, 1, d);
136:   VecAssemblyBegin(a);
137:   VecAssemblyEnd(a);

139:   /* assemble b */
140:   VecNestSetSubVec(b, 0, e);
141:   VecNestSetSubVec(b, 1, f);
142:   VecAssemblyBegin(b);
143:   VecAssemblyEnd(b);

145:   VecDot(X,X, &val);
146:   PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
147:   return 0;
148: }
149: #endif

151: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
152: {
153:   int            size;
154:   Vec            v;
155:   PetscInt       i;
156:   PetscScalar    vx;

158:   MPI_Comm_size(comm, &size);
159:   VecCreate(comm, &v);
160:   VecSetSizes(v, PETSC_DECIDE, length);
161:   if (size == 1) VecSetType(v, VECSEQ);
162:   else VecSetType(v, VECMPI);

164:   for (i=0; i<length; i++) {
165:     vx   = (PetscScalar)(start_value + i * stride);
166:     VecSetValue(v, i, vx, INSERT_VALUES);
167:   }
168:   VecAssemblyBegin(v);
169:   VecAssemblyEnd(v);

171:   *_v = v;
172:   return 0;
173: }

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

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

182: */
183: PetscErrorCode test_axpy_dot_max(void)
184: {
185:   Vec            x1,y1, x2,y2;
186:   Vec            tmp_buf[2];
187:   Vec            X, Y;
188:   PetscReal      real,real2;
189:   PetscScalar    scalar;
190:   PetscInt       index;

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

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

197:   gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
198:   gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);

200:   tmp_buf[0] = x1;
201:   tmp_buf[1] = x2;
202:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
203:   VecAssemblyBegin(X);
204:   VecAssemblyEnd(X);
205:   VecDestroy(&x1);
206:   VecDestroy(&x2);

208:   tmp_buf[0] = y1;
209:   tmp_buf[1] = y2;
210:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
211:   VecAssemblyBegin(Y);
212:   VecAssemblyEnd(Y);
213:   VecDestroy(&y1);
214:   VecDestroy(&y2);

216:   PetscPrintf(PETSC_COMM_WORLD, "VecAXPY \n");
217:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
218:   VecNestGetSubVec(Y, 0, &y1);
219:   VecNestGetSubVec(Y, 1, &y2);
220:   PetscPrintf(PETSC_COMM_WORLD, "(1) y1 = \n");
221:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
222:   PetscPrintf(PETSC_COMM_WORLD, "(1) y2 = \n");
223:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
224:   VecDot(X,Y, &scalar);

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

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

231:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
232:   VecNestGetSubVec(Y, 0, &y1);
233:   VecNestGetSubVec(Y, 1, &y2);
234:   PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
235:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
236:   PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
237:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
238:   VecDot(X,Y, &scalar);

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

244:   VecMax(X, &index, &real);
245:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);
246:   VecMin(X, &index, &real);
247:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %" PetscInt_FMT " \n",(double) real, index);

249:   VecDestroy(&X);
250:   VecDestroy(&Y);
251:   return 0;
252: }

254: int main(int argc, char **args)
255: {

257:   PetscInitialize(&argc, &args,(char*)0, help);
258:   test_view();
259:   test_axpy_dot_max();
260:   PetscFinalize();
261:   return 0;
262: }

264: /*TEST

266:    test:
267:       args: -explicit_is 0

269:    test:
270:       suffix: 2
271:       args: -explicit_is 1
272:       output_file: output/ex37_1.out

274:    test:
275:       suffix: 3
276:       nsize: 2
277:       args: -explicit_is 0

279:    testset:
280:       nsize: 2
281:       args: -explicit_is 1
282:       output_file: output/ex37_4.out
283:       filter: grep -v -e "type: mpi" -e "type=mpi"

285:       test:
286:         suffix: 4

288:       test:
289:         requires: kokkos_kernels
290:         suffix: kokkos
291:         args: -vec_type kokkos

293:       test:
294:         requires: hip
295:         suffix: hip
296:         args: -vec_type hip

298: TEST*/