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