Actual source code: ex2.c
petsc-3.6.1 2015-08-06
1: static char help[] = "Test file for the PCFactorSetShiftType()\n";
2: /*
3: * Test file for the PCFactorSetShiftType() routine or -pc_factor_shift_type POSITIVE_DEFINITE option.
4: * The test matrix is the example from Kershaw's paper [J.Comp.Phys 1978]
5: * of a positive definite matrix for which ILU(0) will give a negative pivot.
6: * This means that the CG method will break down; the Manteuffel shift
7: * [Math. Comp. 1980] repairs this.
8: *
9: * Run the executable twice:
10: * 1/ without options: the iterative method diverges because of an
11: * indefinite preconditioner
12: * 2/ with -pc_factor_shift_positive_definite option (or comment in the PCFactorSetShiftType() line below):
13: * the method will now successfully converge.
14: */
16: #include <petscksp.h>
20: int main(int argc,char **argv)
21: {
22: KSP ksp;
23: PC pc;
24: Mat A,M;
25: Vec X,B,D;
26: MPI_Comm comm;
27: PetscScalar v;
28: KSPConvergedReason reason;
29: PetscInt i,j,its;
30: PetscErrorCode ierr;
33: PetscInitialize(&argc,&argv,0,help);
34: PetscOptionsSetValue("-options_left",NULL);
35: comm = MPI_COMM_SELF;
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;
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,&ksp);
75: KSPSetOperators(ksp,A,A);
77: KSPSetType(ksp,KSPCG);
78: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
80: /*
81: * ILU preconditioner;
82: * The iterative method will break down unless you comment in the SetShift
83: * line below, or use the -pc_factor_shift_positive_definite option.
84: * Run the code twice: once as given to see the negative pivot and the
85: * divergence behaviour, then comment in the Shift line, or add the
86: * command line option, and see that the pivots are all positive and
87: * the method converges.
88: */
89: KSPGetPC(ksp,&pc);
90: PCSetType(pc,PCICC);
91: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
93: KSPSetFromOptions(ksp);
94: KSPSetUp(ksp);
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: PCFactorGetMatrix(pc,&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(ksp,B,X);
112: KSPGetConvergedReason(ksp,&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_factor_shift_positive_definite option.\n");
116: } else if (reason<0) {
117: printf("\nOther kind of divergence: this should not happen.\n");
118: } else {
119: KSPGetIterationNumber(ksp,&its);
120: printf("\nConvergence in %d iterations.\n",(int)its);
121: }
122: printf("\n");
124: KSPDestroy(&ksp);
125: MatDestroy(&A);
126: VecDestroy(&B);
127: VecDestroy(&X);
128: VecDestroy(&D);
129: PetscFinalize();
130: return(0);
131: }