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 petscsles.h
20: int main(int argc,char **argv)
21: {
22: SLES solver; KSP itmeth; PC prec; Mat A,M; Vec X,B,D;
23: MPI_Comm comm;
24: PetscScalar v; KSPConvergedReason reason;
25: int i,j,its,ierr;
28: PetscInitialize(&argc,&argv,0,0);
29: PetscOptionsSetValue("-options_left",PETSC_NULL);
30: comm = MPI_COMM_SELF;
31:
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;
56: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
57: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
58: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
60: VecAssemblyBegin(B);
61: VecAssemblyEnd(B);
62: printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);
64: /*
65: * A Conjugate Gradient method
66: * with ILU(0) preconditioning
67: */
68: SLESCreate(comm,&solver);
69: SLESSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN);
71: SLESGetKSP(solver,&itmeth);
72: KSPSetType(itmeth,KSPCG);
73: KSPSetInitialGuessNonzero(itmeth,PETSC_TRUE);
75: /*
76: * ILU preconditioner;
77: * The iterative method will break down unless you comment in the SetShift
78: * line below, or use the -pc_ilu_shift option.
79: * Run the code twice: once as given to see the negative pivot and the
80: * divergence behaviour, then comment in the Shift line, or add the
81: * command line option, and see that the pivots are all positive and
82: * the method converges.
83: */
84: SLESGetPC(solver,&prec);
85: PCSetType(prec,PCICC);
86: /* PCICCSetShift(prec,PETSC_TRUE); */
88: SLESSetFromOptions(solver);
89: SLESSetUp(solver,B,X);
91: /*
92: * Now that the factorisation is done, show the pivots;
93: * note that the last one is negative. This in itself is not an error,
94: * but it will make the iterative method diverge.
95: */
96: PCGetFactoredMatrix(prec,&M);
97: VecDuplicate(B,&D);
98: MatGetDiagonal(M,D);
99: printf("\nPivots:\n\n"); VecView(D,0);
101: /*
102: * Solve the system;
103: * without the shift this will diverge with
104: * an indefinite preconditioner
105: */
106: SLESSolve(solver,B,X,&its);
107: KSPGetConvergedReason(itmeth,&reason);
108: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
109: printf("\nDivergence because of indefinite preconditioner;\n");
110: printf("Run the executable again but with -pc_icc_shift option.\n");
111: } else if (reason<0) {
112: printf("\nOther kind of divergence: this should not happen.\n");
113: } else {
114: printf("\nConvergence in %d iterations.\n",its);
115: }
116: printf("\n");
118: PetscFinalize();
119: return(0);
120: }