Actual source code: ex1.c
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;
32: PetscInitialize(&argc, &argv, 0, help);
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;
51: v = 0;
52: VecSetValues(X, 1, &i, &v, INSERT_VALUES);
54: for (i = 0; i < 3; i++) {
55: v = -2;
56: 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;
61: j = 3;
62: v = 2;
64: MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES);
65: MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES);
66: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
68: VecAssemblyBegin(B);
69: VecAssemblyEnd(B);
71: /*
72: * A Conjugate Gradient method
73: * with ILU(0) preconditioning
74: */
75: KSPCreate(comm, &solver);
76: KSPSetOperators(solver, A, A);
78: KSPSetType(solver, KSPCG);
79: KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
81: /*
82: * ILU preconditioner;
83: * this will break down unless you add the Shift line,
84: * or use the -pc_factor_shift_positive_definite option */
85: KSPGetPC(solver, &prec);
86: PCSetType(prec, PCILU);
87: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
89: KSPSetFromOptions(solver);
90: KSPSetUp(solver);
92: /*
93: * Now that the factorisation is done, show the pivots;
94: * note that the last one is negative. This in itself is not an error,
95: * but it will make the iterative method diverge.
96: */
97: PCFactorGetMatrix(prec, &M);
98: VecDuplicate(B, &D);
99: MatGetDiagonal(M, D);
101: /*
102: * Solve the system;
103: * without the shift this will diverge with
104: * an indefinite preconditioner
105: */
106: KSPSolve(solver, B, X);
107: KSPGetConvergedReason(solver, &reason);
108: if (reason == KSP_DIVERGED_INDEFINITE_PC) {
109: PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n");
110: PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");
111: } else if (reason < 0) {
112: PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n");
113: } else {
114: KSPGetIterationNumber(solver, &its);
115: }
117: VecDestroy(&X);
118: VecDestroy(&B);
119: VecDestroy(&D);
120: MatDestroy(&A);
121: KSPDestroy(&solver);
122: PetscFinalize();
123: return 0;
124: }
126: /*TEST
128: test:
129: args: -pc_factor_shift_type positive_definite
131: TEST*/