Actual source code: ex214.c


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

  5: #include <petscmat.h>

  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

 21:   PetscInitialize(&argc, &args, (char *)0, help);
 22:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

 31:   PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL);

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

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

 45:   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:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 53:     }
 54:     if (i < m - 1) {
 55:       J = Ii + n;
 56:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 57:     }
 58:     if (j > 0) {
 59:       J = Ii - 1;
 60:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 61:     }
 62:     if (j < n - 1) {
 63:       J = Ii + 1;
 64:       MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES);
 65:     }
 66:     v = 4.0;
 67:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES);
 68:   }
 69:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 70:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 72:   MatGetLocalSize(A, &m, &n);
 73:   MatGetSize(A, &M, &N);

 76:   /* Create dense matrix C and X; C holds true solution with identical columns */
 77:   nrhs = N;
 78:   PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL);
 79:   MatCreate(PETSC_COMM_WORLD, &C);
 80:   MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs);
 81:   MatSetType(C, MATDENSE);
 82:   MatSetFromOptions(C);
 83:   MatSetUp(C);

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

 90:   PetscStrcpy(solver, MATSOLVERMUMPS);
 91:   if (rank == 0 && displ) PetscPrintf(PETSC_COMM_SELF, "Solving with %s: nrhs %" PetscInt_FMT ", size mat %" PetscInt_FMT " x %" PetscInt_FMT "\n", solver, nrhs, M, N);

 93:   for (test = 0; test < 2; test++) {
 94:     if (test == 0) {
 95:       /* Test LU Factorization */
 96:       PetscPrintf(PETSC_COMM_WORLD, "test LU factorization\n");
 97:       MatGetFactor(A, solver, MAT_FACTOR_LU, &F);
 98:       MatLUFactorSymbolic(F, A, NULL, NULL, NULL);
 99:       MatLUFactorNumeric(F, A, NULL);
100:     } else {
101:       /* Test Cholesky Factorization */
102:       PetscBool flg;
103:       MatIsSymmetric(A, 0.0, &flg);

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

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

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

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:     MatZeroEntries(RHS);
132:     for (i = 0; i < nrhs; i++) {
133:       v = 1.0;
134:       MatSetValues(RHS, 1, &i, 1, &i, &v, INSERT_VALUES);
135:     }
136:     MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY);
137:     MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY);

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

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

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:     MatCreate(PETSC_COMM_WORLD, &spRHST);
158:     if (rank == 0) {
159:       /* MUMPS requirs RHS be centralized on the host! */
160:       MatSetSizes(spRHST, nrhs, M, PETSC_DECIDE, PETSC_DECIDE);
161:     } else {
162:       MatSetSizes(spRHST, 0, 0, PETSC_DECIDE, PETSC_DECIDE);
163:     }
164:     MatSetType(spRHST, MATAIJ);
165:     MatSetFromOptions(spRHST);
166:     MatSetUp(spRHST);
167:     if (rank == 0) {
168:       v = 1.0;
169:       for (i = 0; i < nrhs; i++) MatSetValues(spRHST, 1, &i, 1, &i, &v, INSERT_VALUES);
170:     }
171:     MatAssemblyBegin(spRHST, MAT_FINAL_ASSEMBLY);
172:     MatAssemblyEnd(spRHST, MAT_FINAL_ASSEMBLY);

174:     MatMatTransposeSolve(F, spRHST, X);

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

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

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

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;

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

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

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

223: /*TEST

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

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

233: TEST*/