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) {
399: PetscCall(MatMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
400: } else {
401: PetscCall(MatMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS));
402: }
403: for (nsolve = 0; nsolve < 2; nsolve++) {
404: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
405: PetscCall(MatMatSolve(F, RHS, X));
407: /* Check the error */
408: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
409: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
410: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
411: }
413: if (matsolvexx) {
414: /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
415: PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
416: PetscCall(MatMatSolve(F, X, 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, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
421: }
423: if (ipack == 2 && size == 1) {
424: Mat spRHS, spRHST, RHST;
426: PetscCall(MatTranspose(RHS, 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(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
431: PetscCall(MatMatSolve(F, spRHS, X));
433: /* Check the error */
434: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
435: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
436: 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));
437: }
438: PetscCall(MatDestroy(&spRHST));
439: PetscCall(MatDestroy(&spRHS));
440: PetscCall(MatDestroy(&RHST));
441: }
442: }
444: /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
445: if (testMatMatSolveTranspose) {
446: if (!nfact) {
447: PetscCall(MatTransposeMatMult(Ae, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
448: } else {
449: PetscCall(MatTransposeMatMult(Ae, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
450: }
452: for (nsolve = 0; nsolve < 2; nsolve++) {
453: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
454: PetscCall(MatMatSolveTranspose(F, RHS1, X));
456: /* Check the error */
457: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
458: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
459: if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
460: }
462: if (ipack == 2 && size == 1) {
463: Mat spRHS, spRHST, RHST;
465: PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
466: PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
467: PetscCall(MatCreateTranspose(spRHST, &spRHS));
468: for (nsolve = 0; nsolve < 2; nsolve++) {
469: PetscCall(MatMatSolveTranspose(F, spRHS, X));
471: /* Check the error */
472: PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
473: PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
474: 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));
475: }
476: PetscCall(MatDestroy(&spRHST));
477: PetscCall(MatDestroy(&spRHS));
478: PetscCall(MatDestroy(&RHST));
479: }
480: }
482: /* Test MatSolve() */
483: if (testMatSolve) {
484: for (nsolve = 0; nsolve < 2; nsolve++) {
485: PetscCall(VecSetRandom(x, rand));
486: PetscCall(VecCopy(x, u));
487: PetscCall(MatMult(Ae, x, b));
489: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolve \n", nsolve));
490: PetscCall(MatSolve(F, b, x));
492: /* Check the error */
493: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
494: PetscCall(VecNorm(u, NORM_2, &norm));
495: if (norm > tol) {
496: PetscReal resi;
497: PetscCall(MatMult(Ae, x, u)); /* u = A*x */
498: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
499: PetscCall(VecNorm(u, NORM_2, &resi));
500: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
501: }
502: }
503: }
505: /* Test MatSolveTranspose() */
506: if (testMatSolveTranspose) {
507: for (nsolve = 0; nsolve < 2; nsolve++) {
508: PetscCall(VecSetRandom(x, rand));
509: PetscCall(VecCopy(x, u));
510: PetscCall(MatMultTranspose(Ae, x, b));
512: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
513: PetscCall(MatSolveTranspose(F, b, x));
515: /* Check the error */
516: PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
517: PetscCall(VecNorm(u, NORM_2, &norm));
518: if (norm > tol) {
519: PetscReal resi;
520: PetscCall(MatMultTranspose(Ae, x, u)); /* u = A*x */
521: PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
522: PetscCall(VecNorm(u, NORM_2, &resi));
523: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
524: }
525: }
526: }
527: }
529: /* Free data structures */
530: PetscCall(MatDestroy(&Ae));
531: PetscCall(MatDestroy(&A));
532: PetscCall(MatDestroy(&C));
533: PetscCall(MatDestroy(&F));
534: PetscCall(MatDestroy(&X));
535: PetscCall(MatDestroy(&RHS));
536: PetscCall(MatDestroy(&RHS1));
538: PetscCall(PetscRandomDestroy(&rand));
539: PetscCall(ISDestroy(&perm));
540: PetscCall(ISDestroy(&iperm));
541: PetscCall(VecDestroy(&x));
542: PetscCall(VecDestroy(&b));
543: PetscCall(VecDestroy(&u));
544: PetscCall(PetscFinalize());
545: return 0;
546: }
548: /*TEST
550: test:
551: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
552: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
553: output_file: output/ex125.out
555: test:
556: suffix: 2
557: args: -mat_solver_type petsc
558: output_file: output/ex125.out
560: test:
561: suffix: mkl_pardiso
562: requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
563: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso
565: test:
566: suffix: mkl_pardiso_2
567: requires: mkl_pardiso
568: args: -mat_solver_type mkl_pardiso
569: output_file: output/ex125_mkl_pardiso.out
571: testset:
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: requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
578: suffix: mumps_single
579: args: -pc_precision single -tol 1e-5
580: test:
581: suffix: mumps_double
582: args: -pc_precision double
584: testset:
585: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
586: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
587: output_file: output/ex125_mumps_seq.out
589: test:
590: requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
591: suffix: mumps_nest_single
592: args: -pc_precision single -tol 1e-4
593: test:
594: suffix: mumps_nest_double
595: args: -pc_precision double
597: testset:
598: nsize: 3
599: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
600: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
601: output_file: output/ex125_mumps_par.out
603: test:
604: requires: defined(PETSC_HAVE_MUMPS_MIXED_PRECISION)
605: suffix: mumps_2_single
606: args: -pc_precision single -tol 1e-5
607: test:
608: suffix: mumps_2_double
609: args: -pc_precision double
611: test:
612: suffix: mumps_2_nest
613: nsize: 3
614: requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
615: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
616: output_file: output/ex125_mumps_par.out
618: test:
619: suffix: mumps_3
620: requires: mumps
621: args: -mat_solver_type mumps
622: output_file: output/ex125_mumps_seq.out
624: testset:
625: requires: mumps
626: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
627: output_file: output/ex125_mumps_seq.out
629: test:
630: requires: !__float128
631: suffix: mumps_3_nest
632: test:
633: suffix: mumps_3_nest_fp128
634: requires: __float128
635: args: -tol 1e-8
637: test:
638: suffix: mumps_4
639: nsize: 3
640: requires: mumps
641: args: -mat_solver_type mumps
642: output_file: output/ex125_mumps_par.out
644: testset:
645: nsize: 3
646: requires: mumps
647: args: -mat_solver_type mumps -test_nest -test_nest_bordered {{0 1}}
648: output_file: output/ex125_mumps_par.out
650: test:
651: requires: !__float128
652: suffix: mumps_4_nest
653: test:
654: suffix: mumps_4_nest_fp128
655: requires: __float128
656: args: -tol 1e-8
658: test:
659: suffix: mumps_5
660: nsize: 3
661: requires: mumps
662: args: -mat_solver_type mumps -cholesky
663: output_file: output/ex125_mumps_par_cholesky.out
665: testset:
666: nsize: 3
667: requires: mumps
668: args: -mat_solver_type mumps -cholesky -test_nest -test_nest_bordered {{0 1}}
669: output_file: output/ex125_mumps_par_cholesky.out
671: test:
672: requires: !__float128
673: suffix: mumps_5_nest
674: test:
675: suffix: mumps_5_nest_fp128
676: requires: __float128
677: args: -tol 1e-8
679: test:
680: nsize: 2
681: requires: mumps
682: args: -mat_solver_type mumps -test_nest -test_nest_bordered -m 13 -n 13
683: output_file: output/ex125_mumps_par.out
685: test:
686: requires: !__float128
687: suffix: mumps_6
688: test:
689: suffix: mumps_6_fp128
690: requires: __float128
691: args: -tol 1e-8
693: test:
694: suffix: superlu
695: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
696: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
697: output_file: output/ex125_superlu.out
699: test:
700: suffix: superlu_dist
701: nsize: {{1 3}}
702: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
703: args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
704: output_file: output/ex125_superlu_dist.out
706: test:
707: suffix: superlu_dist_2
708: nsize: {{1 3}}
709: requires: superlu_dist !complex
710: args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
711: output_file: output/ex125_superlu_dist.out
713: test:
714: suffix: superlu_dist_3
715: nsize: {{1 3}}
716: requires: superlu_dist !complex
717: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
718: args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
719: output_file: output/ex125_superlu_dist_nonsymmetric.out
721: test:
722: suffix: superlu_dist_complex
723: nsize: 3
724: requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
725: args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
726: output_file: output/ex125_superlu_dist_complex.out
728: test:
729: suffix: superlu_dist_complex_2
730: nsize: 3
731: requires: superlu_dist complex
732: args: -mat_solver_type superlu_dist
733: output_file: output/ex125_superlu_dist_complex_2.out
735: test:
736: suffix: cusparse
737: requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
738: #TODO: fix the bug with cholesky
739: #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
740: args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}
742: test:
743: suffix: cusparse_2
744: requires: cuda
745: args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}
747: testset:
748: nsize: {{1 2}separate output}
749: requires: double !defined(PETSC_USE_64BIT_INDICES) datafilespath !complex
750: args: -f ${DATAFILESPATH}/matrices/mixed_poisson
751: test:
752: requires: superlu_dist TODO # superlu_dist is broken
753: suffix: saddle_point_superlu_dist
754: args: -mat_solver_type superlu_dist -mat_superlu_dist_rowperm {{norowperm largediag_mc64}} -test_inertia 0
755: test:
756: requires: mumps
757: suffix: saddle_point_mumps_lu
758: args: -mat_solver_type mumps -mat_mumps_icntl_14 100
759: test:
760: requires: mumps
761: suffix: saddle_point_mumps_cholesky
762: args: -cholesky -mat_solver_type mumps
764: TEST*/