1 2 static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() for computing inv(A) with MUMPS.\n\ 3 Example: mpiexec -n <np> ./ex214 -displ \n\n"; 4 5 #include <petscmat.h> 6 7 int main(int argc, char **args) 8 { 9 PetscMPIInt size, rank; 10 #if defined(PETSC_HAVE_MUMPS) 11 Mat A, RHS, C, F, X, AX, spRHST; 12 PetscInt m, n, nrhs, M, N, i, Istart, Iend, Ii, j, J, test; 13 PetscScalar v; 14 PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON; 15 PetscRandom rand; 16 PetscBool displ = PETSC_FALSE; 17 char solver[256]; 18 #endif 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 24 25 #if !defined(PETSC_HAVE_MUMPS) 26 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "This example requires MUMPS, exit...\n")); 27 PetscCall(PetscFinalize()); 28 return 0; 29 #else 30 31 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL)); 32 33 /* Create matrix A */ 34 m = 4; 35 n = 4; 36 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 37 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 38 39 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 40 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n)); 41 PetscCall(MatSetFromOptions(A)); 42 PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL)); 43 PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL)); 44 45 PetscCall(MatGetOwnershipRange(A, &Istart, &Iend)); 46 for (Ii = Istart; Ii < Iend; Ii++) { 47 v = -1.0; 48 i = Ii / n; 49 j = Ii - i * n; 50 if (i > 0) { 51 J = Ii - n; 52 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 53 } 54 if (i < m - 1) { 55 J = Ii + n; 56 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 57 } 58 if (j > 0) { 59 J = Ii - 1; 60 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 61 } 62 if (j < n - 1) { 63 J = Ii + 1; 64 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 65 } 66 v = 4.0; 67 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES)); 68 } 69 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 70 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 71 72 PetscCall(MatGetLocalSize(A, &m, &n)); 73 PetscCall(MatGetSize(A, &M, &N)); 74 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 75 76 /* Create dense matrix C and X; C holds true solution with identical columns */ 77 nrhs = N; 78 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 79 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 80 PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 81 PetscCall(MatSetType(C, MATDENSE)); 82 PetscCall(MatSetFromOptions(C)); 83 PetscCall(MatSetUp(C)); 84 85 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 86 PetscCall(PetscRandomSetFromOptions(rand)); 87 PetscCall(MatSetRandom(C, rand)); 88 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X)); 89 90 PetscCall(PetscStrncpy(solver, MATSOLVERMUMPS, sizeof(solver))); 91 if (rank == 0 && displ) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Solving with %s: nrhs %" PetscInt_FMT ", size mat %" PetscInt_FMT " x %" PetscInt_FMT "\n", solver, nrhs, M, N)); 92 93 for (test = 0; test < 2; test++) { 94 if (test == 0) { 95 /* Test LU Factorization */ 96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "test LU factorization\n")); 97 PetscCall(MatGetFactor(A, solver, MAT_FACTOR_LU, &F)); 98 PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL)); 99 PetscCall(MatLUFactorNumeric(F, A, NULL)); 100 } else { 101 /* Test Cholesky Factorization */ 102 PetscBool flg; 103 PetscCall(MatIsSymmetric(A, 0.0, &flg)); 104 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A must be symmetric for Cholesky factorization"); 105 106 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "test Cholesky factorization\n")); 107 PetscCall(MatGetFactor(A, solver, MAT_FACTOR_CHOLESKY, &F)); 108 PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL)); 109 PetscCall(MatCholeskyFactorNumeric(F, A, NULL)); 110 } 111 112 /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */ 113 /* ---------------------------------------------------------- */ 114 PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS)); 115 PetscCall(MatMatSolve(F, RHS, X)); 116 /* Check the error */ 117 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 118 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 119 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(1) MatMatSolve: Norm of error %g\n", norm)); 120 121 /* Test X=RHS */ 122 PetscCall(MatMatSolve(F, RHS, RHS)); 123 /* Check the error */ 124 PetscCall(MatAXPY(RHS, -1.0, C, SAME_NONZERO_PATTERN)); 125 PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm)); 126 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(1.1) MatMatSolve: Norm of error %g\n", norm)); 127 128 /* (2) Test MatMatSolve() for inv(A) with dense RHS: 129 RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */ 130 /* -------------------------------------------------------------------- */ 131 PetscCall(MatZeroEntries(RHS)); 132 for (i = 0; i < nrhs; i++) { 133 v = 1.0; 134 PetscCall(MatSetValues(RHS, 1, &i, 1, &i, &v, INSERT_VALUES)); 135 } 136 PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY)); 137 PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY)); 138 139 PetscCall(MatMatSolve(F, RHS, X)); 140 if (displ) { 141 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " \n(2) first %" PetscInt_FMT " columns of inv(A) with dense RHS:\n", nrhs)); 142 PetscCall(MatView(X, PETSC_VIEWER_STDOUT_WORLD)); 143 } 144 145 /* Check the residual */ 146 PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, 2.0, &AX)); 147 PetscCall(MatAXPY(AX, -1.0, RHS, SAME_NONZERO_PATTERN)); 148 PetscCall(MatNorm(AX, NORM_INFINITY, &norm)); 149 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(2) MatMatSolve: Norm of residual %g\n", norm)); 150 PetscCall(MatZeroEntries(X)); 151 152 /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host: 153 spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */ 154 /* --------------------------------------------------------------------------- */ 155 /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix, 156 thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */ 157 PetscCall(MatCreate(PETSC_COMM_WORLD, &spRHST)); 158 if (rank == 0) { 159 /* MUMPS requirs RHS be centralized on the host! */ 160 PetscCall(MatSetSizes(spRHST, nrhs, M, PETSC_DECIDE, PETSC_DECIDE)); 161 } else { 162 PetscCall(MatSetSizes(spRHST, 0, 0, PETSC_DECIDE, PETSC_DECIDE)); 163 } 164 PetscCall(MatSetType(spRHST, MATAIJ)); 165 PetscCall(MatSetFromOptions(spRHST)); 166 PetscCall(MatSetUp(spRHST)); 167 if (rank == 0) { 168 v = 1.0; 169 for (i = 0; i < nrhs; i++) PetscCall(MatSetValues(spRHST, 1, &i, 1, &i, &v, INSERT_VALUES)); 170 } 171 PetscCall(MatAssemblyBegin(spRHST, MAT_FINAL_ASSEMBLY)); 172 PetscCall(MatAssemblyEnd(spRHST, MAT_FINAL_ASSEMBLY)); 173 174 PetscCall(MatMatTransposeSolve(F, spRHST, X)); 175 176 if (displ) { 177 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " \n(3) first %" PetscInt_FMT " columns of inv(A) with sparse RHS:\n", nrhs)); 178 PetscCall(MatView(X, PETSC_VIEWER_STDOUT_WORLD)); 179 } 180 181 /* Check the residual: R = A*X - RHS */ 182 PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, 2.0, &AX)); 183 184 PetscCall(MatAXPY(AX, -1.0, RHS, SAME_NONZERO_PATTERN)); 185 PetscCall(MatNorm(AX, NORM_INFINITY, &norm)); 186 if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(3) MatMatSolve: Norm of residual %g\n", norm)); 187 188 /* (4) Test MatMatSolve() for inv(A) with selected entries: 189 input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */ 190 /* --------------------------------------------------------------------------------- */ 191 if (nrhs == N) { /* mumps requires nrhs = n */ 192 /* Create spRHS on proc[0] */ 193 Mat spRHS = NULL; 194 195 /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */ 196 PetscCall(MatCreateTranspose(spRHST, &spRHS)); 197 PetscCall(MatMumpsGetInverse(F, spRHS)); 198 PetscCall(MatDestroy(&spRHS)); 199 200 PetscCall(MatMumpsGetInverseTranspose(F, spRHST)); 201 if (displ) { 202 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSelected entries of inv(A^T):\n")); 203 PetscCall(MatView(spRHST, PETSC_VIEWER_STDOUT_WORLD)); 204 } 205 PetscCall(MatDestroy(&spRHS)); 206 } 207 PetscCall(MatDestroy(&AX)); 208 PetscCall(MatDestroy(&F)); 209 PetscCall(MatDestroy(&RHS)); 210 PetscCall(MatDestroy(&spRHST)); 211 } 212 213 /* Free data structures */ 214 PetscCall(MatDestroy(&A)); 215 PetscCall(MatDestroy(&C)); 216 PetscCall(MatDestroy(&X)); 217 PetscCall(PetscRandomDestroy(&rand)); 218 PetscCall(PetscFinalize()); 219 return 0; 220 #endif 221 } 222 223 /*TEST 224 225 test: 226 requires: mumps double !complex 227 228 test: 229 suffix: 2 230 requires: mumps double !complex 231 nsize: 2 232 233 TEST*/ 234