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 KSP solver; 22 PC prec; 23 Mat A, M; 24 Vec X, B, D; 25 MPI_Comm comm; 26 PetscScalar v; 27 KSPConvergedReason reason; 28 PetscInt i, j, its; 29 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 70 /* 71 * A Conjugate Gradient method 72 * with ILU(0) preconditioning 73 */ 74 PetscCall(KSPCreate(comm, &solver)); 75 PetscCall(KSPSetOperators(solver, A, A)); 76 77 PetscCall(KSPSetType(solver, KSPCG)); 78 PetscCall(KSPSetInitialGuessNonzero(solver, PETSC_TRUE)); 79 80 /* 81 * ILU preconditioner; 82 * this will break down unless you add the Shift line, 83 * or use the -pc_factor_shift_positive_definite option */ 84 PetscCall(KSPGetPC(solver, &prec)); 85 PetscCall(PCSetType(prec, PCILU)); 86 /* PetscCall(PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE)); */ 87 88 PetscCall(KSPSetFromOptions(solver)); 89 PetscCall(KSPSetUp(solver)); 90 91 /* 92 * Now that the factorisation is done, show the pivots; 93 * note that the last one is negative. This in itself is not an error, 94 * but it will make the iterative method diverge. 95 */ 96 PetscCall(PCFactorGetMatrix(prec, &M)); 97 PetscCall(VecDuplicate(B, &D)); 98 PetscCall(MatGetDiagonal(M, D)); 99 100 /* 101 * Solve the system; 102 * without the shift this will diverge with 103 * an indefinite preconditioner 104 */ 105 PetscCall(KSPSolve(solver, B, X)); 106 PetscCall(KSPGetConvergedReason(solver, &reason)); 107 if (reason == KSP_DIVERGED_INDEFINITE_PC) { 108 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n")); 109 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n")); 110 } else if (reason < 0) { 111 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n")); 112 } else { 113 PetscCall(KSPGetIterationNumber(solver, &its)); 114 } 115 116 PetscCall(VecDestroy(&X)); 117 PetscCall(VecDestroy(&B)); 118 PetscCall(VecDestroy(&D)); 119 PetscCall(MatDestroy(&A)); 120 PetscCall(KSPDestroy(&solver)); 121 PetscCall(PetscFinalize()); 122 return 0; 123 } 124 125 /*TEST 126 127 test: 128 args: -pc_factor_shift_type positive_definite 129 130 TEST*/ 131