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