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 */
17
18 #include <petscksp.h>
19
main(int argc,char ** argv)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
31 PetscFunctionBeginUser;
32 PetscCall(PetscInitialize(&argc, &argv, 0, help));
33 comm = MPI_COMM_SELF;
34
35 /*
36 * Construct the Kershaw matrix
37 * and a suitable rhs / initial guess
38 */
39 PetscCall(MatCreateSeqAIJ(comm, 4, 4, 4, 0, &A));
40 PetscCall(VecCreateSeq(comm, 4, &B));
41 PetscCall(VecDuplicate(B, &X));
42 for (i = 0; i < 4; i++) {
43 v = 3;
44 PetscCall(MatSetValues(A, 1, &i, 1, &i, &v, INSERT_VALUES));
45 v = 1;
46 PetscCall(VecSetValues(B, 1, &i, &v, INSERT_VALUES));
47 PetscCall(VecSetValues(X, 1, &i, &v, INSERT_VALUES));
48 }
49
50 i = 0;
51 v = 0;
52 PetscCall(VecSetValues(X, 1, &i, &v, INSERT_VALUES));
53
54 for (i = 0; i < 3; i++) {
55 v = -2;
56 j = i + 1;
57 PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
58 PetscCall(MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES));
59 }
60 i = 0;
61 j = 3;
62 v = 2;
63
64 PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
65 PetscCall(MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES));
66 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
67 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
68 PetscCall(VecAssemblyBegin(B));
69 PetscCall(VecAssemblyEnd(B));
70
71 /*
72 * A Conjugate Gradient method
73 * with ILU(0) preconditioning
74 */
75 PetscCall(KSPCreate(comm, &solver));
76 PetscCall(KSPSetOperators(solver, A, A));
77
78 PetscCall(KSPSetType(solver, KSPCG));
79 PetscCall(KSPSetInitialGuessNonzero(solver, PETSC_TRUE));
80
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 PetscCall(KSPGetPC(solver, &prec));
86 PetscCall(PCSetType(prec, PCILU));
87 /* PetscCall(PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE)); */
88
89 PetscCall(KSPSetFromOptions(solver));
90 PetscCall(KSPSetUp(solver));
91
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 PetscCall(PCFactorGetMatrix(prec, &M));
98 PetscCall(VecDuplicate(B, &D));
99 PetscCall(MatGetDiagonal(M, D));
100
101 /*
102 * Solve the system;
103 * without the shift this will diverge with
104 * an indefinite preconditioner
105 */
106 PetscCall(KSPSolve(solver, B, X));
107 PetscCall(KSPGetConvergedReason(solver, &reason));
108 if (reason == KSP_DIVERGED_INDEFINITE_PC) {
109 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n"));
110 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n"));
111 } else if (reason < 0) {
112 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n"));
113 } else {
114 PetscCall(KSPGetIterationNumber(solver, &its));
115 }
116
117 PetscCall(VecDestroy(&X));
118 PetscCall(VecDestroy(&B));
119 PetscCall(VecDestroy(&D));
120 PetscCall(MatDestroy(&A));
121 PetscCall(KSPDestroy(&solver));
122 PetscCall(PetscFinalize());
123 return 0;
124 }
125
126 /*TEST
127
128 test:
129 args: -pc_factor_shift_type positive_definite
130 output_file: output/empty.out
131
132 TEST*/
133