Actual source code: ex2.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_positive_definite option (or comment in the PCFactorSetShiftType() line below):
13: * the method will now successfully converge.
14: */
16: #include <petscksp.h>
18: int main(int argc, char **argv)
19: {
20: KSP ksp;
21: PC pc;
22: Mat A, M;
23: Vec X, B, D;
24: MPI_Comm comm;
25: PetscScalar v;
26: KSPConvergedReason reason;
27: PetscInt i, j, its;
30: PetscInitialize(&argc, &argv, 0, help);
31: comm = MPI_COMM_SELF;
33: /*
34: * Construct the Kershaw matrix
35: * and a suitable rhs / initial guess
36: */
37: MatCreateSeqAIJ(comm, 4, 4, 4, 0, &A);
38: VecCreateSeq(comm, 4, &B);
39: VecDuplicate(B, &X);
40: for (i = 0; i < 4; i++) {
41: v = 3;
42: MatSetValues(A, 1, &i, 1, &i, &v, INSERT_VALUES);
43: v = 1;
44: VecSetValues(B, 1, &i, &v, INSERT_VALUES);
45: VecSetValues(X, 1, &i, &v, INSERT_VALUES);
46: }
48: i = 0;
49: v = 0;
50: VecSetValues(X, 1, &i, &v, INSERT_VALUES);
52: for (i = 0; i < 3; i++) {
53: v = -2;
54: 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;
59: j = 3;
60: 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: PetscPrintf(PETSC_COMM_WORLD, "\nThe Kershaw matrix:\n\n");
69: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
71: /*
72: * A Conjugate Gradient method
73: * with ILU(0) preconditioning
74: */
75: KSPCreate(comm, &ksp);
76: KSPSetOperators(ksp, A, A);
78: KSPSetType(ksp, KSPCG);
79: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
81: /*
82: * ILU preconditioner;
83: * The iterative method will break down unless you comment in the SetShift
84: * line below, or use the -pc_factor_shift_positive_definite option.
85: * Run the code twice: once as given to see the negative pivot and the
86: * divergence behaviour, then comment in the Shift line, or add the
87: * command line option, and see that the pivots are all positive and
88: * the method converges.
89: */
90: KSPGetPC(ksp, &pc);
91: PCSetType(pc, PCICC);
92: /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */
94: KSPSetFromOptions(ksp);
95: KSPSetUp(ksp);
97: /*
98: * Now that the factorisation is done, show the pivots;
99: * note that the last one is negative. This in itself is not an error,
100: * but it will make the iterative method diverge.
101: */
102: PCFactorGetMatrix(pc, &M);
103: VecDuplicate(B, &D);
104: MatGetDiagonal(M, D);
105: PetscPrintf(PETSC_COMM_WORLD, "\nPivots:\n\n");
106: VecView(D, 0);
108: /*
109: * Solve the system;
110: * without the shift this will diverge with
111: * an indefinite preconditioner
112: */
113: KSPSolve(ksp, B, X);
114: KSPGetConvergedReason(ksp, &reason);
115: if (reason == KSP_DIVERGED_INDEFINITE_PC) {
116: PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n");
117: PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with -pc_factor_shift_positive_definite option.\n");
118: } else if (reason < 0) {
119: PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n");
120: } else {
121: KSPGetIterationNumber(ksp, &its);
122: PetscPrintf(PETSC_COMM_WORLD, "\nConvergence in %d iterations.\n", (int)its);
123: }
124: PetscPrintf(PETSC_COMM_WORLD, "\n");
126: KSPDestroy(&ksp);
127: MatDestroy(&A);
128: VecDestroy(&B);
129: VecDestroy(&X);
130: VecDestroy(&D);
131: PetscFinalize();
132: return 0;
133: }
135: /*TEST
137: test:
138: filter: sed -e "s/in 5 iterations/in 4 iterations/g"
140: TEST*/