Actual source code: ex96.c
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/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: {
41: AppCtx user;
42: PetscInt Npx = PETSC_DECIDE, Npy = PETSC_DECIDE, Npz = PETSC_DECIDE;
43: PetscMPIInt size, rank;
44: PetscInt m, n, M, N, i, nrows;
45: PetscScalar one = 1.0;
46: PetscReal fill = 2.0;
47: Mat A, A_tmp, P, C, C1, C2;
48: PetscScalar *array, none = -1.0, alpha;
49: Vec x, v1, v2, v3, v4;
50: PetscReal norm, norm_tmp, norm_tmp1, tol = 100. * PETSC_MACHINE_EPSILON;
51: PetscRandom rdm;
52: PetscBool Test_MatMatMult = PETSC_TRUE, Test_MatPtAP = PETSC_TRUE, Test_3D = PETSC_TRUE, flg;
53: const PetscInt *ia, *ja;
56: PetscInitialize(&argc, &argv, NULL, help);
57: PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL);
59: user.ratio = 2;
60: user.coarse.mx = 20;
61: user.coarse.my = 20;
62: user.coarse.mz = 20;
64: PetscOptionsGetInt(NULL, NULL, "-Mx", &user.coarse.mx, NULL);
65: PetscOptionsGetInt(NULL, NULL, "-My", &user.coarse.my, NULL);
66: PetscOptionsGetInt(NULL, NULL, "-Mz", &user.coarse.mz, NULL);
67: PetscOptionsGetInt(NULL, NULL, "-ratio", &user.ratio, NULL);
69: if (user.coarse.mz) Test_3D = PETSC_TRUE;
71: MPI_Comm_size(PETSC_COMM_WORLD, &size);
72: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
73: PetscOptionsGetInt(NULL, NULL, "-Npx", &Npx, NULL);
74: PetscOptionsGetInt(NULL, NULL, "-Npy", &Npy, NULL);
75: PetscOptionsGetInt(NULL, NULL, "-Npz", &Npz, NULL);
77: /* Set up distributed array for fine grid */
78: if (!Test_3D) {
79: 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);
80: } else {
81: 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);
82: }
83: DMSetFromOptions(user.coarse.da);
84: DMSetUp(user.coarse.da);
86: /* This makes sure the coarse DMDA has the same partition as the fine DMDA */
87: DMRefine(user.coarse.da, PetscObjectComm((PetscObject)user.coarse.da), &user.fine.da);
89: /* Test DMCreateMatrix() */
90: /*------------------------------------------------------------*/
91: DMSetMatType(user.fine.da, MATAIJ);
92: DMCreateMatrix(user.fine.da, &A);
93: DMSetMatType(user.fine.da, MATBAIJ);
94: DMCreateMatrix(user.fine.da, &C);
96: MatConvert(C, MATAIJ, MAT_INITIAL_MATRIX, &A_tmp); /* not work for mpisbaij matrix! */
97: MatEqual(A, A_tmp, &flg);
99: MatDestroy(&C);
100: MatDestroy(&A_tmp);
102: /*------------------------------------------------------------*/
104: MatGetLocalSize(A, &m, &n);
105: MatGetSize(A, &M, &N);
106: /* if (rank == 0) printf("A %d, %d\n",M,N); */
108: /* set val=one to A */
109: if (size == 1) {
110: MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
111: if (flg) {
112: MatSeqAIJGetArray(A, &array);
113: for (i = 0; i < ia[nrows]; i++) array[i] = one;
114: MatSeqAIJRestoreArray(A, &array);
115: }
116: MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
117: } else {
118: Mat AA, AB;
119: MatMPIAIJGetSeqAIJ(A, &AA, &AB, NULL);
120: MatGetRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
121: if (flg) {
122: MatSeqAIJGetArray(AA, &array);
123: for (i = 0; i < ia[nrows]; i++) array[i] = one;
124: MatSeqAIJRestoreArray(AA, &array);
125: }
126: MatRestoreRowIJ(AA, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
127: MatGetRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
128: if (flg) {
129: MatSeqAIJGetArray(AB, &array);
130: for (i = 0; i < ia[nrows]; i++) array[i] = one;
131: MatSeqAIJRestoreArray(AB, &array);
132: }
133: MatRestoreRowIJ(AB, 0, PETSC_FALSE, PETSC_FALSE, &nrows, &ia, &ja, &flg);
134: }
135: /* MatView(A, PETSC_VIEWER_STDOUT_WORLD); */
137: /* Create interpolation between the fine and coarse grids */
138: DMCreateInterpolation(user.coarse.da, user.fine.da, &P, NULL);
139: MatGetLocalSize(P, &m, &n);
140: MatGetSize(P, &M, &N);
141: /* if (rank == 0) printf("P %d, %d\n",M,N); */
143: /* Create vectors v1 and v2 that are compatible with A */
144: VecCreate(PETSC_COMM_WORLD, &v1);
145: MatGetLocalSize(A, &m, NULL);
146: VecSetSizes(v1, m, PETSC_DECIDE);
147: VecSetFromOptions(v1);
148: VecDuplicate(v1, &v2);
149: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
150: PetscRandomSetFromOptions(rdm);
152: /* Test MatMatMult(): C = A*P */
153: /*----------------------------*/
154: if (Test_MatMatMult) {
155: MatDuplicate(A, MAT_COPY_VALUES, &A_tmp);
156: MatMatMult(A_tmp, P, MAT_INITIAL_MATRIX, fill, &C);
158: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
159: alpha = 1.0;
160: for (i = 0; i < 2; i++) {
161: alpha -= 0.1;
162: MatScale(A_tmp, alpha);
163: MatMatMult(A_tmp, P, MAT_REUSE_MATRIX, fill, &C);
164: }
165: /* Free intermediate data structures created for reuse of C=Pt*A*P */
166: MatProductClear(C);
168: /* Test MatDuplicate() */
169: /*----------------------------*/
170: MatDuplicate(C, MAT_COPY_VALUES, &C1);
171: MatDuplicate(C1, MAT_COPY_VALUES, &C2);
172: MatDestroy(&C1);
173: MatDestroy(&C2);
175: /* Create vector x that is compatible with P */
176: VecCreate(PETSC_COMM_WORLD, &x);
177: MatGetLocalSize(P, NULL, &n);
178: VecSetSizes(x, n, PETSC_DECIDE);
179: VecSetFromOptions(x);
181: norm = 0.0;
182: for (i = 0; i < 10; i++) {
183: VecSetRandom(x, rdm);
184: MatMult(P, x, v1);
185: MatMult(A_tmp, v1, v2); /* v2 = A*P*x */
186: MatMult(C, x, v1); /* v1 = C*x */
187: VecAXPY(v1, none, v2);
188: VecNorm(v1, NORM_1, &norm_tmp);
189: VecNorm(v2, NORM_1, &norm_tmp1);
190: norm_tmp /= norm_tmp1;
191: if (norm_tmp > norm) norm = norm_tmp;
192: }
193: if (norm >= tol && rank == 0) PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult(), |v1 - v2|/|v2|: %g\n", (double)norm);
195: VecDestroy(&x);
196: MatDestroy(&C);
197: MatDestroy(&A_tmp);
198: }
200: /* Test P^T * A * P - MatPtAP() */
201: /*------------------------------*/
202: if (Test_MatPtAP) {
203: MatPtAP(A, P, MAT_INITIAL_MATRIX, fill, &C);
204: MatGetLocalSize(C, &m, &n);
206: /* Test MAT_REUSE_MATRIX - reuse symbolic C */
207: alpha = 1.0;
208: for (i = 0; i < 1; i++) {
209: alpha -= 0.1;
210: MatScale(A, alpha);
211: MatPtAP(A, P, MAT_REUSE_MATRIX, fill, &C);
212: }
214: /* Free intermediate data structures created for reuse of C=Pt*A*P */
215: MatProductClear(C);
217: /* Test MatDuplicate() */
218: /*----------------------------*/
219: MatDuplicate(C, MAT_COPY_VALUES, &C1);
220: MatDuplicate(C1, MAT_COPY_VALUES, &C2);
221: MatDestroy(&C1);
222: MatDestroy(&C2);
224: /* Create vector x that is compatible with P */
225: VecCreate(PETSC_COMM_WORLD, &x);
226: MatGetLocalSize(P, &m, &n);
227: VecSetSizes(x, n, PETSC_DECIDE);
228: VecSetFromOptions(x);
230: VecCreate(PETSC_COMM_WORLD, &v3);
231: VecSetSizes(v3, n, PETSC_DECIDE);
232: VecSetFromOptions(v3);
233: VecDuplicate(v3, &v4);
235: norm = 0.0;
236: for (i = 0; i < 10; i++) {
237: VecSetRandom(x, rdm);
238: MatMult(P, x, v1);
239: MatMult(A, v1, v2); /* v2 = A*P*x */
241: MatMultTranspose(P, v2, v3); /* v3 = Pt*A*P*x */
242: MatMult(C, x, v4); /* v3 = C*x */
243: VecAXPY(v4, none, v3);
244: VecNorm(v4, NORM_1, &norm_tmp);
245: VecNorm(v3, NORM_1, &norm_tmp1);
247: norm_tmp /= norm_tmp1;
248: if (norm_tmp > norm) norm = norm_tmp;
249: }
250: if (norm >= tol && rank == 0) PetscPrintf(PETSC_COMM_SELF, "Error: MatPtAP(), |v3 - v4|/|v3|: %g\n", (double)norm);
251: MatDestroy(&C);
252: VecDestroy(&v3);
253: VecDestroy(&v4);
254: VecDestroy(&x);
255: }
257: /* Clean up */
258: MatDestroy(&A);
259: PetscRandomDestroy(&rdm);
260: VecDestroy(&v1);
261: VecDestroy(&v2);
262: DMDestroy(&user.fine.da);
263: DMDestroy(&user.coarse.da);
264: MatDestroy(&P);
265: PetscFinalize();
266: return 0;
267: }
269: /*TEST
271: test:
272: args: -Mx 10 -My 5 -Mz 10
273: output_file: output/ex96_1.out
275: test:
276: suffix: nonscalable
277: nsize: 3
278: args: -Mx 10 -My 5 -Mz 10
279: output_file: output/ex96_1.out
281: test:
282: suffix: scalable
283: nsize: 3
284: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via scalable
285: output_file: output/ex96_1.out
287: test:
288: suffix: seq_scalable
289: nsize: 3
290: 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
291: output_file: output/ex96_1.out
293: test:
294: suffix: seq_sorted
295: nsize: 3
296: 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
297: output_file: output/ex96_1.out
299: test:
300: suffix: seq_scalable_fast
301: nsize: 3
302: 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
303: output_file: output/ex96_1.out
305: test:
306: suffix: seq_heap
307: nsize: 3
308: 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
309: output_file: output/ex96_1.out
311: test:
312: suffix: seq_btheap
313: nsize: 3
314: 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
315: output_file: output/ex96_1.out
317: test:
318: suffix: seq_llcondensed
319: nsize: 3
320: 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
321: output_file: output/ex96_1.out
323: test:
324: suffix: seq_rowmerge
325: nsize: 3
326: 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
327: output_file: output/ex96_1.out
329: test:
330: suffix: allatonce
331: nsize: 3
332: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce
333: output_file: output/ex96_1.out
335: test:
336: suffix: allatonce_merged
337: nsize: 3
338: args: -Mx 10 -My 5 -Mz 10 -matmatmult_via scalable -matptap_via allatonce_merged
339: output_file: output/ex96_1.out
341: TEST*/