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";
8
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 */
13
14 #include <petscdm.h>
15 #include <petscdmda.h>
16
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;
31
32 #define COARSE_LEVEL 0
33 #define FINE_LEVEL 1
34
35 /*
36 Mm_ratio - ration of grid lines between fine and coarse grids.
37 */
main(int argc,char ** argv)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;
53
54 PetscFunctionBeginUser;
55 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
56 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
57
58 user.ratio = 2;
59 user.coarse.mx = 20;
60 user.coarse.my = 20;
61 user.coarse.mz = 20;
62
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));
67
68 if (user.coarse.mz) Test_3D = PETSC_TRUE;
69
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));
75
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));
84
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));
87
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));
94
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));
100
101 /*------------------------------------------------------------*/
102
103 PetscCall(MatGetLocalSize(A, &m, &n));
104 PetscCall(MatGetSize(A, &M, &N));
105 /* if (rank == 0) printf("A %d, %d\n",M,N); */
106
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)); */
135
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); */
141
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));
150
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));
156
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));
166
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));
173
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));
179
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));
193
194 PetscCall(VecDestroy(&x));
195 PetscCall(MatDestroy(&C));
196 PetscCall(MatDestroy(&A_tmp));
197 }
198
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));
204
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 }
212
213 /* Free intermediate data structures created for reuse of C=Pt*A*P */
214 PetscCall(MatProductClear(C));
215
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));
222
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));
228
229 PetscCall(VecCreate(PETSC_COMM_WORLD, &v3));
230 PetscCall(VecSetSizes(v3, n, PETSC_DECIDE));
231 PetscCall(VecSetFromOptions(v3));
232 PetscCall(VecDuplicate(v3, &v4));
233
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 */
239
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));
245
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 }
255
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 }
267
268 /*TEST
269
270 test:
271 args: -Mx 10 -My 5 -Mz 10
272 output_file: output/empty.out
273
274 test:
275 suffix: nonscalable
276 nsize: 3
277 args: -Mx 10 -My 5 -Mz 10
278 output_file: output/empty.out
279
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/empty.out
285
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/empty.out
291
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/empty.out
297
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/empty.out
303
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/empty.out
309
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/empty.out
315
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/empty.out
321
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/empty.out
327
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/empty.out
333
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/empty.out
339
340 TEST*/
341