1 static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() for computing inv(A) with MUMPS.\n\
2 Example: mpiexec -n <np> ./ex214 -displ \n\n";
3
4 #include <petscmat.h>
5
main(int argc,char ** args)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
18
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));
23
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
29
30 PetscCall(PetscOptionsGetBool(NULL, NULL, "-displ", &displ, NULL));
31
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));
37
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));
43
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));
70
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);
74
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));
83
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));
88
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));
91
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");
104
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 }
110
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));
119
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));
126
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));
137
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 }
143
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));
150
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));
172
173 PetscCall(MatMatTransposeSolve(F, spRHST, X));
174
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 }
179
180 /* Check the residual: R = A*X - RHS */
181 PetscCall(MatMatMult(A, X, MAT_REUSE_MATRIX, 2.0, &AX));
182
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));
186
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;
193
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));
198
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 }
211
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 }
221
222 /*TEST
223
224 test:
225 requires: mumps double !complex
226
227 test:
228 suffix: 2
229 requires: mumps double !complex
230 nsize: 2
231
232 TEST*/
233