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