Actual source code: ex2.c
1: static char help[] = "Test file for the PCFactorSetShiftType()\n";
2: /*
3: * Test file for the PCFactorSetShiftType() routine or -pc_factor_shift_type POSITIVE_DEFINITE option.
4: * The test matrix is the example from Kershaw's paper [J.Comp.Phys 1978]
5: * of a positive definite matrix for which ILU(0) will give a negative pivot.
6: * This means that the CG method will break down; the Manteuffel shift
7: * [Math. Comp. 1980] repairs this.
8: *
9: * Run the executable twice:
10: * 1/ without options: the iterative method diverges because of an
11: * indefinite preconditioner
12: * 2/ with -pc_factor_shift_positive_definite option (or comment in the PCFactorSetShiftType() line below):
13: * the method will now successfully converge.
14: */
16: #include <petscksp.h>
18: int main(int argc,char **argv)
19: {
20: KSP ksp;
21: PC pc;
22: Mat A,M;
23: Vec X,B,D;
24: MPI_Comm comm;
25: PetscScalar v;
26: KSPConvergedReason reason;
27: PetscInt i,j,its;
29: PetscInitialize(&argc,&argv,0,help);
30: comm = MPI_COMM_SELF;
32: /*
33: * Construct the Kershaw matrix
34: * and a suitable rhs / initial guess
35: */
36: MatCreateSeqAIJ(comm,4,4,4,0,&A);
37: VecCreateSeq(comm,4,&B);
38: VecDuplicate(B,&X);
39: for (i=0; i<4; i++) {
40: v = 3;
41: MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
42: v = 1;
43: VecSetValues(B,1,&i,&v,INSERT_VALUES);
44: VecSetValues(X,1,&i,&v,INSERT_VALUES);
45: }
47: i =0; v=0;
48: VecSetValues(X,1,&i,&v,INSERT_VALUES);
50: for (i=0; i<3; i++) {
51: v = -2; j=i+1;
52: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
53: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
54: }
55: i=0; j=3; v=2;
57: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
58: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
59: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
60: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
61: VecAssemblyBegin(B);
62: VecAssemblyEnd(B);
63: PetscPrintf(PETSC_COMM_WORLD,"\nThe Kershaw matrix:\n\n");
64: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
66: /*
67: * A Conjugate Gradient method
68: * with ILU(0) preconditioning
69: */
70: KSPCreate(comm,&ksp);
71: KSPSetOperators(ksp,A,A);
73: KSPSetType(ksp,KSPCG);
74: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
76: /*
77: * ILU preconditioner;
78: * The iterative method will break down unless you comment in the SetShift
79: * line below, or use the -pc_factor_shift_positive_definite option.
80: * Run the code twice: once as given to see the negative pivot and the
81: * divergence behaviour, then comment in the Shift line, or add the
82: * command line option, and see that the pivots are all positive and
83: * the method converges.
84: */
85: KSPGetPC(ksp,&pc);
86: PCSetType(pc,PCICC);
87: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
89: KSPSetFromOptions(ksp);
90: KSPSetUp(ksp);
92: /*
93: * Now that the factorisation is done, show the pivots;
94: * note that the last one is negative. This in itself is not an error,
95: * but it will make the iterative method diverge.
96: */
97: PCFactorGetMatrix(pc,&M);
98: VecDuplicate(B,&D);
99: MatGetDiagonal(M,D);
100: PetscPrintf(PETSC_COMM_WORLD,"\nPivots:\n\n");
101: VecView(D,0);
103: /*
104: * Solve the system;
105: * without the shift this will diverge with
106: * an indefinite preconditioner
107: */
108: KSPSolve(ksp,B,X);
109: KSPGetConvergedReason(ksp,&reason);
110: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
111: PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
112: PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with -pc_factor_shift_positive_definite option.\n");
113: } else if (reason<0) {
114: PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
115: } else {
116: KSPGetIterationNumber(ksp,&its);
117: PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n",(int)its);
118: }
119: PetscPrintf(PETSC_COMM_WORLD,"\n");
121: KSPDestroy(&ksp);
122: MatDestroy(&A);
123: VecDestroy(&B);
124: VecDestroy(&X);
125: VecDestroy(&D);
126: PetscFinalize();
127: return 0;
128: }
130: /*TEST
132: test:
133: filter: sed -e "s/in 5 iterations/in 4 iterations/g"
135: TEST*/