Actual source code: ex2.c
2: static char help[] = "Tests repeated solving linear system on 2 by 2 matrix provided by MUMPS developer, Dec 17, 2012.\n\n";
3: /*
4: We have investigated the problem further, and we have
5: been able to reproduce it and obtain an erroneous
6: solution with an even smaller, 2x2, matrix:
7: [1 2]
8: [2 3]
9: and a right-hand side vector with all ones (1,1)
10: The correct solution is the vector (-1,1), in both solves.
12: mpiexec -n 2 ./ex2 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_7 6 -mat_mumps_cntl_1 0.99
14: With this combination of options, I get off-diagonal pivots during the
15: factorization, which is the cause of the problem (different isol_loc
16: returned in the second solve, whereas, as I understand it, Petsc expects
17: isol_loc not to change between successive solves).
18: */
20: #include <petscksp.h>
22: int main(int argc, char **args)
23: {
24: Mat C;
25: PetscInt N = 2, rowidx, colidx;
26: Vec u, b, r;
27: KSP ksp;
28: PetscReal norm;
29: PetscMPIInt rank, size;
30: PetscScalar v;
33: PetscInitialize(&argc, &args, (char *)0, help);
34: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
35: MPI_Comm_size(PETSC_COMM_WORLD, &size);
37: /* create stiffness matrix C = [1 2; 2 3] */
38: MatCreate(PETSC_COMM_WORLD, &C);
39: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
40: MatSetFromOptions(C);
41: MatSetUp(C);
42: if (rank == 0) {
43: rowidx = 0;
44: colidx = 0;
45: v = 1.0;
46: MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES);
47: rowidx = 0;
48: colidx = 1;
49: v = 2.0;
50: MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES);
52: rowidx = 1;
53: colidx = 0;
54: v = 2.0;
55: MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES);
56: rowidx = 1;
57: colidx = 1;
58: v = 3.0;
59: MatSetValues(C, 1, &rowidx, 1, &colidx, &v, INSERT_VALUES);
60: }
61: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
64: /* create right hand side and solution */
65: VecCreate(PETSC_COMM_WORLD, &u);
66: VecSetSizes(u, PETSC_DECIDE, N);
67: VecSetFromOptions(u);
68: VecDuplicate(u, &b);
69: VecDuplicate(u, &r);
70: VecSet(u, 0.0);
71: VecSet(b, 1.0);
73: /* solve linear system C*u = b */
74: KSPCreate(PETSC_COMM_WORLD, &ksp);
75: KSPSetOperators(ksp, C, C);
76: KSPSetFromOptions(ksp);
77: KSPSolve(ksp, b, u);
79: /* check residual r = C*u - b */
80: MatMult(C, u, r);
81: VecAXPY(r, -1.0, b);
82: VecNorm(r, NORM_2, &norm);
83: PetscPrintf(PETSC_COMM_WORLD, "|| C*u - b|| = %g\n", (double)norm);
85: /* solve C^T*u = b twice */
86: KSPSolveTranspose(ksp, b, u);
87: /* check residual r = C^T*u - b */
88: MatMultTranspose(C, u, r);
89: VecAXPY(r, -1.0, b);
90: VecNorm(r, NORM_2, &norm);
91: PetscPrintf(PETSC_COMM_WORLD, "|| C^T*u - b|| = %g\n", (double)norm);
93: KSPSolveTranspose(ksp, b, u);
94: MatMultTranspose(C, u, r);
95: VecAXPY(r, -1.0, b);
96: VecNorm(r, NORM_2, &norm);
97: PetscPrintf(PETSC_COMM_WORLD, "|| C^T*u - b|| = %g\n", (double)norm);
99: /* solve C*u = b again */
100: KSPSolve(ksp, b, u);
101: MatMult(C, u, r);
102: VecAXPY(r, -1.0, b);
103: VecNorm(r, NORM_2, &norm);
104: PetscPrintf(PETSC_COMM_WORLD, "|| C*u - b|| = %g\n", (double)norm);
106: KSPDestroy(&ksp);
107: VecDestroy(&u);
108: VecDestroy(&r);
109: VecDestroy(&b);
110: MatDestroy(&C);
111: PetscFinalize();
112: return 0;
113: }