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