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 PetscErrorCode ierr; 31 32 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 33 comm = MPI_COMM_SELF; 34 35 /* 36 * Construct the Kershaw matrix 37 * and a suitable rhs / initial guess 38 */ 39 ierr = MatCreateSeqAIJ(comm,4,4,4,0,&A);CHKERRQ(ierr); 40 ierr = VecCreateSeq(comm,4,&B);CHKERRQ(ierr); 41 ierr = VecDuplicate(B,&X);CHKERRQ(ierr); 42 for (i=0; i<4; i++) { 43 v = 3; 44 ierr = MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 45 v = 1; 46 ierr = VecSetValues(B,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 47 ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 48 } 49 50 i=0; v=0; 51 ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 52 53 for (i=0; i<3; i++) { 54 v = -2; j=i+1; 55 ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 56 ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 57 } 58 i=0; j=3; v=2; 59 60 ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr); 61 ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 62 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 63 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 64 ierr = VecAssemblyBegin(B);CHKERRQ(ierr); 65 ierr = VecAssemblyEnd(B);CHKERRQ(ierr); 66 67 /* 68 * A Conjugate Gradient method 69 * with ILU(0) preconditioning 70 */ 71 ierr = KSPCreate(comm,&solver);CHKERRQ(ierr); 72 ierr = KSPSetOperators(solver,A,A);CHKERRQ(ierr); 73 74 ierr = KSPSetType(solver,KSPCG);CHKERRQ(ierr); 75 ierr = KSPSetInitialGuessNonzero(solver,PETSC_TRUE);CHKERRQ(ierr); 76 77 /* 78 * ILU preconditioner; 79 * this will break down unless you add the Shift line, 80 * or use the -pc_factor_shift_positive_definite option */ 81 ierr = KSPGetPC(solver,&prec);CHKERRQ(ierr); 82 ierr = PCSetType(prec,PCILU);CHKERRQ(ierr); 83 /* ierr = PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE);CHKERRQ(ierr); */ 84 85 ierr = KSPSetFromOptions(solver);CHKERRQ(ierr); 86 ierr = KSPSetUp(solver);CHKERRQ(ierr); 87 88 /* 89 * Now that the factorisation is done, show the pivots; 90 * note that the last one is negative. This in itself is not an error, 91 * but it will make the iterative method diverge. 92 */ 93 ierr = PCFactorGetMatrix(prec,&M);CHKERRQ(ierr); 94 ierr = VecDuplicate(B,&D);CHKERRQ(ierr); 95 ierr = MatGetDiagonal(M,D);CHKERRQ(ierr); 96 97 /* 98 * Solve the system; 99 * without the shift this will diverge with 100 * an indefinite preconditioner 101 */ 102 ierr = KSPSolve(solver,B,X);CHKERRQ(ierr); 103 ierr = KSPGetConvergedReason(solver,&reason);CHKERRQ(ierr); 104 if (reason==KSP_DIVERGED_INDEFINITE_PC) { 105 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");CHKERRQ(ierr); 106 ierr = PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");CHKERRQ(ierr); 107 } else if (reason<0) { 108 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");CHKERRQ(ierr); 109 } else { 110 ierr = KSPGetIterationNumber(solver,&its);CHKERRQ(ierr); 111 } 112 113 ierr = VecDestroy(&X);CHKERRQ(ierr); 114 ierr = VecDestroy(&B);CHKERRQ(ierr); 115 ierr = VecDestroy(&D);CHKERRQ(ierr); 116 ierr = MatDestroy(&A);CHKERRQ(ierr); 117 ierr = KSPDestroy(&solver);CHKERRQ(ierr); 118 ierr = PetscFinalize(); 119 return ierr; 120 } 121 122 123 /*TEST 124 125 test: 126 args: -pc_factor_shift_type positive_definite 127 128 TEST*/ 129