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