xref: /petsc/src/mat/tests/ex125.c (revision 78bc9606edda02d8c6ac97f0eac63545386596c2)
1c4762a1bSJed Brown static char help[] = "Tests MatSolve() and MatMatSolve() (interface to superlu_dist, mumps and mkl_pardiso).\n\
2*78bc9606SHong Zhang Example: mpiexec -n <np> ./ex125 -f <matrix binary file> -nrhs 4 -mat_solver_type <>\n\n";
3*78bc9606SHong Zhang 
4*78bc9606SHong Zhang /*
5*78bc9606SHong Zhang -mat_solver_type:
6*78bc9606SHong Zhang   superlu
7*78bc9606SHong Zhang   superlu_dist
8*78bc9606SHong Zhang   mumps
9*78bc9606SHong Zhang   mkl_pardiso
10*78bc9606SHong Zhang   cusparse
11*78bc9606SHong Zhang   petsc
12*78bc9606SHong Zhang */
13c4762a1bSJed Brown 
14c4762a1bSJed Brown #include <petscmat.h>
15c4762a1bSJed Brown 
16d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
17d71ae5a4SJacob Faibussowitsch {
18b18964edSHong Zhang   Mat           A, RHS = NULL, RHS1 = NULL, C, F, X;
19c4762a1bSJed Brown   Vec           u, x, b;
20c4762a1bSJed Brown   PetscMPIInt   size;
21*78bc9606SHong Zhang   PetscInt      m, n, nfact, nsolve, nrhs, ipack = 5;
22c4762a1bSJed Brown   PetscReal     norm, tol = 1.e-10;
23c4762a1bSJed Brown   IS            perm, iperm;
24c4762a1bSJed Brown   MatFactorInfo info;
25c4762a1bSJed Brown   PetscRandom   rand;
26*78bc9606SHong Zhang   PetscBool     flg, symm, testMatSolve = PETSC_TRUE, testMatMatSolve = PETSC_TRUE, testMatMatSolveTranspose = PETSC_TRUE, testMatSolveTranspose = PETSC_TRUE, match = PETSC_FALSE;
27c4762a1bSJed Brown   PetscBool     chol = PETSC_FALSE, view = PETSC_FALSE, matsolvexx = PETSC_FALSE;
28c4762a1bSJed Brown #if defined(PETSC_HAVE_MUMPS)
29c4762a1bSJed Brown   PetscBool test_mumps_opts = PETSC_FALSE;
30c4762a1bSJed Brown #endif
31c4762a1bSJed Brown   PetscViewer fd;                       /* viewer */
32c4762a1bSJed Brown   char        file[PETSC_MAX_PATH_LEN]; /* input file name */
33*78bc9606SHong Zhang   char        pack[PETSC_MAX_PATH_LEN];
34c4762a1bSJed Brown 
35327415f7SBarry Smith   PetscFunctionBeginUser;
369566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
379566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
38c4762a1bSJed Brown 
39c4762a1bSJed Brown   /* Determine file from which we read the matrix A */
409566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
419a14fc28SStefano Zampini   if (flg) { /* Load matrix A */
429566063dSJacob Faibussowitsch     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
439566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
449566063dSJacob Faibussowitsch     PetscCall(MatSetFromOptions(A));
459566063dSJacob Faibussowitsch     PetscCall(MatLoad(A, fd));
469566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&fd));
479a14fc28SStefano Zampini   } else {
489a14fc28SStefano Zampini     n = 13;
499566063dSJacob Faibussowitsch     PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
509566063dSJacob Faibussowitsch     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
519566063dSJacob Faibussowitsch     PetscCall(MatSetType(A, MATAIJ));
529566063dSJacob Faibussowitsch     PetscCall(MatSetFromOptions(A));
539566063dSJacob Faibussowitsch     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
549566063dSJacob Faibussowitsch     PetscCall(MatSetUp(A));
559566063dSJacob Faibussowitsch     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
569566063dSJacob Faibussowitsch     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
579566063dSJacob Faibussowitsch     PetscCall(MatShift(A, 1.0));
589a14fc28SStefano Zampini   }
599566063dSJacob Faibussowitsch   PetscCall(MatGetLocalSize(A, &m, &n));
60e00437b9SBarry Smith   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
61c4762a1bSJed Brown 
62c4762a1bSJed Brown   /* if A is symmetric, set its flag -- required by MatGetInertia() */
63*78bc9606SHong Zhang   PetscCall(MatIsSymmetric(A, 0.0, &symm));
64*78bc9606SHong Zhang   PetscCall(MatSetOption(A, MAT_SYMMETRIC, symm));
65c4762a1bSJed Brown 
669566063dSJacob Faibussowitsch   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
67c4762a1bSJed Brown 
68a5b23f4aSJose E. Roman   /* Create dense matrix C and X; C holds true solution with identical columns */
69c4762a1bSJed Brown   nrhs = 2;
709566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
719566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ex125: nrhs %" PetscInt_FMT "\n", nrhs));
729566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
739566063dSJacob Faibussowitsch   PetscCall(MatSetOptionsPrefix(C, "rhs_"));
749566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
759566063dSJacob Faibussowitsch   PetscCall(MatSetType(C, MATDENSE));
769566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(C));
779566063dSJacob Faibussowitsch   PetscCall(MatSetUp(C));
78c4762a1bSJed Brown 
799566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_factor", &view, NULL));
809566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolve", &testMatMatSolve, NULL));
81*78bc9606SHong Zhang   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matmatsolvetranspose", &testMatMatSolveTranspose, NULL));
82*78bc9606SHong Zhang   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matsolvetranspose", &testMatSolveTranspose, NULL));
839566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-cholesky", &chol, NULL));
84c4762a1bSJed Brown #if defined(PETSC_HAVE_MUMPS)
859566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_mumps_opts", &test_mumps_opts, NULL));
86c4762a1bSJed Brown #endif
87c4762a1bSJed Brown 
889566063dSJacob Faibussowitsch   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
899566063dSJacob Faibussowitsch   PetscCall(PetscRandomSetFromOptions(rand));
909566063dSJacob Faibussowitsch   PetscCall(MatSetRandom(C, rand));
919566063dSJacob Faibussowitsch   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
92c4762a1bSJed Brown 
93c4762a1bSJed Brown   /* Create vectors */
949566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(A, &x, &b));
959566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(x, &u)); /* save the true solution */
96c4762a1bSJed Brown 
97c4762a1bSJed Brown   /* Test Factorization */
989566063dSJacob Faibussowitsch   PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
99c4762a1bSJed Brown 
100*78bc9606SHong Zhang   PetscCall(PetscOptionsGetString(NULL, NULL, "-mat_solver_type", pack, sizeof(pack), NULL));
101c4762a1bSJed Brown #if defined(PETSC_HAVE_SUPERLU)
102*78bc9606SHong Zhang   PetscCall(PetscStrcmp(MATSOLVERSUPERLU, pack, &match));
103*78bc9606SHong Zhang   if (match) {
10428b400f6SJacob Faibussowitsch     PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
1059566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU LU:\n"));
1069566063dSJacob Faibussowitsch     PetscCall(MatGetFactor(A, MATSOLVERSUPERLU, MAT_FACTOR_LU, &F));
107*78bc9606SHong Zhang     matsolvexx = PETSC_FALSE; /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix, need further work */
108*78bc9606SHong Zhang     ipack      = 0;
109*78bc9606SHong Zhang     goto skipoptions;
110*78bc9606SHong Zhang   }
111c4762a1bSJed Brown #endif
112c4762a1bSJed Brown #if defined(PETSC_HAVE_SUPERLU_DIST)
113*78bc9606SHong Zhang   PetscCall(PetscStrcmp(MATSOLVERSUPERLU_DIST, pack, &match));
114*78bc9606SHong Zhang   if (match) {
11528b400f6SJacob Faibussowitsch     PetscCheck(!chol, PETSC_COMM_WORLD, PETSC_ERR_SUP, "SuperLU does not provide Cholesky!");
1169566063dSJacob Faibussowitsch     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SUPERLU_DIST LU:\n"));
1179566063dSJacob Faibussowitsch     PetscCall(MatGetFactor(A, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &F));
118c4762a1bSJed Brown     matsolvexx = PETSC_TRUE;
119*78bc9606SHong Zhang     if (symm) { /* A is symmetric */
120*78bc9606SHong Zhang       testMatMatSolveTranspose = PETSC_TRUE;
121*78bc9606SHong Zhang       testMatSolveTranspose    = PETSC_TRUE;
122*78bc9606SHong Zhang     } else { /* superlu_dist does not support solving A^t x = rhs yet */
123*78bc9606SHong Zhang       testMatMatSolveTranspose = PETSC_FALSE;
124*78bc9606SHong Zhang       testMatSolveTranspose    = PETSC_FALSE;
125*78bc9606SHong Zhang     }
126*78bc9606SHong Zhang     ipack = 1;
127*78bc9606SHong Zhang     goto skipoptions;
128*78bc9606SHong Zhang   }
129c4762a1bSJed Brown #endif
130c4762a1bSJed Brown #if defined(PETSC_HAVE_MUMPS)
131*78bc9606SHong Zhang   PetscCall(PetscStrcmp(MATSOLVERMUMPS, pack, &match));
132*78bc9606SHong Zhang   if (match) {
133c4762a1bSJed Brown     if (chol) {
1349566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS CHOLESKY:\n"));
1359566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_CHOLESKY, &F));
136c4762a1bSJed Brown     } else {
1379566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MUMPS LU:\n"));
1389566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERMUMPS, MAT_FACTOR_LU, &F));
139c4762a1bSJed Brown     }
140c4762a1bSJed Brown     matsolvexx = PETSC_TRUE;
141c4762a1bSJed Brown     if (test_mumps_opts) {
142c4762a1bSJed Brown       /* test mumps options */
143c4762a1bSJed Brown       PetscInt  icntl;
144c4762a1bSJed Brown       PetscReal cntl;
145c4762a1bSJed Brown 
146c4762a1bSJed Brown       icntl = 2; /* sequential matrix ordering */
1479566063dSJacob Faibussowitsch       PetscCall(MatMumpsSetIcntl(F, 7, icntl));
148c4762a1bSJed Brown 
149c4762a1bSJed Brown       cntl = 1.e-6; /* threshold for row pivot detection */
1509566063dSJacob Faibussowitsch       PetscCall(MatMumpsSetIcntl(F, 24, 1));
1519566063dSJacob Faibussowitsch       PetscCall(MatMumpsSetCntl(F, 3, cntl));
152c4762a1bSJed Brown     }
153*78bc9606SHong Zhang     ipack = 2;
154*78bc9606SHong Zhang     goto skipoptions;
155*78bc9606SHong Zhang   }
156c4762a1bSJed Brown #endif
157c4762a1bSJed Brown #if defined(PETSC_HAVE_MKL_PARDISO)
158*78bc9606SHong Zhang   PetscCall(PetscStrcmp(MATSOLVERMKL_PARDISO, pack, &match));
159*78bc9606SHong Zhang   if (match) {
160c4762a1bSJed Brown     if (chol) {
1619566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO CHOLESKY:\n"));
1629566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_CHOLESKY, &F));
163c4762a1bSJed Brown     } else {
1649566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " MKL_PARDISO LU:\n"));
1659566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERMKL_PARDISO, MAT_FACTOR_LU, &F));
166c4762a1bSJed Brown     }
167*78bc9606SHong Zhang     ipack = 3;
168*78bc9606SHong Zhang     goto skipoptions;
169*78bc9606SHong Zhang   }
170c4762a1bSJed Brown #endif
17138a8e8c1SStefano Zampini #if defined(PETSC_HAVE_CUDA)
172*78bc9606SHong Zhang   PetscCall(PetscStrcmp(MATSOLVERCUSPARSE, pack, &match));
173*78bc9606SHong Zhang   if (match) {
17438a8e8c1SStefano Zampini     if (chol) {
1759566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE CHOLESKY:\n"));
1769566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_CHOLESKY, &F));
17738a8e8c1SStefano Zampini     } else {
1789566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " CUSPARSE LU:\n"));
1799566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERCUSPARSE, MAT_FACTOR_LU, &F));
18038a8e8c1SStefano Zampini     }
181*78bc9606SHong Zhang     testMatSolveTranspose    = PETSC_FALSE;
182*78bc9606SHong Zhang     testMatMatSolveTranspose = PETSC_FALSE;
183*78bc9606SHong Zhang     ipack                    = 4;
184*78bc9606SHong Zhang     goto skipoptions;
185*78bc9606SHong Zhang   }
18638a8e8c1SStefano Zampini #endif
187*78bc9606SHong Zhang   /* PETSc */
188*78bc9606SHong Zhang   match = PETSC_TRUE;
189*78bc9606SHong Zhang   if (match) {
190c4762a1bSJed Brown     if (chol) {
1919566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC CHOLESKY:\n"));
1929566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
193c4762a1bSJed Brown     } else {
1949566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PETSC LU:\n"));
1959566063dSJacob Faibussowitsch       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
196c4762a1bSJed Brown     }
197c4762a1bSJed Brown     matsolvexx = PETSC_TRUE;
198*78bc9606SHong Zhang     ipack      = 5;
199*78bc9606SHong Zhang     goto skipoptions;
200c4762a1bSJed Brown   }
201c4762a1bSJed Brown 
202*78bc9606SHong Zhang skipoptions:
2039566063dSJacob Faibussowitsch   PetscCall(MatFactorInfoInitialize(&info));
204c4762a1bSJed Brown   info.fill      = 5.0;
205c4762a1bSJed Brown   info.shifttype = (PetscReal)MAT_SHIFT_NONE;
206c4762a1bSJed Brown   if (chol) {
2079566063dSJacob Faibussowitsch     PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
208c4762a1bSJed Brown   } else {
2099566063dSJacob Faibussowitsch     PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
210c4762a1bSJed Brown   }
211c4762a1bSJed Brown 
212c4762a1bSJed Brown   for (nfact = 0; nfact < 2; nfact++) {
213c4762a1bSJed Brown     if (chol) {
2149566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the CHOLESKY numfactorization \n", nfact));
2159566063dSJacob Faibussowitsch       PetscCall(MatCholeskyFactorNumeric(F, A, &info));
216c4762a1bSJed Brown     } else {
2179566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %" PetscInt_FMT "-the LU numfactorization \n", nfact));
2189566063dSJacob Faibussowitsch       PetscCall(MatLUFactorNumeric(F, A, &info));
219c4762a1bSJed Brown     }
220c4762a1bSJed Brown     if (view) {
2219566063dSJacob Faibussowitsch       PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO));
2229566063dSJacob Faibussowitsch       PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
2239566063dSJacob Faibussowitsch       PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
224c4762a1bSJed Brown       view = PETSC_FALSE;
225c4762a1bSJed Brown     }
226c4762a1bSJed Brown 
227c4762a1bSJed Brown #if defined(PETSC_HAVE_SUPERLU_DIST)
228c4762a1bSJed Brown     if (ipack == 1) { /* Test MatSuperluDistGetDiagU()
229c4762a1bSJed Brown        -- input: matrix factor F; output: main diagonal of matrix U on all processes */
230c4762a1bSJed Brown       PetscInt     M;
231c4762a1bSJed Brown       PetscScalar *diag;
232c4762a1bSJed Brown   #if !defined(PETSC_USE_COMPLEX)
233c4762a1bSJed Brown       PetscInt nneg, nzero, npos;
234c4762a1bSJed Brown   #endif
235c4762a1bSJed Brown 
2369566063dSJacob Faibussowitsch       PetscCall(MatGetSize(F, &M, NULL));
2379566063dSJacob Faibussowitsch       PetscCall(PetscMalloc1(M, &diag));
2389566063dSJacob Faibussowitsch       PetscCall(MatSuperluDistGetDiagU(F, diag));
2399566063dSJacob Faibussowitsch       PetscCall(PetscFree(diag));
240c4762a1bSJed Brown 
241c4762a1bSJed Brown   #if !defined(PETSC_USE_COMPLEX)
242c4762a1bSJed Brown       /* Test MatGetInertia() */
243*78bc9606SHong Zhang       if (symm) { /* A is symmetric */
2449566063dSJacob Faibussowitsch         PetscCall(MatGetInertia(F, &nneg, &nzero, &npos));
2459566063dSJacob Faibussowitsch         PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, " MatInertia: nneg: %" PetscInt_FMT ", nzero: %" PetscInt_FMT ", npos: %" PetscInt_FMT "\n", nneg, nzero, npos));
246*78bc9606SHong Zhang       }
247c4762a1bSJed Brown   #endif
248c4762a1bSJed Brown     }
249c4762a1bSJed Brown #endif
250c4762a1bSJed Brown 
251d47f36abSHong Zhang #if defined(PETSC_HAVE_MUMPS)
252d47f36abSHong Zhang     /* mumps interface allows repeated call of MatCholeskyFactorSymbolic(), while the succession calls do nothing */
253d47f36abSHong Zhang     if (ipack == 2) {
254d47f36abSHong Zhang       if (chol) {
2559566063dSJacob Faibussowitsch         PetscCall(MatCholeskyFactorSymbolic(F, A, perm, &info));
2569566063dSJacob Faibussowitsch         PetscCall(MatCholeskyFactorNumeric(F, A, &info));
257d47f36abSHong Zhang       } else {
2589566063dSJacob Faibussowitsch         PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, &info));
2599566063dSJacob Faibussowitsch         PetscCall(MatLUFactorNumeric(F, A, &info));
260d47f36abSHong Zhang       }
261d47f36abSHong Zhang     }
262d47f36abSHong Zhang #endif
263d47f36abSHong Zhang 
264b18964edSHong Zhang     /* Test MatMatSolve(), A X = B, where B can be dense or sparse */
265c4762a1bSJed Brown     if (testMatMatSolve) {
266c4762a1bSJed Brown       if (!nfact) {
2679566063dSJacob Faibussowitsch         PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
268c4762a1bSJed Brown       } else {
2699566063dSJacob Faibussowitsch         PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS));
270c4762a1bSJed Brown       }
271c4762a1bSJed Brown       for (nsolve = 0; nsolve < 2; nsolve++) {
2729566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatMatSolve \n", nsolve));
2739566063dSJacob Faibussowitsch         PetscCall(MatMatSolve(F, RHS, X));
274c4762a1bSJed Brown 
275c4762a1bSJed Brown         /* Check the error */
2769566063dSJacob Faibussowitsch         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
2779566063dSJacob Faibussowitsch         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
278b18964edSHong Zhang         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolve: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
279c4762a1bSJed Brown       }
280b18964edSHong Zhang 
281c4762a1bSJed Brown       if (matsolvexx) {
282c4762a1bSJed Brown         /* Test MatMatSolve(F,RHS,RHS), RHS is a dense matrix */
2839566063dSJacob Faibussowitsch         PetscCall(MatCopy(RHS, X, SAME_NONZERO_PATTERN));
2849566063dSJacob Faibussowitsch         PetscCall(MatMatSolve(F, X, X));
285c4762a1bSJed Brown         /* Check the error */
2869566063dSJacob Faibussowitsch         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
2879566063dSJacob Faibussowitsch         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
288b18964edSHong Zhang         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatMatSolve(F,RHS,RHS): Norm of error %g\n", (double)norm));
289c4762a1bSJed Brown       }
290c4762a1bSJed Brown 
291c4762a1bSJed Brown       if (ipack == 2 && size == 1) {
292c4762a1bSJed Brown         Mat spRHS, spRHST, RHST;
293c4762a1bSJed Brown 
2949566063dSJacob Faibussowitsch         PetscCall(MatTranspose(RHS, MAT_INITIAL_MATRIX, &RHST));
2959566063dSJacob Faibussowitsch         PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
2969566063dSJacob Faibussowitsch         PetscCall(MatCreateTranspose(spRHST, &spRHS));
297c4762a1bSJed Brown         for (nsolve = 0; nsolve < 2; nsolve++) {
2989566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the sparse MatMatSolve \n", nsolve));
2999566063dSJacob Faibussowitsch           PetscCall(MatMatSolve(F, spRHS, X));
300c4762a1bSJed Brown 
301c4762a1bSJed Brown           /* Check the error */
3029566063dSJacob Faibussowitsch           PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
3039566063dSJacob Faibussowitsch           PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
304b18964edSHong Zhang           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));
305c4762a1bSJed Brown         }
306b18964edSHong Zhang         PetscCall(MatDestroy(&spRHST));
307b18964edSHong Zhang         PetscCall(MatDestroy(&spRHS));
308b18964edSHong Zhang         PetscCall(MatDestroy(&RHST));
309b18964edSHong Zhang       }
310b18964edSHong Zhang     }
311b18964edSHong Zhang 
312b18964edSHong Zhang     /* Test testMatMatSolveTranspose(), A^T X = B, where B can be dense or sparse */
313b18964edSHong Zhang     if (testMatMatSolveTranspose) {
314b18964edSHong Zhang       if (!nfact) {
315b18964edSHong Zhang         PetscCall(MatTransposeMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS1));
316b18964edSHong Zhang       } else {
317b18964edSHong Zhang         PetscCall(MatTransposeMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS1));
318b18964edSHong Zhang       }
319b18964edSHong Zhang 
320b18964edSHong Zhang       for (nsolve = 0; nsolve < 2; nsolve++) {
321*78bc9606SHong Zhang         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatMatSolveTranspose\n", nsolve));
322b18964edSHong Zhang         PetscCall(MatMatSolveTranspose(F, RHS1, X));
323b18964edSHong Zhang 
324b18964edSHong Zhang         /* Check the error */
325b18964edSHong Zhang         PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
326b18964edSHong Zhang         PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
327b18964edSHong Zhang         if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "-the MatMatSolveTranspose: Norm of error %g, nsolve %" PetscInt_FMT "\n", nsolve, (double)norm, nsolve));
328b18964edSHong Zhang       }
329b18964edSHong Zhang 
330b18964edSHong Zhang       if (ipack == 2 && size == 1) {
331b18964edSHong Zhang         Mat spRHS, spRHST, RHST;
332b18964edSHong Zhang 
333b18964edSHong Zhang         PetscCall(MatTranspose(RHS1, MAT_INITIAL_MATRIX, &RHST));
334b18964edSHong Zhang         PetscCall(MatConvert(RHST, MATAIJ, MAT_INITIAL_MATRIX, &spRHST));
335b18964edSHong Zhang         PetscCall(MatCreateTranspose(spRHST, &spRHS));
336b18964edSHong Zhang         for (nsolve = 0; nsolve < 2; nsolve++) {
337b18964edSHong Zhang           PetscCall(MatMatSolveTranspose(F, spRHS, X));
338b18964edSHong Zhang 
339b18964edSHong Zhang           /* Check the error */
340b18964edSHong Zhang           PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
341b18964edSHong Zhang           PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
342b18964edSHong Zhang           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));
343c4762a1bSJed Brown         }
3449566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&spRHST));
3459566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&spRHS));
3469566063dSJacob Faibussowitsch         PetscCall(MatDestroy(&RHST));
347c4762a1bSJed Brown       }
348c4762a1bSJed Brown     }
349c4762a1bSJed Brown 
350c4762a1bSJed Brown     /* Test MatSolve() */
351c4762a1bSJed Brown     if (testMatSolve) {
352c4762a1bSJed Brown       for (nsolve = 0; nsolve < 2; nsolve++) {
3539566063dSJacob Faibussowitsch         PetscCall(VecSetRandom(x, rand));
3549566063dSJacob Faibussowitsch         PetscCall(VecCopy(x, u));
3559566063dSJacob Faibussowitsch         PetscCall(MatMult(A, x, b));
356c4762a1bSJed Brown 
3579566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatSolve \n", nsolve));
3589566063dSJacob Faibussowitsch         PetscCall(MatSolve(F, b, x));
359c4762a1bSJed Brown 
360c4762a1bSJed Brown         /* Check the error */
3619566063dSJacob Faibussowitsch         PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
3629566063dSJacob Faibussowitsch         PetscCall(VecNorm(u, NORM_2, &norm));
363c4762a1bSJed Brown         if (norm > tol) {
364c4762a1bSJed Brown           PetscReal resi;
3659566063dSJacob Faibussowitsch           PetscCall(MatMult(A, x, u));    /* u = A*x */
3669566063dSJacob Faibussowitsch           PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
3679566063dSJacob Faibussowitsch           PetscCall(VecNorm(u, NORM_2, &resi));
3689566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolve: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
369c4762a1bSJed Brown         }
370c4762a1bSJed Brown       }
371c4762a1bSJed Brown     }
372*78bc9606SHong Zhang 
373*78bc9606SHong Zhang     /* Test MatSolveTranspose() */
374*78bc9606SHong Zhang     if (testMatSolveTranspose) {
375*78bc9606SHong Zhang       for (nsolve = 0; nsolve < 2; nsolve++) {
376*78bc9606SHong Zhang         PetscCall(VecSetRandom(x, rand));
377*78bc9606SHong Zhang         PetscCall(VecCopy(x, u));
378*78bc9606SHong Zhang         PetscCall(MatMultTranspose(A, x, b));
379*78bc9606SHong Zhang 
380*78bc9606SHong Zhang         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "   %" PetscInt_FMT "-the MatSolveTranspose\n", nsolve));
381*78bc9606SHong Zhang         PetscCall(MatSolveTranspose(F, b, x));
382*78bc9606SHong Zhang 
383*78bc9606SHong Zhang         /* Check the error */
384*78bc9606SHong Zhang         PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
385*78bc9606SHong Zhang         PetscCall(VecNorm(u, NORM_2, &norm));
386*78bc9606SHong Zhang         if (norm > tol) {
387*78bc9606SHong Zhang           PetscReal resi;
388*78bc9606SHong Zhang           PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
389*78bc9606SHong Zhang           PetscCall(VecNorm(u, NORM_2, &resi));
390*78bc9606SHong Zhang           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSolveTranspose: Norm of error %g, resi %g, numfact %" PetscInt_FMT "\n", (double)norm, (double)resi, nfact));
391*78bc9606SHong Zhang         }
392*78bc9606SHong Zhang       }
393*78bc9606SHong Zhang     }
394c4762a1bSJed Brown   }
395c4762a1bSJed Brown 
396c4762a1bSJed Brown   /* Free data structures */
3979566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
3989566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&C));
3999566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&F));
4009566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&X));
4019566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&RHS));
402b18964edSHong Zhang   PetscCall(MatDestroy(&RHS1));
403c4762a1bSJed Brown 
4049566063dSJacob Faibussowitsch   PetscCall(PetscRandomDestroy(&rand));
4059566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&perm));
4069566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&iperm));
4079566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&x));
4089566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&b));
4099566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&u));
4109566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
411b122ec5aSJacob Faibussowitsch   return 0;
412c4762a1bSJed Brown }
413c4762a1bSJed Brown 
414c4762a1bSJed Brown /*TEST
415c4762a1bSJed Brown 
416c4762a1bSJed Brown    test:
417dfd57a17SPierre Jolivet       requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
418*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type petsc
419c4762a1bSJed Brown       output_file: output/ex125.out
420c4762a1bSJed Brown 
421c4762a1bSJed Brown    test:
4229a14fc28SStefano Zampini       suffix: 2
423*78bc9606SHong Zhang       args: -mat_solver_type petsc
4249a14fc28SStefano Zampini       output_file: output/ex125.out
4259a14fc28SStefano Zampini 
4269a14fc28SStefano Zampini    test:
427c4762a1bSJed Brown       suffix: mkl_pardiso
428dfd57a17SPierre Jolivet       requires: mkl_pardiso datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
429*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mkl_pardiso
430c4762a1bSJed Brown 
431c4762a1bSJed Brown    test:
4329a14fc28SStefano Zampini       suffix: mkl_pardiso_2
4339a14fc28SStefano Zampini       requires: mkl_pardiso
434*78bc9606SHong Zhang       args: -mat_solver_type mkl_pardiso
4359a14fc28SStefano Zampini       output_file: output/ex125_mkl_pardiso.out
4369a14fc28SStefano Zampini 
4379a14fc28SStefano Zampini    test:
438c4762a1bSJed Brown       suffix: mumps
439dfd57a17SPierre Jolivet       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
440*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
441c4762a1bSJed Brown       output_file: output/ex125_mumps_seq.out
442c4762a1bSJed Brown 
443c4762a1bSJed Brown    test:
444c4762a1bSJed Brown       suffix: mumps_2
445c4762a1bSJed Brown       nsize: 3
446dfd57a17SPierre Jolivet       requires: mumps datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
447*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type mumps
448c4762a1bSJed Brown       output_file: output/ex125_mumps_par.out
449c4762a1bSJed Brown 
450c4762a1bSJed Brown    test:
4519a14fc28SStefano Zampini       suffix: mumps_3
4529a14fc28SStefano Zampini       requires: mumps
453*78bc9606SHong Zhang       args: -mat_solver_type mumps
4549a14fc28SStefano Zampini       output_file: output/ex125_mumps_seq.out
4559a14fc28SStefano Zampini 
4569a14fc28SStefano Zampini    test:
4579a14fc28SStefano Zampini       suffix: mumps_4
4589a14fc28SStefano Zampini       nsize: 3
4599a14fc28SStefano Zampini       requires: mumps
460*78bc9606SHong Zhang       args: -mat_solver_type mumps
4619a14fc28SStefano Zampini       output_file: output/ex125_mumps_par.out
4629a14fc28SStefano Zampini 
4639a14fc28SStefano Zampini    test:
464d47f36abSHong Zhang       suffix: mumps_5
465d47f36abSHong Zhang       nsize: 3
466d47f36abSHong Zhang       requires: mumps
467*78bc9606SHong Zhang       args: -mat_solver_type mumps -cholesky
468d47f36abSHong Zhang       output_file: output/ex125_mumps_par_cholesky.out
469d47f36abSHong Zhang 
470d47f36abSHong Zhang    test:
471*78bc9606SHong Zhang       suffix: superlu
472*78bc9606SHong Zhang       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu
473*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu
474*78bc9606SHong Zhang       output_file: output/ex125_superlu.out
475*78bc9606SHong Zhang 
476*78bc9606SHong Zhang    test:
477c4762a1bSJed Brown       suffix: superlu_dist
4789a14fc28SStefano Zampini       nsize: {{1 3}}
479d4783600SBarry Smith       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
480*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/small -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
481*78bc9606SHong Zhang       output_file: output/ex125_superlu_dist.out
482c4762a1bSJed Brown 
483c4762a1bSJed Brown    test:
484c4762a1bSJed Brown       suffix: superlu_dist_2
4859a14fc28SStefano Zampini       nsize: {{1 3}}
4869a14fc28SStefano Zampini       requires: superlu_dist !complex
487*78bc9606SHong Zhang       args: -n 36 -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
488c4762a1bSJed Brown       output_file: output/ex125_superlu_dist.out
489c4762a1bSJed Brown 
490c4762a1bSJed Brown    test:
491*78bc9606SHong Zhang       suffix: superlu_dist_3
492*78bc9606SHong Zhang       nsize: {{1 3}}
493*78bc9606SHong Zhang       requires: superlu_dist !complex
494*78bc9606SHong Zhang       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) superlu_dist
495*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/medium -mat_solver_type superlu_dist -mat_superlu_dist_rowperm NOROWPERM
496*78bc9606SHong Zhang       output_file: output/ex125_superlu_dist_nonsymmetric.out
497*78bc9606SHong Zhang 
498*78bc9606SHong Zhang    test:
499c4762a1bSJed Brown       suffix: superlu_dist_complex
500c4762a1bSJed Brown       nsize: 3
501d4783600SBarry Smith       requires: datafilespath double superlu_dist complex !defined(PETSC_USE_64BIT_INDICES)
502*78bc9606SHong Zhang       args: -f ${DATAFILESPATH}/matrices/farzad_B_rhs -mat_solver_type superlu_dist
503c4762a1bSJed Brown       output_file: output/ex125_superlu_dist_complex.out
504c4762a1bSJed Brown 
50538a8e8c1SStefano Zampini    test:
5069a14fc28SStefano Zampini       suffix: superlu_dist_complex_2
5079a14fc28SStefano Zampini       nsize: 3
5089a14fc28SStefano Zampini       requires: superlu_dist complex
509*78bc9606SHong Zhang       args: -mat_solver_type superlu_dist
510*78bc9606SHong Zhang       output_file: output/ex125_superlu_dist_complex_2.out
5119a14fc28SStefano Zampini 
5129a14fc28SStefano Zampini    test:
51338a8e8c1SStefano Zampini       suffix: cusparse
514dfd57a17SPierre Jolivet       requires: cuda datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
51573ff1584SJunchao Zhang       #todo: fix the bug with cholesky
516*78bc9606SHong Zhang       #args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0 1}separate output}
517*78bc9606SHong Zhang       args: -mat_type aijcusparse -f ${DATAFILESPATH}/matrices/small -mat_solver_type cusparse -cholesky {{0}separate output}
51838a8e8c1SStefano Zampini 
5199a14fc28SStefano Zampini    test:
5209a14fc28SStefano Zampini       suffix: cusparse_2
5219a14fc28SStefano Zampini       requires: cuda
522*78bc9606SHong Zhang       args: -mat_type aijcusparse -mat_solver_type cusparse -cholesky {{0 1}separate output}
5239a14fc28SStefano Zampini 
524c4762a1bSJed Brown TEST*/
525