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