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, test_inertia;
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, NULL, 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: test_inertia = symm;
94: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_inertia", &test_inertia, NULL));
96: PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL));
98: /* test MATNEST support */
99: flg = PETSC_FALSE;
100: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest", &flg, NULL));
101: if (flg) {
102: Mat B;
104: flg = PETSC_FALSE;
105: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest_bordered", &flg, NULL));
106: if (!flg) {
107: Mat mats[9] = {NULL, NULL, A, NULL, A, NULL, A, NULL, NULL};
109: /* Create a nested matrix representing
110: | 0 0 A |
111: | 0 A 0 |
112: | A 0 0 |
113: */
114: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 3, NULL, 3, NULL, mats, &B));
115: flg = PETSC_TRUE;
116: } else {
117: Mat mats[4];
119: /* Create a nested matrix representing
120: | Id R |
121: | R^t A |
122: */
123: PetscCall(MatGetSize(A, NULL, &n));
124: m = n + 12;
125: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
126: PetscCall(CreateIdentity(m, &mats[0]));
127: PetscCall(CreateRandom(m, n, &mats[1]));
128: mats[3] = A;
130: /* use CreateTranspose/CreateHermitianTranspose or explicit matrix for debugging purposes */
131: flg = PETSC_FALSE;
132: PetscCall(PetscOptionsGetBool(NULL, NULL, "-expl", &flg, NULL));
133: #if PetscDefined(USE_COMPLEX)
134: if (chol) { /* Hermitian transpose not supported by MUMPS Cholesky factor */
135: if (!flg) PetscCall(MatCreateTranspose(mats[1], &mats[2]));
136: else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
137: flg = PETSC_TRUE;
138: } else {
139: if (!flg) {
140: Mat B;
142: PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B));
143: PetscCall(MatCreateHermitianTranspose(B, &mats[2]));
144: PetscCall(MatDestroy(&B));
145: if (n == m) {
146: PetscCall(MatScale(mats[2], PetscCMPLX(4.0, -2.0)));
147: PetscCall(MatShift(mats[2], PetscCMPLX(-2.0, 1.0))); // mats[2] = (4 - 2i) B* - (2 - i) I
148: PetscCall(MatCreateHermitianTranspose(mats[2], &B));
149: PetscCall(MatDestroy(mats + 2));
150: PetscCall(MatScale(B, 0.5));
151: PetscCall(MatShift(B, PetscCMPLX(1.0, 0.5))); // B = 0.5 mats[2]* - (1 - 0.5i) I = (2 + i) B - (1 + 0.5i) I + (1 + 0.5i) I = (2 + i) B
152: PetscCall(MatCreateHermitianTranspose(B, &mats[2])); // mats[2] = B* = (2 - i) B*
153: PetscCall(MatDestroy(&B));
154: PetscCall(MatScale(mats[1], PetscCMPLX(2.0, 1.0))); // mats[1] = (2 + i) B = mats[2]*
155: } else flg = PETSC_TRUE;
156: } else PetscCall(MatHermitianTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
157: }
158: #else
159: if (!flg) {
160: Mat B;
162: PetscCall(MatDuplicate(mats[1], MAT_COPY_VALUES, &B));
163: PetscCall(MatCreateTranspose(B, &mats[2]));
164: PetscCall(MatDestroy(&B));
165: if (n == m) {
166: PetscCall(MatScale(mats[2], 4.0));
167: PetscCall(MatShift(mats[2], -2.0)); // mats[2] = 4 B' - 2 I
168: PetscCall(MatCreateTranspose(mats[2], &B));
169: PetscCall(MatDestroy(mats + 2));
170: PetscCall(MatScale(B, 0.5));
171: PetscCall(MatShift(B, 1.0)); // B = 0.5 mats[2]' + I = 0.5 (4 B' - 2 I)' + I = 2 B
172: PetscCall(MatCreateTranspose(B, &mats[2])); // mats[2] = B' = 2 B'
173: PetscCall(MatDestroy(&B));
174: PetscCall(MatScale(mats[1], 2.0)); // mats[1] = 2 B = mats[2]'
175: } else flg = PETSC_TRUE;
176: } else PetscCall(MatTranspose(mats[1], MAT_INITIAL_MATRIX, &mats[2]));
177: #endif
178: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, mats, &B));
179: PetscCall(MatDestroy(&mats[0]));
180: PetscCall(MatDestroy(&mats[1]));
181: PetscCall(MatDestroy(&mats[2]));
182: }
183: PetscCall(MatDestroy(&A));
184: A = B;
185: PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
187: /* not all the combinations of MatMat operations are supported by MATNEST. */
188: PetscCall(MatComputeOperator(A, MATAIJ, &Ae));
189: } else {
190: PetscCall(PetscObjectReference((PetscObject)A));
191: Ae = A;
192: flg = PETSC_TRUE;
193: }
194: PetscCall(MatGetLocalSize(A, &m, &n));
195: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
197: PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
198: PetscCall(MatViewFromOptions(Ae, NULL, "-A_view_expl"));
200: /* Create dense matrix C and X; C holds true solution with identical columns */
201: nrhs = 2;
202: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
203: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs));
204: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
205: PetscCall(MatSetOptionsPrefix(C, "rhs_"));
206: PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
207: PetscCall(MatSetType(C, MATDENSE));
208: PetscCall(MatSetFromOptions(C));
209: PetscCall(MatSetUp(C));
211: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL));
212: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL));
213: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL));
214: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL));
215: #if defined(PETSC_HAVE_MUMPS)
216: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL));
217: #endif
219: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
220: PetscCall(PetscRandomSetFromOptions(rand));
221: PetscCall(MatSetRandom(C, rand));
222: PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
224: /* Create vectors */
225: PetscCall(MatCreateVecs(A, &x, &b));
226: PetscCall(VecDuplicate(x, &u)); /* save the true solution */
228: /* Test Factorization */
229: if (flg) PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm)); // TODO FIXME: MatConvert_Nest_AIJ() does not support chained MatCreate[Hermitian]Transpose()
231: PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
232: #if defined(PETSC_HAVE_SUPERLU)
233: PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match));
234: if (match) {
235: PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
236: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n"));
237: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
238: matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */
239: ipack = 0;
240: goto skipoptions;
241: }
242: #endif
243: #if defined(PETSC_HAVE_SUPERLU_DIST)
244: PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match));
245: if (match) {
246: PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
247: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n"));
248: PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
249: matsolvexx = PETSC_TRUE;
250: if (symm) { /* A is symmetric */
251: testMatMatSolveTranspose = PETSC_TRUE;
252: testMatSolveTranspose = PETSC_TRUE;
253: } else { /* superlu_dist does not support solving A^t x = rhs yet */
254: testMatMatSolveTranspose = PETSC_FALSE;
255: testMatSolveTranspose = PETSC_FALSE;
256: }
257: ipack = 1;
258: goto skipoptions;
259: }
260: #endif
261: #if defined(PETSC_HAVE_MUMPS)
262: PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match));
263: if (match) {
264: if (chol) {
265: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n"));
266: PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F));
267: } else {
268: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n"));
269: PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
270: }
271: matsolvexx = PETSC_TRUE;
272: if (test_mumps_opts) {
273: /* test mumps options */
274: PetscInt icntl;
275: PetscReal cntl;
277: icntl = 2; /* sequential matrix ordering */
278: PetscCall(MatMumpsSetIcntl(F, 7, icntl));
280: cntl = 1.e-6; /* threshold for row pivot detection */
281: PetscCall(MatMumpsSetIcntl(F, 24, 1));
282: PetscCall(MatMumpsSetCntl(F, 3, cntl));
283: }
284: ipack = 2;
285: goto skipoptions;
286: }
287: #endif
288: #if defined(PETSC_HAVE_MKL_PARDISO)
289: PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match));
290: if (match) {
291: if (chol) {
292: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n"));
293: PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F));
294: } else {
295: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n"));
296: PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F));
297: }
298: ipack = 3;
299: goto skipoptions;
300: }
301: #endif
302: #if defined(PETSC_HAVE_CUDA)
303: PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match));
304: if (match) {
305: if (chol) {
306: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n"));
307: PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F));
308: } else {
309: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n"));
310: PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F));
311: }
312: testMatSolveTranspose = PETSC_FALSE;
313: testMatMatSolveTranspose = PETSC_FALSE;
314: ipack = 4;
315: goto skipoptions;
316: }
317: #endif
318: /* PETSc */
319: match = PETSC_TRUE;
320: if (match) {
321: if (chol) {
322: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n"));
323: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
324: } else {
325: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n"));
326: PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
327: }
328: matsolvexx = PETSC_TRUE;
329: ipack = 5;
330: goto skipoptions;
331: }
333: skipoptions:
334: PetscCall(MatFactorInfoInitialize(&info));
335: info.fill = 5.0;
336: info.shifttype = (PetscReal)MAT_SHIFT_NONE;
337: if (chol) {
338: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
339: } else {
340: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
341: }
343: for (nfact = 0; nfact < 2; nfact++) {
344: if (chol) {
345: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact));
346: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
347: } else {
348: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact));
349: PetscCall(MatLUFactorNumeric(F, A, &info));
350: }
351: if (view) {
352: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
353: PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
354: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
355: view = PETSC_FALSE;
356: }
358: #if defined(PETSC_HAVE_SUPERLU_DIST)
359: if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
360: -- input: matrix factor F; output: main diagonal of matrix U on all processes */
361: PetscInt M;
362: PetscScalar *diag;
363: #if !defined(PETSC_USE_COMPLEX)
364: PetscInt nneg, nzero, npos;
365: #endif
367: PetscCall(MatGetSize(F, &M, NULL));
368: PetscCall(PetscMalloc1(M, &diag));
369: PetscCall(MatSuperluDistGetDiagU(F, diag));
370: PetscCall(PetscFree(diag));
372: #if !defined(PETSC_USE_COMPLEX)
373: /* Test MatGetInertia() */
374: if (test_inertia) { /* A is symmetric */
375: PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
376: PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
377: }
378: #endif
379: }
380: #endif
382: #if defined(PETSC_HAVE_MUMPS)
383: /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
384: if (ipack == 2) {
385: if (chol) {
386: PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
387: PetscCall(MatCholeskyFactorNumeric(F, A, &info));
388: } else {
389: PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
390: PetscCall(MatLUFactorNumeric(F, A, &info));
391: }
392: }
393: #endif
395: /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
396: if (testMatMatSolve) {
397: if (!nfact) {
398: PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
399: } else {
400: PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS));
401: }
402: for (nsolve = 0; nsolve < 2; nsolve++) {
403: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
404: PetscCall(MatMatSolve(F, RHS, X));
406: /* Check the error */
407: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
408: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
409: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
410: }
412: if (matsolvexx) {
413: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
414: PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
415: PetscCall(MatMatSolve(F, X, X));
416: /* Check the error */
417: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
418: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
419: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
420: }
422: if (ipack == 2 && size == 1) {
423: Mat spRHS, spRHST, RHST;
425: PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST));
426: PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
427: PetscCall(MatCreateTranspose(spRHST, &spRHS));
428: for (nsolve = 0; nsolve < 2; nsolve++) {
429: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
430: PetscCall(MatMatSolve(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 MatMatSolve: 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 testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
444: if (testMatMatSolveTranspose) {
445: if (!nfact) {
446: PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
447: } else {
448: PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
449: }
451: for (nsolve = 0; nsolve < 2; nsolve++) {
452: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
453: PetscCall(MatMatSolveTranspose(F, RHS1, X));
455: /* Check the error */
456: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
457: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
458: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
459: }
461: if (ipack == 2 && size == 1) {
462: Mat spRHS, spRHST, RHST;
464: PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
465: PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
466: PetscCall(MatCreateTranspose(spRHST, &spRHS));
467: for (nsolve = 0; nsolve < 2; nsolve++) {
468: PetscCall(MatMatSolveTranspose(F, spRHS, X));
470: /* Check the error */
471: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
472: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
473: 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));
474: }
475: PetscCall(MatDestroy(&spRHST));
476: PetscCall(MatDestroy(&spRHS));
477: PetscCall(MatDestroy(&RHST));
478: }
479: }
481: /* Test MatSolve() */
482: if (testMatSolve) {
483: for (nsolve = 0; nsolve < 2; nsolve++) {
484: PetscCall(VecSetRandom(x, rand));
485: PetscCall(VecCopy(x, u));
486: PetscCall(MatMult(Ae, x, b));
488: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve));
489: PetscCall(MatSolve(F, b, x));
491: /* Check the error */
492: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
493: PetscCall(VecNorm(u, NORM_2, &norm));
494: if (norm > tol) {
495: PetscReal resi;
496: PetscCall(MatMult(Ae, x, u)); /* u = A*x */
497: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
498: PetscCall(VecNorm(u, NORM_2, &resi));
499: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
500: }
501: }
502: }
504: /* Test MatSolveTranspose() */
505: if (testMatSolveTranspose) {
506: for (nsolve = 0; nsolve < 2; nsolve++) {
507: PetscCall(VecSetRandom(x, rand));
508: PetscCall(VecCopy(x, u));
509: PetscCall(MatMultTranspose(Ae, x, b));
511: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
512: PetscCall(MatSolveTranspose(F, b, x));
514: /* Check the error */
515: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
516: PetscCall(VecNorm(u, NORM_2, &norm));
517: if (norm > tol) {
518: PetscReal resi;
519: PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */
520: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
521: PetscCall(VecNorm(u, NORM_2, &resi));
522: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
523: }
524: }
525: }
526: }
528: /* Free data structures */
529: PetscCall(MatDestroy(&Ae));
530: PetscCall(MatDestroy(&A));
531: PetscCall(MatDestroy(&C));
532: PetscCall(MatDestroy(&F));
533: PetscCall(MatDestroy(&X));
534: PetscCall(MatDestroy(&RHS));
535: PetscCall(MatDestroy(&RHS1));
537: PetscCall(PetscRandomDestroy(&rand));
538: PetscCall(ISDestroy(&perm));
539: PetscCall(ISDestroy(&iperm));
540: PetscCall(VecDestroy(&x));
541: PetscCall(VecDestroy(&b));
542: PetscCall(VecDestroy(&u));
543: PetscCall(PetscFinalize());
544: return 0;
545: }
547: /*TEST
549: test:
550: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
551: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
552: output_file: output/ex125.out
554: test:
555: suffix: 2
556: args: -mat_solver_type petsc
557: output_file: output/ex125.out
559: test:
560: suffix: mkl_pardiso
561: requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
562: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso
564: test:
565: suffix: mkl_pardiso_2
566: requires: mkl_pardiso
567: args: -mat_solver_type mkl_pardiso
568: output_file: output/ex125_mkl_pardiso.out
570: test:
571: suffix: mumps
572: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
573: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
574: output_file: output/ex125_mumps_seq.out
576: test:
577: suffix: mumps_nest
578: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
579: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
580: output_file: output/ex125_mumps_seq.out
582: test:
583: suffix: mumps_2
584: nsize: 3
585: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
586: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
587: output_file: output/ex125_mumps_par.out
589: test:
590: suffix: mumps_2_nest
591: nsize: 3
592: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
593: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
594: output_file: output/ex125_mumps_par.out
596: test:
597: suffix: mumps_3
598: requires: mumps
599: args: -mat_solver_type mumps
600: output_file: output/ex125_mumps_seq.out
602: test:
603: suffix: mumps_3_nest
604: requires: mumps
605: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
606: output_file: output/ex125_mumps_seq.out
608: test:
609: suffix: mumps_4
610: nsize: 3
611: requires: mumps
612: args: -mat_solver_type mumps
613: output_file: output/ex125_mumps_par.out
615: test:
616: suffix: mumps_4_nest
617: nsize: 3
618: requires: mumps
619: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
620: output_file: output/ex125_mumps_par.out
622: test:
623: suffix: mumps_5
624: nsize: 3
625: requires: mumps
626: args: -mat_solver_type mumps -cholesky
627: output_file: output/ex125_mumps_par_cholesky.out
629: test:
630: suffix: mumps_5_nest
631: nsize: 3
632: requires: mumps
633: args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}}
634: output_file: output/ex125_mumps_par_cholesky.out
636: test:
637: suffix: mumps_6
638: nsize: 2
639: requires: mumps
640: args: -mat_solver_type mumps -test_nest -test_nest_bordered -m 13 -n 13
641: output_file: output/ex125_mumps_par.out
643: test:
644: suffix: superlu
645: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
646: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
647: output_file: output/ex125_superlu.out
649: test:
650: suffix: superlu_dist
651: nsize: {{1 3}}
652: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
653: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
654: output_file: output/ex125_superlu_dist.out
656: test:
657: suffix: superlu_dist_2
658: nsize: {{1 3}}
659: requires: superlu_dist !complex
660: args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
661: output_file: output/ex125_superlu_dist.out
663: test:
664: suffix: superlu_dist_3
665: nsize: {{1 3}}
666: requires: superlu_dist !complex
667: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
668: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
669: output_file: output/ex125_superlu_dist_nonsymmetric.out
671: test:
672: suffix: superlu_dist_complex
673: nsize: 3
674: requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
675: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
676: output_file: output/ex125_superlu_dist_complex.out
678: test:
679: suffix: superlu_dist_complex_2
680: nsize: 3
681: requires: superlu_dist complex
682: args: -mat_solver_type superlu_dist
683: output_file: output/ex125_superlu_dist_complex_2.out
685: test:
686: suffix: cusparse
687: requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
688: #TODO: fix the bug with cholesky
689: #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
690: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}
692: test:
693: suffix: cusparse_2
694: requires: cuda
695: args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}
697: testset:
698: nsize: {{1 2}separate output}
699: requires: double !defined(PETSC_USE_64BIT_INDICES) datafilespath !complex
700: args: -f ${DATAFILESPATH}/matrices/mixed_poisson
701: test:
702: requires: superlu_dist TODO # superlu_dist is broken
703: suffix: saddle_point_superlu_dist
704: args: -mat_solver_type superlu_dist -mat_superlu_dist_rowperm {{norowperm largediag_mc64}} -test_inertia 0
705: test:
706: requires: mumps
707: suffix: saddle_point_mumps_lu
708: args: -mat_solver_type mumps -mat_mumps_icntl_14 100
709: test:
710: requires: mumps
711: suffix: saddle_point_mumps_cholesky
712: args: -cholesky -mat_solver_type mumps
714: TEST*/