Actual source code: ex214.c

  1: static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() for computing inv(A) with MUMPS.\n\
  2: Example: mpiexec -n <np> ./ex214 -displ \n\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   PetscMPIInt size, rank;
  9: #if defined(PETSC_HAVE_MUMPS)
 10:   Mat         A, RHS, C, F, X, AX, spRHST;
 11:   PetscInt    m, n, nrhs, M, N, i, Istart, Iend, Ii, j, J, test;
 12:   PetscScalar v;
 13:   PetscReal   norm, tol = PETSC_SQRT_MACHINE_EPSILON;
 14:   PetscRandom rand;
 15:   PetscBool   displ = PETSC_FALSE;
 16:   char        solver[256];
 17: #endif

 19:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 22:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 24: #if !defined(PETSC_HAVE_MUMPS)
 25:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "This example requires MUMPS, exit...\n"));
 26:   PetscCall(PetscFinalize());
 27:   return 0;
 28: #else

 30:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));

 32:   /* Create matrix A */
 33:   m = 4;
 34:   n = 4;
 35:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 36:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 38:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 39:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 40:   PetscCall(MatSetFromOptions(A));
 41:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
 42:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));

 44:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
 45:   for (Ii = Istart; Ii < Iend; Ii++) {
 46:     v = -1.0;
 47:     i = Ii / n;
 48:     j = Ii - i * n;
 49:     if (i > 0) {
 50:       J = Ii - n;
 51:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 52:     }
 53:     if (i < m - 1) {
 54:       J = Ii + n;
 55:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 56:     }
 57:     if (j > 0) {
 58:       J = Ii - 1;
 59:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 60:     }
 61:     if (j < n - 1) {
 62:       J = Ii + 1;
 63:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES));
 64:     }
 65:     v = 4.0;
 66:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES));
 67:   }
 68:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 69:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 71:   PetscCall(MatGetLocalSize(A, &m, &n));
 72:   PetscCall(MatGetSize(A, &M, &N));
 73:   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:   /* Create dense matrix C and X; C holds true solution with identical columns */
 76:   nrhs = N;
 77:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
 78:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 79:   PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
 80:   PetscCall(MatSetType(C, MATDENSE));
 81:   PetscCall(MatSetFromOptions(C));
 82:   PetscCall(MatSetUp(C));

 84:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 85:   PetscCall(PetscRandomSetFromOptions(rand));
 86:   PetscCall(MatSetRandom(C, rand));
 87:   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));

 89:   PetscCall(PetscStrncpy(solver, MATSOLVERMUMPS, sizeof(solver)));
 90:   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:   for (test = 0; test < 2; test++) {
 93:     if (test == 0) {
 94:       /* Test LU Factorization */
 95:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "test LU factorization\n"));
 96:       PetscCall(MatGetFactor(A, solver, MAT_FACTOR_LU, &F));
 97:       PetscCall(MatLUFactorSymbolic(F, A, NULL, NULL, NULL));
 98:       PetscCall(MatLUFactorNumeric(F, A, NULL));
 99:     } else {
100:       /* Test Cholesky Factorization */
101:       PetscBool flg;
102:       PetscCall(MatIsSymmetric(A, 0.0, &flg));
103:       PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "A must be symmetric for Cholesky factorization");

105:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "test Cholesky factorization\n"));
106:       PetscCall(MatGetFactor(A, solver, MAT_FACTOR_CHOLESKY, &F));
107:       PetscCall(MatCholeskyFactorSymbolic(F, A, NULL, NULL));
108:       PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
109:     }

111:     /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
112:     /* ---------------------------------------------------------- */
113:     PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, 2.0, &RHS));
114:     PetscCall(MatMatSolve(F, RHS, X));
115:     /* Check the error */
116:     PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
117:     PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
118:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(1) MatMatSolve: Norm of error %g\n", norm));

120:     /* Test X=RHS */
121:     PetscCall(MatMatSolve(F, RHS, RHS));
122:     /* Check the error */
123:     PetscCall(MatAXPY(RHS, -1.0, C, SAME_NONZERO_PATTERN));
124:     PetscCall(MatNorm(RHS, NORM_FROBENIUS, &norm));
125:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(1.1) MatMatSolve: Norm of error %g\n", norm));

127:     /* (2) Test MatMatSolve() for inv(A) with dense RHS:
128:      RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
129:     /* -------------------------------------------------------------------- */
130:     PetscCall(MatZeroEntries(RHS));
131:     for (i = 0; i < nrhs; i++) {
132:       v = 1.0;
133:       PetscCall(MatSetValues(RHS, 1, &i, 1, &i, &v, INSERT_VALUES));
134:     }
135:     PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY));
136:     PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY));

138:     PetscCall(MatMatSolve(F, RHS, X));
139:     if (displ) {
140:       if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " \n(2) first %" PetscInt_FMT " columns of inv(A) with dense RHS:\n", nrhs));
141:       PetscCall(MatView(X, PETSC_VIEWER_STDOUT_WORLD));
142:     }

144:     /* Check the residual */
145:     PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, 2.0, &AX));
146:     PetscCall(MatAXPY(AX, -1.0, RHS, SAME_NONZERO_PATTERN));
147:     PetscCall(MatNorm(AX, NORM_INFINITY, &norm));
148:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(2) MatMatSolve: Norm of residual %g\n", norm));
149:     PetscCall(MatZeroEntries(X));

151:     /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
152:      spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
153:     /* --------------------------------------------------------------------------- */
154:     /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
155:      thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
156:     PetscCall(MatCreate(PETSC_COMM_WORLD, &spRHST));
157:     if (rank == 0) {
158:       /* MUMPS requires RHS be centralized on the host! */
159:       PetscCall(MatSetSizes(spRHST, nrhs, M, PETSC_DECIDE, PETSC_DECIDE));
160:     } else {
161:       PetscCall(MatSetSizes(spRHST, 0, 0, PETSC_DECIDE, PETSC_DECIDE));
162:     }
163:     PetscCall(MatSetType(spRHST, MATAIJ));
164:     PetscCall(MatSetFromOptions(spRHST));
165:     PetscCall(MatSetUp(spRHST));
166:     if (rank == 0) {
167:       v = 1.0;
168:       for (i = 0; i < nrhs; i++) PetscCall(MatSetValues(spRHST, 1, &i, 1, &i, &v, INSERT_VALUES));
169:     }
170:     PetscCall(MatAssemblyBegin(spRHST, MAT_FINAL_ASSEMBLY));
171:     PetscCall(MatAssemblyEnd(spRHST, MAT_FINAL_ASSEMBLY));

173:     PetscCall(MatMatTransposeSolve(F, spRHST, X));

175:     if (displ) {
176:       if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " \n(3) first %" PetscInt_FMT " columns of inv(A) with sparse RHS:\n", nrhs));
177:       PetscCall(MatView(X, PETSC_VIEWER_STDOUT_WORLD));
178:     }

180:     /* Check the residual: R = A*X - RHS */
181:     PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, 2.0, &AX));

183:     PetscCall(MatAXPY(AX, -1.0, RHS, SAME_NONZERO_PATTERN));
184:     PetscCall(MatNorm(AX, NORM_INFINITY, &norm));
185:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "(3) MatMatSolve: Norm of residual %g\n", norm));

187:     /* (4) Test MatMatSolve() for inv(A) with selected entries:
188:      input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
189:     /* --------------------------------------------------------------------------------- */
190:     if (nrhs == N) { /* mumps requires nrhs = n */
191:       /* Create spRHS on proc[0] */
192:       Mat spRHS = NULL;

194:       /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
195:       PetscCall(MatCreateTranspose(spRHST, &spRHS));
196:       PetscCall(MatMumpsGetInverse(F, spRHS));
197:       PetscCall(MatDestroy(&spRHS));

199:       PetscCall(MatMumpsGetInverseTranspose(F, spRHST));
200:       if (displ) {
201:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSelected entries of inv(A^T):\n"));
202:         PetscCall(MatView(spRHST, PETSC_VIEWER_STDOUT_WORLD));
203:       }
204:       PetscCall(MatDestroy(&spRHS));
205:     }
206:     PetscCall(MatDestroy(&AX));
207:     PetscCall(MatDestroy(&F));
208:     PetscCall(MatDestroy(&RHS));
209:     PetscCall(MatDestroy(&spRHST));
210:   }

212:   /* Free data structures */
213:   PetscCall(MatDestroy(&A));
214:   PetscCall(MatDestroy(&C));
215:   PetscCall(MatDestroy(&X));
216:   PetscCall(PetscRandomDestroy(&rand));
217:   PetscCall(PetscFinalize());
218:   return 0;
219: #endif
220: }

222: /*TEST

224:    test:
225:      requires: mumps double !complex

227:    test:
228:      suffix: 2
229:      requires: mumps double !complex
230:      nsize: 2

232: TEST*/