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 -mat_solver_type <>\n\n";
4: /*
5: -mat_solver_type:
6: superlu
7: superlu_dist
8: mumps
9: mkl_pardiso
10: cusparse
11: petsc
12: */
14: #include <petscmat.h>
16: PetscErrorCode CreateRandom(PetscInt n, PetscInt m, Mat *A)
17: {
18: PetscFunctionBeginUser;
19: PetscCall(MatCreate(PETSC_COMM_WORLD, A));
20: PetscCall(MatSetType(*A, MATAIJ));
21: PetscCall(MatSetFromOptions(*A));
22: PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, m));
23: PetscCall(MatSeqAIJSetPreallocation(*A, 5, NULL));
24: PetscCall(MatMPIAIJSetPreallocation(*A, 5, NULL, 5, NULL));
25: PetscCall(MatSetRandom(*A, NULL));
26: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
27: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: PetscErrorCode CreateIdentity(PetscInt n, Mat *A)
32: {
33: PetscFunctionBeginUser;
34: PetscCall(MatCreate(PETSC_COMM_WORLD, A));
35: PetscCall(MatSetType(*A, MATAIJ));
36: PetscCall(MatSetFromOptions(*A));
37: PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, n, n));
38: PetscCall(MatSetUp(*A));
39: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
40: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
41: PetscCall(MatShift(*A, 1.0));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: int main(int argc, char **args)
46: {
47: Mat A, Ae, RHS = NULL, RHS1 = NULL, C, F, X;
48: Vec u, x, b;
49: PetscMPIInt size;
50: PetscInt m, n, nfact, nsolve, nrhs, ipack = 5;
51: PetscReal norm, tol = 10 * PETSC_SQRT_MACHINE_EPSILON;
52: IS perm = NULL, iperm = NULL;
53: MatFactorInfo info;
54: PetscRandom rand;
55: PetscBool flg, symm, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE, testMatSolveTranspose = PETSC_TRUE, match = PETSC_FALSE;
56: PetscBool chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE;
57: #if defined(PETSC_HAVE_MUMPS)
58: PetscBool test_mumps_opts = PETSC_FALSE;
59: #endif
60: PetscViewer fd; /* viewer */
61: char file[PETSC_MAX_PATH_LEN]; /* input file name */
62: char pack[PETSC_MAX_PATH_LEN];
64: PetscFunctionBeginUser;
65: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
66: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
68: /* Determine file from which we read the matrix A */
69: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
70: if (flg) { /* Load matrix A */
71: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
72: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
73: PetscCall(MatSetFromOptions(A));
74: PetscCall(MatLoad(A, fd));
75: PetscCall(PetscViewerDestroy(&fd));
76: } else {
77: n = 13;
78: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
79: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
80: PetscCall(MatSetType(A, MATAIJ));
81: PetscCall(MatSetFromOptions(A));
82: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
83: PetscCall(MatSetUp(A));
84: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatShift(A, 1.0));
87: }
89: /* if A is symmetric, set its flag -- required by MatGetInertia() */
90: PetscCall(MatIsSymmetric(A, 0.0, &symm));
91: PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
93: PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL));
95: /* test MATNEST support */
96: flg = PETSC_FALSE;
97: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest", &flg, NULL));
98: if (flg) {
99: Mat B;
101: flg = PETSC_FALSE;
102: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest_bordered", &flg, NULL));
103: if (!flg) {
104: Mat mats[9] = {NULL, NULL, A, NULL, A, NULL, A, NULL, NULL};
106: /* Create a nested matrix representing
107: | 0 0 A |
108: | 0 A 0 |
109: | A 0 0 |
110: */
111: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, mats, &B));
112: } else {
113: Mat mats[4];
115: /* Create a nested matrix representing
116: | Id R |
117: | R^t A |
118: */
119: PetscCall(MatGetSize(A, NULL, &n));
120: m = n + 12;
121: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
122: PetscCall(CreateIdentity(m, &mats[0]));
123: PetscCall(CreateRandom(m, n, &mats[1]));
124: mats[3] = A;
126: /* use CreateTranspose/CreateHermitianTranspose or explicit matrix for debugging purposes */
127: flg = PETSC_FALSE;
128: PetscCall(PetscOptionsGetBool(NULL, NULL, "-expl", &flg, NULL));
129: if (PetscDefined(USE_COMPLEX)) {
130: if (chol) { /* Hermitian transpose not supported by MUMPS Cholesky factor */
131: if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
132: else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
133: } else {
134: if (!flg) PetscCall(MatCreateHermitianTranspose(mats[1], &mats[2]));
135: else PetscCall(MatHermitianTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
136: }
137: } else {
138: if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
139: else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
140: }
141: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &B));
142: PetscCall(MatDestroy(&mats[0]));
143: PetscCall(MatDestroy(&mats[1]));
144: PetscCall(MatDestroy(&mats[2]));
145: }
146: PetscCall(MatDestroy(&A));
147: A = B;
148: PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
150: /* not all the combinations of MatMat operations are supported by MATNEST. */
151: PetscCall(MatComputeOperator(A, MATAIJ, &Ae));
152: } else {
153: PetscCall(PetscObjectReference((PetscObject)A));
154: Ae = A;
155: }
156: PetscCall(MatGetLocalSize(A, &m, &n));
157: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
159: PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
160: PetscCall(MatViewFromOptions(Ae, NULL, "-A_view_expl"));
162: /* Create dense matrix C and X; C holds true solution with identical columns */
163: nrhs = 2;
164: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
165: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs));
166: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
167: PetscCall(MatSetOptionsPrefix(C, "rhs_"));
168: PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
169: PetscCall(MatSetType(C, MATDENSE));
170: PetscCall(MatSetFromOptions(C));
171: PetscCall(MatSetUp(C));
173: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL));
174: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL));
175: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL));
176: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL));
177: #if defined(PETSC_HAVE_MUMPS)
178: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL));
179: #endif
181: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
182: PetscCall(PetscRandomSetFromOptions(rand));
183: PetscCall(MatSetRandom(C, rand));
184: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
186: /* Create vectors */
187: PetscCall(MatCreateVecs(A, &x, &b));
188: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
190: /* Test Factorization */
191: PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
193: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
194: #if defined(PETSC_HAVE_SUPERLU)
195: PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match));
196: if (match) {
197: PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
198: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n"));
199: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
200: matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */
201: ipack = 0;
202: goto skipoptions;
203: }
204: #endif
205: #if defined(PETSC_HAVE_SUPERLU_DIST)
206: PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match));
207: if (match) {
208: PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
209: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n"));
210: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
211: matsolvexx = PETSC_TRUE;
212: if (symm) { /* A is symmetric */
213: testMatMatSolveTranspose = PETSC_TRUE;
214: testMatSolveTranspose = PETSC_TRUE;
215: } else { /* superlu_dist does not support solving A^t x = rhs yet */
216: testMatMatSolveTranspose = PETSC_FALSE;
217: testMatSolveTranspose = PETSC_FALSE;
218: }
219: ipack = 1;
220: goto skipoptions;
221: }
222: #endif
223: #if defined(PETSC_HAVE_MUMPS)
224: PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match));
225: if (match) {
226: if (chol) {
227: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n"));
228: PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F));
229: } else {
230: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n"));
231: PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
232: }
233: matsolvexx = PETSC_TRUE;
234: if (test_mumps_opts) {
235: /* test mumps options */
236: PetscInt icntl;
237: PetscReal cntl;
239: icntl = 2; /* sequential matrix ordering */
240: PetscCall(MatMumpsSetIcntl(F, 7, icntl));
242: cntl = 1.e-6; /* threshold for row pivot detection */
243: PetscCall(MatMumpsSetIcntl(F, 24, 1));
244: PetscCall(MatMumpsSetCntl(F, 3, cntl));
245: }
246: ipack = 2;
247: goto skipoptions;
248: }
249: #endif
250: #if defined(PETSC_HAVE_MKL_PARDISO)
251: PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match));
252: if (match) {
253: if (chol) {
254: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n"));
255: PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F));
256: } else {
257: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n"));
258: PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F));
259: }
260: ipack = 3;
261: goto skipoptions;
262: }
263: #endif
264: #if defined(PETSC_HAVE_CUDA)
265: PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match));
266: if (match) {
267: if (chol) {
268: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n"));
269: PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F));
270: } else {
271: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n"));
272: PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F));
273: }
274: testMatSolveTranspose = PETSC_FALSE;
275: testMatMatSolveTranspose = PETSC_FALSE;
276: ipack = 4;
277: goto skipoptions;
278: }
279: #endif
280: /* PETSc */
281: match = PETSC_TRUE;
282: if (match) {
283: if (chol) {
284: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n"));
285: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
286: } else {
287: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n"));
288: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
289: }
290: matsolvexx = PETSC_TRUE;
291: ipack = 5;
292: goto skipoptions;
293: }
295: skipoptions:
296: PetscCall(MatFactorInfoInitialize(&info));
297: info.fill = 5.0;
298: info.shifttype = (PetscReal)MAT_SHIFT_NONE;
299: if (chol) {
300: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
301: } else {
302: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
303: }
305: for (nfact = 0; nfact < 2; nfact++) {
306: if (chol) {
307: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact));
308: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
309: } else {
310: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact));
311: PetscCall(MatLUFactorNumeric(F, A, &info));
312: }
313: if (view) {
314: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
315: PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
316: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
317: view = PETSC_FALSE;
318: }
320: #if defined(PETSC_HAVE_SUPERLU_DIST)
321: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
322: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
323: PetscInt M;
324: PetscScalar *diag;
325: #if !defined(PETSC_USE_COMPLEX)
326: PetscInt nneg, nzero, npos;
327: #endif
329: PetscCall(MatGetSize(F, &M, NULL));
330: PetscCall(PetscMalloc1(M, &diag));
331: PetscCall(MatSuperluDistGetDiagU(F, diag));
332: PetscCall(PetscFree(diag));
334: #if !defined(PETSC_USE_COMPLEX)
335: /* Test MatGetInertia() */
336: if (symm) { /* A is symmetric */
337: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
338: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
339: }
340: #endif
341: }
342: #endif
344: #if defined(PETSC_HAVE_MUMPS)
345: /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
346: if (ipack == 2) {
347: if (chol) {
348: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
349: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
350: } else {
351: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
352: PetscCall(MatLUFactorNumeric(F, A, &info));
353: }
354: }
355: #endif
357: /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
358: if (testMatMatSolve) {
359: if (!nfact) {
360: PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
361: } else {
362: PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS));
363: }
364: for (nsolve = 0; nsolve < 2; nsolve++) {
365: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
366: PetscCall(MatMatSolve(F, RHS, X));
368: /* Check the error */
369: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
370: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
371: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
372: }
374: if (matsolvexx) {
375: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
376: PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
377: PetscCall(MatMatSolve(F, X, X));
378: /* Check the error */
379: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
380: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
381: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
382: }
384: if (ipack == 2 && size == 1) {
385: Mat spRHS, spRHST, RHST;
387: PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST));
388: PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
389: PetscCall(MatCreateTranspose(spRHST, &spRHS));
390: for (nsolve = 0; nsolve < 2; nsolve++) {
391: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
392: PetscCall(MatMatSolve(F, spRHS, X));
394: /* Check the error */
395: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
396: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
397: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
398: }
399: PetscCall(MatDestroy(&spRHST));
400: PetscCall(MatDestroy(&spRHS));
401: PetscCall(MatDestroy(&RHST));
402: }
403: }
405: /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
406: if (testMatMatSolveTranspose) {
407: if (!nfact) {
408: PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
409: } else {
410: PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
411: }
413: for (nsolve = 0; nsolve < 2; nsolve++) {
414: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
415: PetscCall(MatMatSolveTranspose(F, RHS1, X));
417: /* Check the error */
418: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
419: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
420: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
421: }
423: if (ipack == 2 && size == 1) {
424: Mat spRHS, spRHST, RHST;
426: PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
427: PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
428: PetscCall(MatCreateTranspose(spRHST, &spRHS));
429: for (nsolve = 0; nsolve < 2; nsolve++) {
430: PetscCall(MatMatSolveTranspose(F, spRHS, X));
432: /* Check the error */
433: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
434: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
435: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the sparse MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
436: }
437: PetscCall(MatDestroy(&spRHST));
438: PetscCall(MatDestroy(&spRHS));
439: PetscCall(MatDestroy(&RHST));
440: }
441: }
443: /* Test MatSolve() */
444: if (testMatSolve) {
445: for (nsolve = 0; nsolve < 2; nsolve++) {
446: PetscCall(VecSetRandom(x, rand));
447: PetscCall(VecCopy(x, u));
448: PetscCall(MatMult(Ae, x, b));
450: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve));
451: PetscCall(MatSolve(F, b, x));
453: /* Check the error */
454: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
455: PetscCall(VecNorm(u, NORM_2, &norm));
456: if (norm > tol) {
457: PetscReal resi;
458: PetscCall(MatMult(Ae, x, u)); /* u = A*x */
459: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
460: PetscCall(VecNorm(u, NORM_2, &resi));
461: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
462: }
463: }
464: }
466: /* Test MatSolveTranspose() */
467: if (testMatSolveTranspose) {
468: for (nsolve = 0; nsolve < 2; nsolve++) {
469: PetscCall(VecSetRandom(x, rand));
470: PetscCall(VecCopy(x, u));
471: PetscCall(MatMultTranspose(Ae, x, b));
473: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
474: PetscCall(MatSolveTranspose(F, b, x));
476: /* Check the error */
477: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
478: PetscCall(VecNorm(u, NORM_2, &norm));
479: if (norm > tol) {
480: PetscReal resi;
481: PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */
482: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
483: PetscCall(VecNorm(u, NORM_2, &resi));
484: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
485: }
486: }
487: }
488: }
490: /* Free data structures */
491: PetscCall(MatDestroy(&Ae));
492: PetscCall(MatDestroy(&A));
493: PetscCall(MatDestroy(&C));
494: PetscCall(MatDestroy(&F));
495: PetscCall(MatDestroy(&X));
496: PetscCall(MatDestroy(&RHS));
497: PetscCall(MatDestroy(&RHS1));
499: PetscCall(PetscRandomDestroy(&rand));
500: PetscCall(ISDestroy(&perm));
501: PetscCall(ISDestroy(&iperm));
502: PetscCall(VecDestroy(&x));
503: PetscCall(VecDestroy(&b));
504: PetscCall(VecDestroy(&u));
505: PetscCall(PetscFinalize());
506: return 0;
507: }
509: /*TEST
511: test:
512: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
513: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
514: output_file: output/ex125.out
516: test:
517: suffix: 2
518: args: -mat_solver_type petsc
519: output_file: output/ex125.out
521: test:
522: suffix: mkl_pardiso
523: requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
524: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso
526: test:
527: suffix: mkl_pardiso_2
528: requires: mkl_pardiso
529: args: -mat_solver_type mkl_pardiso
530: output_file: output/ex125_mkl_pardiso.out
532: test:
533: suffix: mumps
534: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
535: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
536: output_file: output/ex125_mumps_seq.out
538: test:
539: suffix: mumps_nest
540: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
541: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
542: output_file: output/ex125_mumps_seq.out
544: test:
545: suffix: mumps_2
546: nsize: 3
547: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
548: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
549: output_file: output/ex125_mumps_par.out
551: test:
552: suffix: mumps_2_nest
553: nsize: 3
554: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
555: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
556: output_file: output/ex125_mumps_par.out
558: test:
559: suffix: mumps_3
560: requires: mumps
561: args: -mat_solver_type mumps
562: output_file: output/ex125_mumps_seq.out
564: test:
565: suffix: mumps_3_nest
566: requires: mumps
567: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
568: output_file: output/ex125_mumps_seq.out
570: test:
571: suffix: mumps_4
572: nsize: 3
573: requires: mumps
574: args: -mat_solver_type mumps
575: output_file: output/ex125_mumps_par.out
577: test:
578: suffix: mumps_4_nest
579: nsize: 3
580: requires: mumps
581: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
582: output_file: output/ex125_mumps_par.out
584: test:
585: suffix: mumps_5
586: nsize: 3
587: requires: mumps
588: args: -mat_solver_type mumps -cholesky
589: output_file: output/ex125_mumps_par_cholesky.out
591: test:
592: suffix: mumps_5_nest
593: nsize: 3
594: requires: mumps
595: args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}}
596: output_file: output/ex125_mumps_par_cholesky.out
598: test:
599: suffix: superlu
600: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
601: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
602: output_file: output/ex125_superlu.out
604: test:
605: suffix: superlu_dist
606: nsize: {{1 3}}
607: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
608: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
609: output_file: output/ex125_superlu_dist.out
611: test:
612: suffix: superlu_dist_2
613: nsize: {{1 3}}
614: requires: superlu_dist !complex
615: args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
616: output_file: output/ex125_superlu_dist.out
618: test:
619: suffix: superlu_dist_3
620: nsize: {{1 3}}
621: requires: superlu_dist !complex
622: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
623: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
624: output_file: output/ex125_superlu_dist_nonsymmetric.out
626: test:
627: suffix: superlu_dist_complex
628: nsize: 3
629: requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
630: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
631: output_file: output/ex125_superlu_dist_complex.out
633: test:
634: suffix: superlu_dist_complex_2
635: nsize: 3
636: requires: superlu_dist complex
637: args: -mat_solver_type superlu_dist
638: output_file: output/ex125_superlu_dist_complex_2.out
640: test:
641: suffix: cusparse
642: requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
643: #todo: fix the bug with cholesky
644: #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
645: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}
647: test:
648: suffix: cusparse_2
649: requires: cuda
650: args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}
652: TEST*/