Actual source code: ex3.c
1: /*
2: * Test file for the PCFactorSetShiftType() 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_positive_definite option (or comment in the PCFactorSetShiftType() line below):
12: * the method will now successfully converge.
13: *
14: * Modified from ex1.c by malte.foerster@scai.fraunhofer.de [petsc-maint #42323]
15: * such that the matrix A has inode structure.
16: */
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: PetscInt i,j,its;
32: PetscErrorCode ierr;
33: PetscInt nnu=1000;
34:
35: PetscInitialize(&argc,&argv,0,0);
36: //PetscOptionsSetValue("-options_left",PETSC_NULL);
37: comm = MPI_COMM_SELF;
38:
39:
40: /*
41: * Construct the Kershaw matrix
42: * and a suitable rhs / initial guess
43: */
44: MatCreateSeqAIJ(comm,nnu,nnu,20,0,&A);
45: VecCreateSeq(comm,nnu,&B);
46: VecDuplicate(B,&X);
47: for (i=0; i<nnu; i++) {
48: v=3;
49: MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
50: v=1;
51: VecSetValues(B,1,&i,&v,INSERT_VALUES);
52: VecSetValues(X,1,&i,&v,INSERT_VALUES);
53: }
55: i=0; v=0;
56: VecSetValues(X,1,&i,&v,INSERT_VALUES);
58: for (i=0; i<nnu-1; i+=1) {
59: v=-2; j=i+1;
60: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
61: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
62: if (i>4) i++;
63: }
65: i=0; j=3; v=2;
66: MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
67: MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
70: VecAssemblyBegin(B);
71: VecAssemblyEnd(B);
72: //printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);
74: /*
75: * A Conjugate Gradient method
76: * with ILU(0) preconditioning
77: */
78: KSPCreate(comm,&solver);
79: KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN);
81: KSPSetType(solver,KSPCG);
82: KSPSetInitialGuessNonzero(solver,PETSC_TRUE);
84: /*
85: * ILU preconditioner;
86: * this will break down unless you add the Shift line,
87: * or use the -pc_factor_shift_positive_definite option */
88: KSPGetPC(solver,&prec);
89: PCSetType(prec,PCILU);
90: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
92: KSPSetFromOptions(solver);
93: KSPSetUp(solver);
95: /*
96: * Now that the factorisation is done, show the pivots;
97: * note that the last one is negative. This in itself is not an error,
98: * but it will make the iterative method diverge.
99: */
100: PCFactorGetMatrix(prec,&M);
101: VecDuplicate(B,&D);
102: MatGetDiagonal(M,D);
103: //printf("\nPivots:\n\n"); VecView(D,0);
105: /*
106: * Solve the system;
107: * without the shift this will diverge with
108: * an indefinite preconditioner
109: */
110: KSPSolve(solver,B,X);
111: KSPGetConvergedReason(solver,&reason);
112: if (reason==KSP_DIVERGED_INDEFINITE_PC) {
113: printf("\nDivergence because of indefinite preconditioner;\n");
114: printf("Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");
115: } else if (reason<0) {
116: printf("\nOther kind of divergence: this should not happen.\n");
117: } else {
118: KSPGetIterationNumber(solver,&its);
119: printf("\nConvergence in %d iterations.\n",(int)its);
120: }
121: printf("\n");
123: VecDestroy(X);
124: VecDestroy(B);
125: VecDestroy(D);
126: MatDestroy(A);
127: KSPDestroy(solver);
128: PetscFinalize();
129: return 0;
130: }