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 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