Actual source code: ex1.c
petsc-3.11.4 2019-09-28
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_type POSITIVE_DEFINITE option (or comment in the PCFactorSetShiftType() line below):
13: * the method will now successfully converge.
14: *
15: * Contributed by Victor Eijkhout 2003.
16: */
18: #include <petscksp.h>
20: int main(int argc,char **argv)
21: {
22: KSP solver;
23: PC prec;
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;
32: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
33: comm = MPI_COMM_SELF;
35: /*
36: * Construct the Kershaw matrix
37: * and a suitable rhs / initial guess
38: */
39: MatCreateSeqAIJ(comm,4,4,4,0,&A);
40: VecCreateSeq(comm,4,&B);
41: VecDuplicate(B,&X);
42: for (i=0; i<4; i++) {
43: v = 3;
44: MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
45: v = 1;
46: VecSetValues(B,1,&i,&v,INSERT_VALUES);
47: VecSetValues(X,1,&i,&v,INSERT_VALUES);
48: }
50: i=0; v=0;
51: VecSetValues(X,1,&i,&v,INSERT_VALUES);
53: for (i=0; i<3; i++) {
54: v = -2; j=i+1;
55: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
56: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
57: }
58: i=0; j=3; v=2;
60: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
61: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
62: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
63: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
64: VecAssemblyBegin(B);
65: VecAssemblyEnd(B);
67: /*
68: * A Conjugate Gradient method
69: * with ILU(0) preconditioning
70: */
71: KSPCreate(comm,&solver);
72: KSPSetOperators(solver,A,A);
74: KSPSetType(solver,KSPCG);
75: KSPSetInitialGuessNonzero(solver,PETSC_TRUE);
77: /*
78: * ILU preconditioner;
79: * this will break down unless you add the Shift line,
80: * or use the -pc_factor_shift_positive_definite option */
81: KSPGetPC(solver,&prec);
82: PCSetType(prec,PCILU);
83: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
85: KSPSetFromOptions(solver);
86: KSPSetUp(solver);
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: PCFactorGetMatrix(prec,&M);
94: VecDuplicate(B,&D);
95: MatGetDiagonal(M,D);
97: /*
98: * Solve the system;
99: * without the shift this will diverge with
100: * an indefinite preconditioner
101: */
102: KSPSolve(solver,B,X);
103: KSPGetConvergedReason(solver,&reason);
104: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
105: PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
106: PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");
107: } else if (reason<0) {
108: PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
109: } else {
110: KSPGetIterationNumber(solver,&its);
111: }
113: VecDestroy(&X);
114: VecDestroy(&B);
115: VecDestroy(&D);
116: MatDestroy(&A);
117: KSPDestroy(&solver);
118: PetscFinalize();
119: return ierr;
120: }
123: /*TEST
125: test:
126: args: -pc_factor_shift_type positive_definite
128: TEST*/