Actual source code: ex2.c
petsc-3.12.0 2019-09-29
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;
26: PetscInt N = 2,rowidx,colidx;
27: Vec u,b,r;
28: KSP ksp;
29: PetscReal norm;
30: PetscMPIInt rank,size;
31: PetscScalar v;
33: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
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) {
43: rowidx = 0; colidx = 0; v = 1.0;
44: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
45: rowidx = 0; colidx = 1; v = 2.0;
46: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
48: rowidx = 1; colidx = 0; v = 2.0;
49: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
50: rowidx = 1; colidx = 1; v = 3.0;
51: MatSetValues(C,1,&rowidx,1,&colidx,&v,INSERT_VALUES);
52: }
53: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
54: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
56: /* create right hand side and solution */
57: VecCreate(PETSC_COMM_WORLD,&u);
58: VecSetSizes(u,PETSC_DECIDE,N);
59: VecSetFromOptions(u);
60: VecDuplicate(u,&b);
61: VecDuplicate(u,&r);
62: VecSet(u,0.0);
63: VecSet(b,1.0);
65: /* solve linear system C*u = b */
66: KSPCreate(PETSC_COMM_WORLD,&ksp);
67: KSPSetOperators(ksp,C,C);
68: KSPSetFromOptions(ksp);
69: KSPSolve(ksp,b,u);
71: /* check residual r = C*u - b */
72: MatMult(C,u,r);
73: VecAXPY(r,-1.0,b);
74: VecNorm(r,NORM_2,&norm);
75: PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);
77: /* solve C^T*u = b twice */
78: KSPSolveTranspose(ksp,b,u);
79: /* check residual r = C^T*u - b */
80: MatMultTranspose(C,u,r);
81: VecAXPY(r,-1.0,b);
82: VecNorm(r,NORM_2,&norm);
83: PetscPrintf(PETSC_COMM_WORLD,"|| C^T*u - b|| = %g\n",(double)norm);
85: KSPSolveTranspose(ksp,b,u);
86: MatMultTranspose(C,u,r);
87: VecAXPY(r,-1.0,b);
88: VecNorm(r,NORM_2,&norm);
89: PetscPrintf(PETSC_COMM_WORLD,"|| C^T*u - b|| = %g\n",(double)norm);
91: /* solve C*u = b again */
92: KSPSolve(ksp,b,u);
93: MatMult(C,u,r);
94: VecAXPY(r,-1.0,b);
95: VecNorm(r,NORM_2,&norm);
96: PetscPrintf(PETSC_COMM_WORLD,"|| C*u - b|| = %g\n",(double)norm);
98: KSPDestroy(&ksp);
99: VecDestroy(&u);
100: VecDestroy(&r);
101: VecDestroy(&b);
102: MatDestroy(&C);
103: PetscFinalize();
104: return ierr;
105: }