Actual source code: ex125.c
1: static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
2: Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 \n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A, RHS = NULL, RHS1 = NULL, C, F, X;
9: Vec u, x, b;
10: PetscMPIInt size;
11: PetscInt m, n, nfact, nsolve, nrhs, ipack = 0;
12: PetscReal norm, tol = 1.e-10;
13: IS perm, iperm;
14: MatFactorInfo info;
15: PetscRandom rand;
16: PetscBool flg, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE;
17: PetscBool chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE;
18: #if defined(PETSC_HAVE_MUMPS)
19: PetscBool test_mumps_opts = PETSC_FALSE;
20: #endif
21: PetscViewer fd; /* viewer */
22: char file[PETSC_MAX_PATH_LEN]; /* input file name */
25: PetscInitialize(&argc, &args, (char *)0, help);
26: MPI_Comm_size(PETSC_COMM_WORLD, &size);
28: /* Determine file from which we read the matrix A */
29: PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg);
30: if (flg) { /* Load matrix A */
31: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
32: MatCreate(PETSC_COMM_WORLD, &A);
33: MatSetFromOptions(A);
34: MatLoad(A, fd);
35: PetscViewerDestroy(&fd);
36: } else {
37: n = 13;
38: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
39: MatCreate(PETSC_COMM_WORLD, &A);
40: MatSetType(A, MATAIJ);
41: MatSetFromOptions(A);
42: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
43: MatSetUp(A);
44: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
45: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
46: MatShift(A, 1.0);
47: }
48: MatGetLocalSize(A, &m, &n);
51: /* if A is symmetric, set its flag -- required by MatGetInertia() */
52: MatIsSymmetric(A, 0.0, &flg);
54: MatViewFromOptions(A, NULL, "-A_view");
56: /* Create dense matrix C and X; C holds true solution with identical columns */
57: nrhs = 2;
58: PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL);
59: PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs);
60: MatCreate(PETSC_COMM_WORLD, &C);
61: MatSetOptionsPrefix(C, "rhs_");
62: MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs);
63: MatSetType(C, MATDENSE);
64: MatSetFromOptions(C);
65: MatSetUp(C);
67: PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL);
68: PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL);
69: PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL);
70: #if defined(PETSC_HAVE_MUMPS)
71: PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL);
72: #endif
74: PetscRandomCreate(PETSC_COMM_WORLD, &rand);
75: PetscRandomSetFromOptions(rand);
76: MatSetRandom(C, rand);
77: MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X);
79: /* Create vectors */
80: MatCreateVecs(A, &x, &b);
81: VecDuplicate(x, &u); /* save the true solution */
83: /* Test Factorization */
84: MatGetOrdering(A, MATORDERINGND, &perm, &iperm);
86: PetscOptionsGetInt(NULL, NULL, "-mat_solver_type", &ipack, NULL);
87: switch (ipack) {
88: #if defined(PETSC_HAVE_SUPERLU)
89: case 0:
91: PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n");
92: MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F);
93: matsolvexx = PETSC_TRUE;
94: break;
95: #endif
96: #if defined(PETSC_HAVE_SUPERLU_DIST)
97: case 1:
99: PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n");
100: MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F);
101: matsolvexx = PETSC_TRUE;
102: break;
103: #endif
104: #if defined(PETSC_HAVE_MUMPS)
105: case 2:
106: if (chol) {
107: PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n");
108: MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F);
109: } else {
110: PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n");
111: MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F);
112: }
113: matsolvexx = PETSC_TRUE;
114: if (test_mumps_opts) {
115: /* test mumps options */
116: PetscInt icntl;
117: PetscReal cntl;
119: icntl = 2; /* sequential matrix ordering */
120: MatMumpsSetIcntl(F, 7, icntl);
122: cntl = 1.e-6; /* threshold for row pivot detection */
123: MatMumpsSetIcntl(F, 24, 1);
124: MatMumpsSetCntl(F, 3, cntl);
125: }
126: break;
127: #endif
128: #if defined(PETSC_HAVE_MKL_PARDISO)
129: case 3:
130: if (chol) {
131: PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n");
132: MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F);
133: } else {
134: PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n");
135: MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F);
136: }
137: break;
138: #endif
139: #if defined(PETSC_HAVE_CUDA)
140: case 4:
141: if (chol) {
142: PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n");
143: MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F);
144: } else {
145: PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n");
146: MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F);
147: }
148: break;
149: #endif
150: default:
151: if (chol) {
152: PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n");
153: MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F);
154: } else {
155: PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n");
156: MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F);
157: }
158: matsolvexx = PETSC_TRUE;
159: }
161: MatFactorInfoInitialize(&info);
162: info.fill = 5.0;
163: info.shifttype = (PetscReal)MAT_SHIFT_NONE;
164: if (chol) {
165: MatCholeskyFactorSymbolic(F, A, perm, &info);
166: } else {
167: MatLUFactorSymbolic(F, A, perm, iperm, &info);
168: }
170: for (nfact = 0; nfact < 2; nfact++) {
171: if (chol) {
172: PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact);
173: MatCholeskyFactorNumeric(F, A, &info);
174: } else {
175: PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact);
176: MatLUFactorNumeric(F, A, &info);
177: }
178: if (view) {
179: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
180: MatView(F, PETSC_VIEWER_STDOUT_WORLD);
181: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
182: view = PETSC_FALSE;
183: }
185: #if defined(PETSC_HAVE_SUPERLU_DIST)
186: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
187: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
188: PetscInt M;
189: PetscScalar *diag;
190: #if !defined(PETSC_USE_COMPLEX)
191: PetscInt nneg, nzero, npos;
192: #endif
194: MatGetSize(F, &M, NULL);
195: PetscMalloc1(M, &diag);
196: MatSuperluDistGetDiagU(F, diag);
197: PetscFree(diag);
199: #if !defined(PETSC_USE_COMPLEX)
200: /* Test MatGetInertia() */
201: MatGetInertia(F, &nneg, &nzero, &npos);
202: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos);
203: #endif
204: }
205: #endif
207: #if defined(PETSC_HAVE_MUMPS)
208: /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
209: if (ipack == 2) {
210: if (chol) {
211: MatCholeskyFactorSymbolic(F, A, perm, &info);
212: MatCholeskyFactorNumeric(F, A, &info);
213: } else {
214: MatLUFactorSymbolic(F, A, perm, iperm, &info);
215: MatLUFactorNumeric(F, A, &info);
216: }
217: }
218: #endif
220: /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
221: if (testMatMatSolve) {
222: if (!nfact) {
223: MatMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS);
224: } else {
225: MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS);
226: }
227: for (nsolve = 0; nsolve < 2; nsolve++) {
228: PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve);
229: MatMatSolve(F, RHS, X);
231: /* Check the error */
232: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
233: MatNorm(X, NORM_FROBENIUS, &norm);
234: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
235: }
237: if (matsolvexx) {
238: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
239: MatCopy(RHS, X, SAME_NONZERO_PATTERN);
240: MatMatSolve(F, X, X);
241: /* Check the error */
242: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
243: MatNorm(X, NORM_FROBENIUS, &norm);
244: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm);
245: }
247: if (ipack == 2 && size == 1) {
248: Mat spRHS, spRHST, RHST;
250: MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST);
251: MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST);
252: MatCreateTranspose(spRHST, &spRHS);
253: for (nsolve = 0; nsolve < 2; nsolve++) {
254: PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve);
255: MatMatSolve(F, spRHS, X);
257: /* Check the error */
258: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
259: MatNorm(X, NORM_FROBENIUS, &norm);
260: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
261: }
262: MatDestroy(&spRHST);
263: MatDestroy(&spRHS);
264: MatDestroy(&RHST);
265: }
266: }
268: /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
269: if (testMatMatSolveTranspose) {
270: if (!nfact) {
271: MatTransposeMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS1);
272: } else {
273: MatTransposeMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS1);
274: }
276: for (nsolve = 0; nsolve < 2; nsolve++) {
277: MatMatSolveTranspose(F, RHS1, X);
279: /* Check the error */
280: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
281: MatNorm(X, NORM_FROBENIUS, &norm);
282: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
283: }
285: if (ipack == 2 && size == 1) {
286: Mat spRHS, spRHST, RHST;
288: MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST);
289: MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST);
290: MatCreateTranspose(spRHST, &spRHS);
291: for (nsolve = 0; nsolve < 2; nsolve++) {
292: MatMatSolveTranspose(F, spRHS, X);
294: /* Check the error */
295: MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN);
296: MatNorm(X, NORM_FROBENIUS, &norm);
297: if (norm > tol) PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve);
298: }
299: MatDestroy(&spRHST);
300: MatDestroy(&spRHS);
301: MatDestroy(&RHST);
302: }
303: }
305: /* Test MatSolve() */
306: if (testMatSolve) {
307: for (nsolve = 0; nsolve < 2; nsolve++) {
308: VecSetRandom(x, rand);
309: VecCopy(x, u);
310: MatMult(A, x, b);
312: PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve);
313: MatSolve(F, b, x);
315: /* Check the error */
316: VecAXPY(u, -1.0, x); /* u <- (-1.0)x + u */
317: VecNorm(u, NORM_2, &norm);
318: if (norm > tol) {
319: PetscReal resi;
320: MatMult(A, x, u); /* u = A*x */
321: VecAXPY(u, -1.0, b); /* u <- (-1.0)b + u */
322: VecNorm(u, NORM_2, &resi);
323: PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact);
324: }
325: }
326: }
327: }
329: /* Free data structures */
330: MatDestroy(&A);
331: MatDestroy(&C);
332: MatDestroy(&F);
333: MatDestroy(&X);
334: MatDestroy(&RHS);
335: MatDestroy(&RHS1);
337: PetscRandomDestroy(&rand);
338: ISDestroy(&perm);
339: ISDestroy(&iperm);
340: VecDestroy(&x);
341: VecDestroy(&b);
342: VecDestroy(&u);
343: PetscFinalize();
344: return 0;
345: }
347: /*TEST
349: test:
350: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
351: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 10
352: output_file: output/ex125.out
354: test:
355: suffix: 2
356: args: -mat_solver_type 10
357: output_file: output/ex125.out
359: test:
360: suffix: mkl_pardiso
361: requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
362: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 3
364: test:
365: suffix: mkl_pardiso_2
366: requires: mkl_pardiso
367: args: -mat_solver_type 3
368: output_file: output/ex125_mkl_pardiso.out
370: test:
371: suffix: mumps
372: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
373: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
374: output_file: output/ex125_mumps_seq.out
376: test:
377: suffix: mumps_2
378: nsize: 3
379: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
380: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 2
381: output_file: output/ex125_mumps_par.out
383: test:
384: suffix: mumps_3
385: requires: mumps
386: args: -mat_solver_type 2
387: output_file: output/ex125_mumps_seq.out
389: test:
390: suffix: mumps_4
391: nsize: 3
392: requires: mumps
393: args: -mat_solver_type 2
394: output_file: output/ex125_mumps_par.out
396: test:
397: suffix: mumps_5
398: nsize: 3
399: requires: mumps
400: args: -mat_solver_type 2 -cholesky
401: output_file: output/ex125_mumps_par_cholesky.out
403: test:
404: suffix: superlu_dist
405: nsize: {{1 3}}
406: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
407: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
409: test:
410: suffix: superlu_dist_2
411: nsize: {{1 3}}
412: requires: superlu_dist !complex
413: args: -n 36 -mat_solver_type 1 -mat_superlu_dist_rowperm NOROWPERM
414: output_file: output/ex125_superlu_dist.out
416: test:
417: suffix: superlu_dist_complex
418: nsize: 3
419: requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
420: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type 1
421: output_file: output/ex125_superlu_dist_complex.out
423: test:
424: suffix: superlu_dist_complex_2
425: nsize: 3
426: requires: superlu_dist complex
427: args: -mat_solver_type 1
428: output_file: output/ex125_superlu_dist_complex.out
430: test:
431: suffix: cusparse
432: requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
433: #todo: fix the bug with cholesky
434: #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0 1}separate output}
435: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type 4 -cholesky {{0}separate output}
437: test:
438: suffix: cusparse_2
439: requires: cuda
440: args: -mat_type aijcusparse -mat_solver_type 4 -cholesky {{0 1}separate output}
442: TEST*/