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 petscksp.h
22: int main(int argc,char **argv)
23: {
24: KSP solver;
25: PC prec;
26: Mat A,M;
27: Vec X,B,D;
28: MPI_Comm comm;
29: PetscScalar v;
30: KSPConvergedReason reason;
31: int i,j,its,ierr;
34: PetscInitialize(&argc,&argv,0,0);
35: PetscOptionsSetValue("-options_left",PETSC_NULL);
36: comm = MPI_COMM_SELF;
37:
38: /*
39: * Construct the Kershaw matrix
40: * and a suitable rhs / initial guess
41: */
42: MatCreateSeqAIJ(comm,4,4,4,0,&A);
43: VecCreateSeq(comm,4,&B);
44: VecDuplicate(B,&X);
45: for (i=0; i<4; i++) {
46: v=3;
47: MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
48: v=1;
49: VecSetValues(B,1,&i,&v,INSERT_VALUES);
50: VecSetValues(X,1,&i,&v,INSERT_VALUES);
51: }
53: i=0; v=0;
54: VecSetValues(X,1,&i,&v,INSERT_VALUES);
56: for (i=0; i<3; i++) {
57: v=-2; j=i+1;
58: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
59: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
60: }
61: i=0; j=3; v=2;
62: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
63: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
64: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
66: VecAssemblyBegin(B);
67: VecAssemblyEnd(B);
68: printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);
70: /*
71: * A Conjugate Gradient method
72: * with ILU(0) preconditioning
73: */
74: KSPCreate(comm,&solver);
75: KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN);
77: KSPSetType(solver,KSPCG);
78: KSPSetInitialGuessNonzero(solver,PETSC_TRUE);
80: /*
81: * ILU preconditioner;
82: * this will break down unless you add the Shift line,
83: * or use the -pc_ilu_shift option */
84: KSPGetPC(solver,&prec);
85: PCSetType(prec,PCILU);
86: /* PCILUSetShift(prec,PETSC_TRUE); */
88: KSPSetFromOptions(solver);
89: KSPSetRhs(solver,B);
90: KSPSetSolution(solver,X);
91: KSPSetUp(solver);
93: /*
94: * Now that the factorisation is done, show the pivots;
95: * note that the last one is negative. This in itself is not an error,
96: * but it will make the iterative method diverge.
97: */
98: PCGetFactoredMatrix(prec,&M);
99: VecDuplicate(B,&D);
100: MatGetDiagonal(M,D);
101: printf("\nPivots:\n\n"); VecView(D,0);
103: /*
104: * Solve the system;
105: * without the shift this will diverge with
106: * an indefinite preconditioner
107: */
108: KSPSolve(solver);
109: KSPGetConvergedReason(solver,&reason);
110: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
111: printf("\nDivergence because of indefinite preconditioner;\n");
112: printf("Run the executable again but with -pc_ilu_shift option.\n");
113: } else if (reason<0) {
114: printf("\nOther kind of divergence: this should not happen.\n");
115: } else {
116: KSPGetIterationNumber(solver,&its);
117: printf("\nConvergence in %d iterations.\n",its);
118: }
119: printf("\n");
121: VecDestroy(X);
122: VecDestroy(B);
123: VecDestroy(D);
124: MatDestroy(A);
125: KSPDestroy(solver);
126: PetscFinalize();
127: return(0);
128: }