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

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

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

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

 41:   VecCreate(PETSC_COMM_WORLD, &c);
 42:   VecSetSizes(c, PETSC_DECIDE, 3);
 43:   VecSetFromOptions(c);
 44:   VecDuplicate(c, &d);
 45:   VecDuplicate(c, &e);
 46:   VecDuplicate(c, &f);

 48:   VecSet(c, 1.0);
 49:   VecSet(d, 2.0);
 50:   VecSet(e, 3.0);
 51:   VecSetValues(f,3,list,vals,INSERT_VALUES);
 52:   VecAssemblyBegin(f);
 53:   VecAssemblyEnd(f);
 54:   VecScale(f, 10.0);

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

 66:   tmp_buf[0] = c;
 67:   tmp_buf[1] = d;
 68:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&a);
 69:   VecDestroy(&c);   VecDestroy(&d);

 71:   tmp_buf[0] = a;
 72:   tmp_buf[1] = b;
 73:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
 74:   VecDestroy(&a);

 76:   VecAssemblyBegin(X);
 77:   VecAssemblyEnd(X);

 79:   VecMax(b, &index, &val);
 80:   PetscPrintf(PETSC_COMM_WORLD, "(max-b) = %f : index = %D \n",(double) val, index);

 82:   VecMin(b, &index, &val);
 83:   PetscPrintf(PETSC_COMM_WORLD, "(min-b) = %f : index = %D \n",(double) val, index);

 85:   VecDestroy(&b);

 87:   VecMax(X, &index, &val);
 88:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %D \n",(double) val, index);
 89:   VecMin(X, &index, &val);
 90:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %D \n",(double) val, index);

 92:   VecView(X, PETSC_VIEWER_STDOUT_WORLD);

 94:   VecDestroy(&X);
 95:   return(0);
 96: }

 98: #if 0
 99: PetscErrorCode test_vec_ops(void)
100: {
101:   Vec            X, a,b;
102:   Vec            c,d,e,f;
103:   PetscScalar    val;

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

109:   VecCreate(PETSC_COMM_WORLD, &X);
110:   VecSetSizes(X, 2, 2);
111:   VecSetType(X, VECNEST);

113:   VecCreate(PETSC_COMM_WORLD, &a);
114:   VecSetSizes(a, 2, 2);
115:   VecSetType(a, VECNEST);

117:   VecCreate(PETSC_COMM_WORLD, &b);
118:   VecSetSizes(b, 2, 2);
119:   VecSetType(b, VECNEST);

121:   /* assemble X */
122:   VecNestSetSubVec(X, 0, a);
123:   VecNestSetSubVec(X, 1, b);
124:   VecAssemblyBegin(X);
125:   VecAssemblyEnd(X);

127:   VecCreate(PETSC_COMM_WORLD, &c);
128:   VecSetSizes(c, 3, 3);
129:   VecSetType(c, VECSEQ);
130:   VecDuplicate(c, &d);
131:   VecDuplicate(c, &e);
132:   VecDuplicate(c, &f);

134:   VecSet(c, 1.0);
135:   VecSet(d, 2.0);
136:   VecSet(e, 3.0);
137:   VecSet(f, 4.0);

139:   /* assemble a */
140:   VecNestSetSubVec(a, 0, c);
141:   VecNestSetSubVec(a, 1, d);
142:   VecAssemblyBegin(a);
143:   VecAssemblyEnd(a);

145:   /* assemble b */
146:   VecNestSetSubVec(b, 0, e);
147:   VecNestSetSubVec(b, 1, f);
148:   VecAssemblyBegin(b);
149:   VecAssemblyEnd(b);

151:   VecDot(X,X, &val);
152:   PetscPrintf(PETSC_COMM_WORLD, "X.X = %f \n",(double) val);
153:   return(0);
154: }
155: #endif

157: PetscErrorCode gen_test_vector(MPI_Comm comm, PetscInt length, PetscInt start_value, PetscInt stride, Vec *_v)
158: {
159:   int            size;
160:   Vec            v;
161:   PetscInt       i;
162:   PetscScalar    vx;

166:   MPI_Comm_size(comm, &size);
167:   VecCreate(comm, &v);
168:   VecSetSizes(v, PETSC_DECIDE, length);
169:   if (size == 1) { VecSetType(v, VECSEQ); }
170:   else { VecSetType(v, VECMPI); }

172:   for (i=0; i<length; i++) {
173:     vx   = (PetscScalar)(start_value + i * stride);
174:     VecSetValue(v, i, vx, INSERT_VALUES);
175:   }
176:   VecAssemblyBegin(v);
177:   VecAssemblyEnd(v);

179:   *_v = v;
180:   return(0);
181: }

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

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

190: */
191: PetscErrorCode test_axpy_dot_max(void)
192: {
193:   Vec            x1,y1, x2,y2;
194:   Vec            tmp_buf[2];
195:   Vec            X, Y;
196:   PetscReal      real,real2;
197:   PetscScalar    scalar;
198:   PetscInt       index;

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

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

207:   gen_test_vector(PETSC_COMM_WORLD, 4, 4, 3, &y1);
208:   gen_test_vector(PETSC_COMM_WORLD, 5, 5, 1, &y2);

210:   tmp_buf[0] = x1;
211:   tmp_buf[1] = x2;
212:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&X);
213:   VecAssemblyBegin(X);
214:   VecAssemblyEnd(X);
215:   VecDestroy(&x1);
216:   VecDestroy(&x2);

218:   tmp_buf[0] = y1;
219:   tmp_buf[1] = y2;
220:   VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_buf,&Y);
221:   VecAssemblyBegin(Y);
222:   VecAssemblyEnd(Y);
223:   VecDestroy(&y1);
224:   VecDestroy(&y2);

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

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

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

241:   VecAXPY(Y, 1.0, X); /* Y <- a X + Y */
242:   VecNestGetSubVec(Y, 0, &y1);
243:   VecNestGetSubVec(Y, 1, &y2);
244:   PetscPrintf(PETSC_COMM_WORLD, "(2) y1 = \n");
245:   VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
246:   PetscPrintf(PETSC_COMM_WORLD, "(2) y2 = \n");
247:   VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
248:   VecDot(X,Y, &scalar);

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

254:   VecMax(X, &index, &real);
255:   PetscPrintf(PETSC_COMM_WORLD, "(max-X) = %f : index = %D \n",(double) real, index);
256:   VecMin(X, &index, &real);
257:   PetscPrintf(PETSC_COMM_WORLD, "(min-X) = %f : index = %D \n",(double) real, index);

259:   VecDestroy(&X);
260:   VecDestroy(&Y);
261:   return(0);
262: }

264: int main(int argc, char **args)
265: {

268:   PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
269:   test_view();
270:   test_axpy_dot_max();
271:   PetscFinalize();
272:   return ierr;
273: }

275: /*TEST

277:    test:
278:       args: -explicit_is 0

280:    test:
281:       suffix: 2
282:       args: -explicit_is 1
283:       output_file: output/ex37_1.out

285:    test:
286:       suffix: 3
287:       nsize: 2
288:       args: -explicit_is 0

290:    testset:
291:       nsize: 2
292:       args: -explicit_is 1
293:       output_file: output/ex37_4.out
294:       filter: grep -v -e "type: mpi" -e "type=mpi"

296:       test:
297:         suffix: 4

299:       test:
300:         requires: kokkos_kernels
301:         suffix: kokkos
302:         args: -vec_type kokkos

304:       test:
305:         requires: hip
306:         suffix: hip
307:         args: -vec_type hip

309: TEST*/