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*/