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
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));
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
31: PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));
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));
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));
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));
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);
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));
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));
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));
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");
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: }
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));
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));
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));
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: }
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));
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));
174: PetscCall(MatMatTransposeSolve(F, spRHST, X));
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: }
181: /* Check the residual: R = A*X - RHS */
182: PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, 2.0, &AX));
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));
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: PetscCall(MatCreateTranspose(spRHST, &spRHS));
197: PetscCall(MatMumpsGetInverse(F, spRHS));
198: PetscCall(MatDestroy(&spRHS));
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: }
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: }
223: /*TEST
225: test:
226: requires: mumps double !complex
228: test:
229: suffix: 2
230: requires: mumps double !complex
231: nsize: 2
233: TEST*/