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