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*/