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