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