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