Actual source code: ex1.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: *
14: * Contributed by Victor Eijkhout 2003.
15: */
17: #include <stdlib.h>
18: #include petscsles.h
22: int main(int argc,char **argv)
23: {
24: SLES solver; KSP itmeth; PC prec; Mat A,M; Vec X,B,D;
25: MPI_Comm comm;
26: PetscScalar v; KSPConvergedReason reason;
27: int i,j,its,ierr;
30: PetscInitialize(&argc,&argv,0,0);
31: PetscOptionsSetValue("-options_left",PETSC_NULL);
32: comm = MPI_COMM_SELF;
33:
34: /*
35: * Construct the Kershaw matrix
36: * and a suitable rhs / initial guess
37: */
38: MatCreateSeqAIJ(comm,4,4,4,0,&A);
39: VecCreateSeq(comm,4,&B);
40: VecDuplicate(B,&X);
41: for (i=0; i<4; i++) {
42: v=3;
43: MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
44: v=1;
45: VecSetValues(B,1,&i,&v,INSERT_VALUES);
46: VecSetValues(X,1,&i,&v,INSERT_VALUES);
47: }
49: i=0; v=0;
50: VecSetValues(X,1,&i,&v,INSERT_VALUES);
52: for (i=0; i<3; i++) {
53: v=-2; j=i+1;
54: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
55: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
56: }
57: i=0; j=3; v=2;
58: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
59: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
60: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
61: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
62: VecAssemblyBegin(B);
63: VecAssemblyEnd(B);
64: printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);
66: /*
67: * A Conjugate Gradient method
68: * with ILU(0) preconditioning
69: */
70: SLESCreate(comm,&solver);
71: SLESSetOperators(solver,A,A,SAME_NONZERO_PATTERN);
73: SLESGetKSP(solver,&itmeth);
74: KSPSetType(itmeth,KSPCG);
75: KSPSetInitialGuessNonzero(itmeth,PETSC_TRUE);
77: /*
78: * ILU preconditioner;
79: * this will break down unless you add the Shift line,
80: * or use the -pc_ilu_shift option */
81: SLESGetPC(solver,&prec);
82: PCSetType(prec,PCILU);
83: /* PCILUSetShift(prec,PETSC_TRUE); */
85: SLESSetFromOptions(solver);
86: SLESSetUp(solver,B,X);
88: /*
89: * Now that the factorisation is done, show the pivots;
90: * note that the last one is negative. This in itself is not an error,
91: * but it will make the iterative method diverge.
92: */
93: PCGetFactoredMatrix(prec,&M);
94: VecDuplicate(B,&D);
95: MatGetDiagonal(M,D);
96: printf("\nPivots:\n\n"); VecView(D,0);
98: /*
99: * Solve the system;
100: * without the shift this will diverge with
101: * an indefinite preconditioner
102: */
103: SLESSolve(solver,B,X,&its);
104: KSPGetConvergedReason(itmeth,&reason);
105: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
106: printf("\nDivergence because of indefinite preconditioner;\n");
107: printf("Run the executable again but with -pc_ilu_shift option.\n");
108: } else if (reason<0) {
109: printf("\nOther kind of divergence: this should not happen.\n");
110: } else {
111: printf("\nConvergence in %d iterations.\n",its);
112: }
113: printf("\n");
115: PetscFinalize();
116: return(0);
117: }