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