Actual source code: ex96.c

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

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

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

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

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

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

 54:   PetscFunctionBeginUser;
 55:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 56:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));

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

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

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

 70:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 71:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 72:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Npx", &Npx, NULL));
 73:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Npy", &Npy, NULL));
 74:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Npz", &Npz, NULL));

 76:   /* Set up distributed array for fine grid */
 77:   if (!Test_3D) {
 78:     PetscCall(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));
 79:   } else {
 80:     PetscCall(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, 1, 1, NULL, NULL, NULL, &user.coarse.da));
 81:   }
 82:   PetscCall(DMSetFromOptions(user.coarse.da));
 83:   PetscCall(DMSetUp(user.coarse.da));

 85:   /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
 86:   PetscCall(DMRefine(user.coarse.da, PetscObjectComm((PetscObject)user.coarse.da), &user.fine.da));

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

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

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

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

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

136:   /* Create interpolation between the fine and coarse grids */
137:   PetscCall(DMCreateInterpolation(user.coarse.da, user.fine.da, &P, NULL));
138:   PetscCall(MatGetLocalSize(P, &m, &n));
139:   PetscCall(MatGetSize(P, &M, &N));
140:   /* if (rank == 0) printf("P %d, %d\n",M,N); */

142:   /* Create vectors v1 and v2 that are compatible with A */
143:   PetscCall(VecCreate(PETSC_COMM_WORLD, &v1));
144:   PetscCall(MatGetLocalSize(A, &m, NULL));
145:   PetscCall(VecSetSizes(v1, m, PETSC_DECIDE));
146:   PetscCall(VecSetFromOptions(v1));
147:   PetscCall(VecDuplicate(v1, &v2));
148:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
149:   PetscCall(PetscRandomSetFromOptions(rdm));

151:   /* Test MatMatMult(): C = A*P */
152:   /*----------------------------*/
153:   if (Test_MatMatMult) {
154:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_tmp));
155:     PetscCall(MatMatMult(A_tmp, P, MAT_INITIAL_MATRIX, fill, &C));

157:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
158:     alpha = 1.0;
159:     for (i = 0; i < 2; i++) {
160:       alpha -= 0.1;
161:       PetscCall(MatScale(A_tmp, alpha));
162:       PetscCall(MatMatMult(A_tmp, P, MAT_REUSE_MATRIX, fill, &C));
163:     }
164:     /* Free intermediate data structures created for reuse of C=Pt*A*P */
165:     PetscCall(MatProductClear(C));

167:     /* Test MatDuplicate()        */
168:     /*----------------------------*/
169:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
170:     PetscCall(MatDuplicate(C1, MAT_COPY_VALUES, &C2));
171:     PetscCall(MatDestroy(&C1));
172:     PetscCall(MatDestroy(&C2));

174:     /* Create vector x that is compatible with P */
175:     PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
176:     PetscCall(MatGetLocalSize(P, NULL, &n));
177:     PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
178:     PetscCall(VecSetFromOptions(x));

180:     norm = 0.0;
181:     for (i = 0; i < 10; i++) {
182:       PetscCall(VecSetRandom(x, rdm));
183:       PetscCall(MatMult(P, x, v1));
184:       PetscCall(MatMult(A_tmp, v1, v2)); /* v2 = A*P*x */
185:       PetscCall(MatMult(C, x, v1));      /* v1 = C*x   */
186:       PetscCall(VecAXPY(v1, none, v2));
187:       PetscCall(VecNorm(v1, NORM_1, &norm_tmp));
188:       PetscCall(VecNorm(v2, NORM_1, &norm_tmp1));
189:       norm_tmp /= norm_tmp1;
190:       if (norm_tmp > norm) norm = norm_tmp;
191:     }
192:     if (norm >= tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult(), |v1 - v2|/|v2|: %g\n", (double)norm));

194:     PetscCall(VecDestroy(&x));
195:     PetscCall(MatDestroy(&C));
196:     PetscCall(MatDestroy(&A_tmp));
197:   }

199:   /* Test P^T * A * P - MatPtAP() */
200:   /*------------------------------*/
201:   if (Test_MatPtAP) {
202:     PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C));
203:     PetscCall(MatGetLocalSize(C, &m, &n));

205:     /* Test MAT_REUSE_MATRIX - reuse symbolic C */
206:     alpha = 1.0;
207:     for (i = 0; i < 1; i++) {
208:       alpha -= 0.1;
209:       PetscCall(MatScale(A, alpha));
210:       PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C));
211:     }

213:     /* Free intermediate data structures created for reuse of C=Pt*A*P */
214:     PetscCall(MatProductClear(C));

216:     /* Test MatDuplicate()        */
217:     /*----------------------------*/
218:     PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &C1));
219:     PetscCall(MatDuplicate(C1, MAT_COPY_VALUES, &C2));
220:     PetscCall(MatDestroy(&C1));
221:     PetscCall(MatDestroy(&C2));

223:     /* Create vector x that is compatible with P */
224:     PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
225:     PetscCall(MatGetLocalSize(P, &m, &n));
226:     PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
227:     PetscCall(VecSetFromOptions(x));

229:     PetscCall(VecCreate(PETSC_COMM_WORLD, &v3));
230:     PetscCall(VecSetSizes(v3, n, PETSC_DECIDE));
231:     PetscCall(VecSetFromOptions(v3));
232:     PetscCall(VecDuplicate(v3, &v4));

234:     norm = 0.0;
235:     for (i = 0; i < 10; i++) {
236:       PetscCall(VecSetRandom(x, rdm));
237:       PetscCall(MatMult(P, x, v1));
238:       PetscCall(MatMult(A, v1, v2)); /* v2 = A*P*x */

240:       PetscCall(MatMultTranspose(P, v2, v3)); /* v3 = Pt*A*P*x */
241:       PetscCall(MatMult(C, x, v4));           /* v3 = C*x   */
242:       PetscCall(VecAXPY(v4, none, v3));
243:       PetscCall(VecNorm(v4, NORM_1, &norm_tmp));
244:       PetscCall(VecNorm(v3, NORM_1, &norm_tmp1));

246:       norm_tmp /= norm_tmp1;
247:       if (norm_tmp > norm) norm = norm_tmp;
248:     }
249:     if (norm >= tol && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatPtAP(), |v3 - v4|/|v3|: %g\n", (double)norm));
250:     PetscCall(MatDestroy(&C));
251:     PetscCall(VecDestroy(&v3));
252:     PetscCall(VecDestroy(&v4));
253:     PetscCall(VecDestroy(&x));
254:   }

256:   /* Clean up */
257:   PetscCall(MatDestroy(&A));
258:   PetscCall(PetscRandomDestroy(&rdm));
259:   PetscCall(VecDestroy(&v1));
260:   PetscCall(VecDestroy(&v2));
261:   PetscCall(DMDestroy(&user.fine.da));
262:   PetscCall(DMDestroy(&user.coarse.da));
263:   PetscCall(MatDestroy(&P));
264:   PetscCall(PetscFinalize());
265:   return 0;
266: }

268: /*TEST

270:    test:
271:       args: -Mx 10 -My 5 -Mz 10
272:       output_file: output/ex96_1.out

274:    test:
275:       suffix: nonscalable
276:       nsize: 3
277:       args: -Mx 10 -My 5 -Mz 10
278:       output_file: output/ex96_1.out

280:    test:
281:       suffix: scalable
282:       nsize: 3
283:       args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable
284:       output_file: output/ex96_1.out

286:    test:
287:      suffix: seq_scalable
288:      nsize: 3
289:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm scalable -inner_offdiag_mat_product_algorithm scalable
290:      output_file: output/ex96_1.out

292:    test:
293:      suffix: seq_sorted
294:      nsize: 3
295:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm sorted -inner_offdiag_mat_product_algorithm sorted
296:      output_file: output/ex96_1.out

298:    test:
299:      suffix: seq_scalable_fast
300:      nsize: 3
301:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm scalable_fast -inner_offdiag_mat_product_algorithm scalable_fast
302:      output_file: output/ex96_1.out

304:    test:
305:      suffix: seq_heap
306:      nsize: 3
307:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm heap -inner_offdiag_mat_product_algorithm heap
308:      output_file: output/ex96_1.out

310:    test:
311:      suffix: seq_btheap
312:      nsize: 3
313:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm btheap -inner_offdiag_mat_product_algorithm btheap
314:      output_file: output/ex96_1.out

316:    test:
317:      suffix: seq_llcondensed
318:      nsize: 3
319:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm llcondensed -inner_offdiag_mat_product_algorithm llcondensed
320:      output_file: output/ex96_1.out

322:    test:
323:      suffix: seq_rowmerge
324:      nsize: 3
325:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable -inner_diag_mat_product_algorithm rowmerge -inner_offdiag_mat_product_algorithm rowmerge
326:      output_file: output/ex96_1.out

328:    test:
329:      suffix: allatonce
330:      nsize: 3
331:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce
332:      output_file: output/ex96_1.out

334:    test:
335:      suffix: allatonce_merged
336:      nsize: 3
337:      args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce_merged
338:      output_file: output/ex96_1.out

340: TEST*/