Actual source code: ex96.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  2: static char help[] ="Tests sequential and parallel DMCreateMatrix(), MatMatMult() and MatPtAP()\n\
  3:   -Mx <xg>, where <xg> = number of coarse grid points in the x-direction\n\
  4:   -My <yg>, where <yg> = number of coarse grid points in the y-direction\n\
  5:   -Mz <zg>, where <zg> = number of coarse grid points in the z-direction\n\
  6:   -Npx <npx>, where <npx> = number of processors in the x-direction\n\
  7:   -Npy <npy>, where <npy> = number of processors in the y-direction\n\
  8:   -Npz <npz>, where <npz> = number of processors in the z-direction\n\n";

 10: /*
 11:     This test is modified from ~src/ksp/examples/tests/ex19.c.
 12:     Example of usage: mpiexec -n 3 ./ex96 -Mx 10 -My 10 -Mz 10
 13: */

 15:  #include <petscdm.h>
 16:  #include <petscdmda.h>

 18: /* User-defined application contexts */
 19: typedef struct {
 20:   PetscInt mx,my,mz;            /* number grid points in x, y and z direction */
 21:   Vec      localX,localF;       /* local vectors with ghost region */
 22:   DM       da;
 23:   Vec      x,b,r;               /* global vectors */
 24:   Mat      J;                   /* Jacobian on grid */
 25: } GridCtx;
 26: typedef struct {
 27:   GridCtx  fine;
 28:   GridCtx  coarse;
 29:   PetscInt ratio;
 30:   Mat      Ii;                  /* interpolation from coarse to fine */
 31: } AppCtx;

 33: #define COARSE_LEVEL 0
 34: #define FINE_LEVEL   1

 36: /*
 37:       Mm_ratio - ration of grid lines between fine and coarse grids.
 38: */
 39: int main(int argc,char **argv)
 40: {
 42:   AppCtx         user;
 43:   PetscInt       Npx=PETSC_DECIDE,Npy=PETSC_DECIDE,Npz=PETSC_DECIDE;
 44:   PetscMPIInt    size,rank;
 45:   PetscInt       m,n,M,N,i,nrows;
 46:   PetscScalar    one = 1.0;
 47:   PetscReal      fill=2.0;
 48:   Mat            A,A_tmp,P,C,C1,C2;
 49:   PetscScalar    *array,none = -1.0,alpha;
 50:   Vec            x,v1,v2,v3,v4;
 51:   PetscReal      norm,norm_tmp,norm_tmp1,tol=100.*PETSC_MACHINE_EPSILON;
 52:   PetscRandom    rdm;
 53:   PetscBool      Test_MatMatMult=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_3D=PETSC_TRUE,flg;
 54:   const PetscInt *ia,*ja;

 56:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 57:   PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL);

 59:   user.ratio     = 2;
 60:   user.coarse.mx = 20; user.coarse.my = 20; user.coarse.mz = 20;

 62:   PetscOptionsGetInt(NULL,NULL,"-Mx",&user.coarse.mx,NULL);
 63:   PetscOptionsGetInt(NULL,NULL,"-My",&user.coarse.my,NULL);
 64:   PetscOptionsGetInt(NULL,NULL,"-Mz",&user.coarse.mz,NULL);
 65:   PetscOptionsGetInt(NULL,NULL,"-ratio",&user.ratio,NULL);

 67:   if (user.coarse.mz) Test_3D = PETSC_TRUE;

 69:   user.fine.mx = user.ratio*(user.coarse.mx-1)+1;
 70:   user.fine.my = user.ratio*(user.coarse.my-1)+1;
 71:   user.fine.mz = user.ratio*(user.coarse.mz-1)+1;

 73:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 74:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 75:   PetscOptionsGetInt(NULL,NULL,"-Npx",&Npx,NULL);
 76:   PetscOptionsGetInt(NULL,NULL,"-Npy",&Npy,NULL);
 77:   PetscOptionsGetInt(NULL,NULL,"-Npz",&Npz,NULL);

 79:   /* Set up distributed array for fine grid */
 80:   if (!Test_3D) {
 81:     DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,user.fine.my,Npx,Npy,1,1,NULL,NULL,&user.fine.da);
 82:   } else {
 83:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.fine.mx,user.fine.my,user.fine.mz,Npx,Npy,Npz,1,1,NULL,NULL,NULL,&user.fine.da);
 84:   }
 85:   DMSetFromOptions(user.fine.da);
 86:   DMSetUp(user.fine.da);

 88:   /* Test DMCreateMatrix()                                         */
 89:   /*------------------------------------------------------------*/
 90:   DMSetMatType(user.fine.da,MATAIJ);
 91:   DMCreateMatrix(user.fine.da,&A);
 92:   DMSetMatType(user.fine.da,MATBAIJ);
 93:   DMCreateMatrix(user.fine.da,&C);

 95:   MatConvert(C,MATAIJ,MAT_INITIAL_MATRIX,&A_tmp); /* not work for mpisbaij matrix! */
 96:   MatEqual(A,A_tmp,&flg);
 97:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"A != C");
 98:   MatDestroy(&C);
 99:   MatDestroy(&A_tmp);

101:   /*------------------------------------------------------------*/

103:   MatGetLocalSize(A,&m,&n);
104:   MatGetSize(A,&M,&N);
105:   /* if (!rank) printf("A %d, %d\n",M,N); */

107:   /* set val=one to A */
108:   if (size == 1) {
109:     MatGetRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
110:     if (flg) {
111:       MatSeqAIJGetArray(A,&array);
112:       for (i=0; i<ia[nrows]; i++) array[i] = one;
113:       MatSeqAIJRestoreArray(A,&array);
114:     }
115:     MatRestoreRowIJ(A,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
116:   } else {
117:     Mat AA,AB;
118:     MatMPIAIJGetSeqAIJ(A,&AA,&AB,NULL);
119:     MatGetRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
120:     if (flg) {
121:       MatSeqAIJGetArray(AA,&array);
122:       for (i=0; i<ia[nrows]; i++) array[i] = one;
123:       MatSeqAIJRestoreArray(AA,&array);
124:     }
125:     MatRestoreRowIJ(AA,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
126:     MatGetRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
127:     if (flg) {
128:       MatSeqAIJGetArray(AB,&array);
129:       for (i=0; i<ia[nrows]; i++) array[i] = one;
130:       MatSeqAIJRestoreArray(AB,&array);
131:     }
132:     MatRestoreRowIJ(AB,0,PETSC_FALSE,PETSC_FALSE,&nrows,&ia,&ja,&flg);
133:   }
134:   /* MatView(A, PETSC_VIEWER_STDOUT_WORLD); */

136:   /* Set up distributed array for coarse grid */
137:   if (!Test_3D) {
138:     DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,Npx,Npy,1,1,NULL,NULL,&user.coarse.da);
139:   } else {
140:     DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.coarse.mx,user.coarse.my,user.coarse.mz,Npx,Npy,Npz,
141:                         1,1,NULL,NULL,NULL,&user.coarse.da);
142:   }
143:   DMSetFromOptions(user.coarse.da);
144:   DMSetUp(user.coarse.da);

146:   /* Create interpolation between the fine and coarse grids */
147:   DMCreateInterpolation(user.coarse.da,user.fine.da,&P,NULL);
148:   MatGetLocalSize(P,&m,&n);
149:   MatGetSize(P,&M,&N);
150:   /* if (!rank) printf("P %d, %d\n",M,N); */

152:   /* Create vectors v1 and v2 that are compatible with A */
153:   VecCreate(PETSC_COMM_WORLD,&v1);
154:   MatGetLocalSize(A,&m,NULL);
155:   VecSetSizes(v1,m,PETSC_DECIDE);
156:   VecSetFromOptions(v1);
157:   VecDuplicate(v1,&v2);
158:   PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
159:   PetscRandomSetFromOptions(rdm);

161:   /* Test MatMatMult(): C = A*P */
162:   /*----------------------------*/
163:   if (Test_MatMatMult) {
164:     MatDuplicate(A,MAT_COPY_VALUES,&A_tmp);
165:     MatMatMult(A_tmp,P,MAT_INITIAL_MATRIX,fill,&C);

167:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
168:     alpha=1.0;
169:     for (i=0; i<2; i++) {
170:       alpha -= 0.1;
171:       MatScale(A_tmp,alpha);
172:       MatMatMult(A_tmp,P,MAT_REUSE_MATRIX,fill,&C);
173:     }
174:     /* Free intermediate data structures created for reuse of C=Pt*A*P */
175:     MatFreeIntermediateDataStructures(C);

177:     /* Test MatDuplicate()        */
178:     /*----------------------------*/
179:     MatDuplicate(C,MAT_COPY_VALUES,&C1);
180:     MatDuplicate(C1,MAT_COPY_VALUES,&C2);
181:     MatDestroy(&C1);
182:     MatDestroy(&C2);

184:     /* Create vector x that is compatible with P */
185:     VecCreate(PETSC_COMM_WORLD,&x);
186:     MatGetLocalSize(P,NULL,&n);
187:     VecSetSizes(x,n,PETSC_DECIDE);
188:     VecSetFromOptions(x);

190:     norm = 0.0;
191:     for (i=0; i<10; i++) {
192:       VecSetRandom(x,rdm);
193:       MatMult(P,x,v1);
194:       MatMult(A_tmp,v1,v2); /* v2 = A*P*x */
195:       MatMult(C,x,v1);  /* v1 = C*x   */
196:       VecAXPY(v1,none,v2);
197:       VecNorm(v1,NORM_1,&norm_tmp);
198:       VecNorm(v2,NORM_1,&norm_tmp1);
199:       norm_tmp /= norm_tmp1;
200:       if (norm_tmp > norm) norm = norm_tmp;
201:     }
202:     if (norm >= tol && !rank) {
203:       PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult(), |v1 - v2|/|v2|: %g\n",(double)norm);
204:     }

206:     VecDestroy(&x);
207:     MatDestroy(&C);
208:     MatDestroy(&A_tmp);
209:   }

211:   /* Test P^T * A * P - MatPtAP() */
212:   /*------------------------------*/
213:   if (Test_MatPtAP) {
214:     MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C);
215:     MatGetLocalSize(C,&m,&n);

217:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
218:     alpha=1.0;
219:     for (i=0; i<1; i++) {
220:       alpha -= 0.1;
221:       MatScale(A,alpha);
222:       MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C);
223:     }

225:     /* Free intermediate data structures created for reuse of C=Pt*A*P */
226:     MatFreeIntermediateDataStructures(C);

228:     /* Test MatDuplicate()        */
229:     /*----------------------------*/
230:     MatDuplicate(C,MAT_COPY_VALUES,&C1);
231:     MatDuplicate(C1,MAT_COPY_VALUES,&C2);
232:     MatDestroy(&C1);
233:     MatDestroy(&C2);

235:     /* Create vector x that is compatible with P */
236:     VecCreate(PETSC_COMM_WORLD,&x);
237:     MatGetLocalSize(P,&m,&n);
238:     VecSetSizes(x,n,PETSC_DECIDE);
239:     VecSetFromOptions(x);

241:     VecCreate(PETSC_COMM_WORLD,&v3);
242:     VecSetSizes(v3,n,PETSC_DECIDE);
243:     VecSetFromOptions(v3);
244:     VecDuplicate(v3,&v4);

246:     norm = 0.0;
247:     for (i=0; i<10; i++) {
248:       VecSetRandom(x,rdm);
249:       MatMult(P,x,v1);
250:       MatMult(A,v1,v2);  /* v2 = A*P*x */

252:       MatMultTranspose(P,v2,v3); /* v3 = Pt*A*P*x */
253:       MatMult(C,x,v4);           /* v3 = C*x   */
254:       VecAXPY(v4,none,v3);
255:       VecNorm(v4,NORM_1,&norm_tmp);
256:       VecNorm(v3,NORM_1,&norm_tmp1);

258:       norm_tmp /= norm_tmp1;
259:       if (norm_tmp > norm) norm = norm_tmp;
260:     }
261:     if (norm >= tol && !rank) {
262:       PetscPrintf(PETSC_COMM_SELF,"Error: MatPtAP(), |v3 - v4|/|v3|: %g\n",(double)norm);
263:     }
264:     MatDestroy(&C);
265:     VecDestroy(&v3);
266:     VecDestroy(&v4);
267:     VecDestroy(&x);
268:   }

270:   /* Clean up */
271:   MatDestroy(&A);
272:   PetscRandomDestroy(&rdm);
273:   VecDestroy(&v1);
274:   VecDestroy(&v2);
275:   DMDestroy(&user.fine.da);
276:   DMDestroy(&user.coarse.da);
277:   MatDestroy(&P);
278:   PetscFinalize();
279:   return ierr;
280: }

282: /*TEST

284:    test:
285:       args: -Mx 10 -My 5 -Mz 10
286:       output_file: output/ex96_1.out

288:    test:
289:       suffix: nonscalable
290:       nsize: 3
291:       args: -Mx 10 -My 5 -Mz 10
292:       output_file: output/ex96_1.out

294:    test:
295:       suffix: scalable
296:       nsize: 3
297:       args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable
298:       output_file: output/ex96_1.out

300:    test:
301:      suffix: seq_scalable
302:      nsize: 3
303:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via scalable -inner_offdiag_matmatmult_via scalable
304:      output_file: output/ex96_1.out

306:    test:
307:      suffix: seq_sorted
308:      nsize: 3
309:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via sorted -inner_offdiag_matmatmult_via sorted
310:      output_file: output/ex96_1.out

312:    test:
313:      suffix: seq_scalable_fast
314:      nsize: 3
315:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via scalable_fast -inner_offdiag_matmatmult_via scalable_fast
316:      output_file: output/ex96_1.out

318:    test:
319:      suffix: seq_heap
320:      nsize: 3
321:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via heap -inner_offdiag_matmatmult_via heap
322:      output_file: output/ex96_1.out

324:    test:
325:      suffix: seq_btheap
326:      nsize: 3
327:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via btheap -inner_offdiag_matmatmult_via btheap
328:      output_file: output/ex96_1.out

330:    test:
331:      suffix: seq_llcondensed
332:      nsize: 3
333:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via llcondensed -inner_offdiag_matmatmult_via llcondensed
334:      output_file: output/ex96_1.out

336:    test:
337:      suffix: seq_rowmerge
338:      nsize: 3
339:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_matmatmult_via rowmerge -inner_offdiag_matmatmult_via rowmerge
340:      output_file: output/ex96_1.out

342: TEST*/