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