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 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 29 PetscFunctionBegin; 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 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nThe Kershaw matrix:\n\n")); 70 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 71 72 /* 73 * A Conjugate Gradient method 74 * with ILU(0) preconditioning 75 */ 76 PetscCall(KSPCreate(comm, &ksp)); 77 PetscCall(KSPSetOperators(ksp, A, A)); 78 79 PetscCall(KSPSetType(ksp, KSPCG)); 80 PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); 81 82 /* 83 * ILU preconditioner; 84 * The iterative method will break down unless you comment in the SetShift 85 * line below, or use the -pc_factor_shift_positive_definite option. 86 * Run the code twice: once as given to see the negative pivot and the 87 * divergence behaviour, then comment in the Shift line, or add the 88 * command line option, and see that the pivots are all positive and 89 * the method converges. 90 */ 91 PetscCall(KSPGetPC(ksp, &pc)); 92 PetscCall(PCSetType(pc, PCICC)); 93 /* PetscCall(PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE)); */ 94 95 PetscCall(KSPSetFromOptions(ksp)); 96 PetscCall(KSPSetUp(ksp)); 97 98 /* 99 * Now that the factorisation is done, show the pivots; 100 * note that the last one is negative. This in itself is not an error, 101 * but it will make the iterative method diverge. 102 */ 103 PetscCall(PCFactorGetMatrix(pc, &M)); 104 PetscCall(VecDuplicate(B, &D)); 105 PetscCall(MatGetDiagonal(M, D)); 106 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nPivots:\n\n")); 107 PetscCall(VecView(D, 0)); 108 109 /* 110 * Solve the system; 111 * without the shift this will diverge with 112 * an indefinite preconditioner 113 */ 114 PetscCall(KSPSolve(ksp, B, X)); 115 PetscCall(KSPGetConvergedReason(ksp, &reason)); 116 if (reason == KSP_DIVERGED_INDEFINITE_PC) { 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n")); 118 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with -pc_factor_shift_positive_definite option.\n")); 119 } else if (reason < 0) { 120 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n")); 121 } else { 122 PetscCall(KSPGetIterationNumber(ksp, &its)); 123 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nConvergence in %d iterations.\n", (int)its)); 124 } 125 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n")); 126 127 PetscCall(KSPDestroy(&ksp)); 128 PetscCall(MatDestroy(&A)); 129 PetscCall(VecDestroy(&B)); 130 PetscCall(VecDestroy(&X)); 131 PetscCall(VecDestroy(&D)); 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 138 test: 139 filter: sed -e "s/in 5 iterations/in 4 iterations/g" 140 141 TEST*/ 142