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:  */

 18: #include <petscksp.h>

 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;

 31:   PetscInitialize(&argc,&argv,0,help);
 32:   comm = MPI_COMM_SELF;

 34:   /*
 35:    * Construct the Kershaw matrix
 36:    * and a suitable rhs / initial guess
 37:    */
 38:   MatCreateSeqAIJ(comm,4,4,4,0,&A);
 39:   VecCreateSeq(comm,4,&B);
 40:   VecDuplicate(B,&X);
 41:   for (i=0; i<4; i++) {
 42:     v    = 3;
 43:     MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
 44:     v    = 1;
 45:     VecSetValues(B,1,&i,&v,INSERT_VALUES);
 46:     VecSetValues(X,1,&i,&v,INSERT_VALUES);
 47:   }

 49:   i=0; v=0;
 50:   VecSetValues(X,1,&i,&v,INSERT_VALUES);

 52:   for (i=0; i<3; i++) {
 53:     v    = -2; j=i+1;
 54:     MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
 55:     MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
 56:   }
 57:   i=0; j=3; v=2;

 59:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
 60:   MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
 61:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 62:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 63:   VecAssemblyBegin(B);
 64:   VecAssemblyEnd(B);

 66:   /*
 67:    * A Conjugate Gradient method
 68:    * with ILU(0) preconditioning
 69:    */
 70:   KSPCreate(comm,&solver);
 71:   KSPSetOperators(solver,A,A);

 73:   KSPSetType(solver,KSPCG);
 74:   KSPSetInitialGuessNonzero(solver,PETSC_TRUE);

 76:   /*
 77:    * ILU preconditioner;
 78:    * this will break down unless you add the Shift line,
 79:    * or use the -pc_factor_shift_positive_definite option */
 80:   KSPGetPC(solver,&prec);
 81:   PCSetType(prec,PCILU);
 82:   /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */

 84:   KSPSetFromOptions(solver);
 85:   KSPSetUp(solver);

 87:   /*
 88:    * Now that the factorisation is done, show the pivots;
 89:    * note that the last one is negative. This in itself is not an error,
 90:    * but it will make the iterative method diverge.
 91:    */
 92:   PCFactorGetMatrix(prec,&M);
 93:   VecDuplicate(B,&D);
 94:   MatGetDiagonal(M,D);

 96:   /*
 97:    * Solve the system;
 98:    * without the shift this will diverge with
 99:    * an indefinite preconditioner
100:    */
101:   KSPSolve(solver,B,X);
102:   KSPGetConvergedReason(solver,&reason);
103:   if (reason==KSP_DIVERGED_INDEFINITE_PC) {
104:     PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
105:     PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");
106:   } else if (reason<0) {
107:     PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
108:   } else {
109:     KSPGetIterationNumber(solver,&its);
110:   }

112:   VecDestroy(&X);
113:   VecDestroy(&B);
114:   VecDestroy(&D);
115:   MatDestroy(&A);
116:   KSPDestroy(&solver);
117:   PetscFinalize();
118:   return 0;
119: }

121: /*TEST

123:    test:
124:       args: -pc_factor_shift_type positive_definite

126: TEST*/