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;
32: PetscInitialize(&argc,&args,(char*)0,help);
33: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
34: MPI_Comm_size(PETSC_COMM_WORLD,&size);
36: /* create stiffness matrix C = [1 2; 2 3] */
37: MatCreate(PETSC_COMM_WORLD,&C);
38: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
39: MatSetFromOptions(C);
40: MatSetUp(C);
41: if (rank == 0) {
42: rowidx = 0; colidx = 0; v = 1.0;
43: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
44: rowidx = 0; colidx = 1; v = 2.0;
45: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
47: rowidx = 1; colidx = 0; v = 2.0;
48: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
49: rowidx = 1; colidx = 1; v = 3.0;
50: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
51: }
52: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
53: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
55: /* create right hand side and solution */
56: VecCreate(PETSC_COMM_WORLD,&u);
57: VecSetSizes(u,PETSC_DECIDE,N);
58: VecSetFromOptions(u);
59: VecDuplicate(u,&b);
60: VecDuplicate(u,&r);
61: VecSet(u,0.0);
62: VecSet(b,1.0);
64: /* solve linear system C*u = b */
65: KSPCreate(PETSC_COMM_WORLD,&ksp);
66: KSPSetOperators(ksp,C,C);
67: KSPSetFromOptions(ksp);
68: KSPSolve(ksp,b,u);
70: /* check residual r = C*u - b */
71: MatMult(C,u,r);
72: VecAXPY(r,-1.0,b);
73: VecNorm(r,NORM_2,&norm);
74: PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);
76: /* solve C^T*u = b twice */
77: KSPSolveTranspose(ksp,b,u);
78: /* check residual r = C^T*u - b */
79: MatMultTranspose(C,u,r);
80: VecAXPY(r,-1.0,b);
81: VecNorm(r,NORM_2,&norm);
82: PetscPrintf(PETSC_COMM_WORLD,"|| C^T*u - b|| = %g\n",(double)norm);
84: KSPSolveTranspose(ksp,b,u);
85: MatMultTranspose(C,u,r);
86: VecAXPY(r,-1.0,b);
87: VecNorm(r,NORM_2,&norm);
88: PetscPrintf(PETSC_COMM_WORLD,"|| C^T*u - b|| = %g\n",(double)norm);
90: /* solve C*u = b again */
91: KSPSolve(ksp,b,u);
92: MatMult(C,u,r);
93: VecAXPY(r,-1.0,b);
94: VecNorm(r,NORM_2,&norm);
95: PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);
97: KSPDestroy(&ksp);
98: VecDestroy(&u);
99: VecDestroy(&r);
100: VecDestroy(&b);
101: MatDestroy(&C);
102: PetscFinalize();
103: return 0;
104: }