1: /*
  2:  * Test file for the PCFactorSetShiftPd routine or -pc_factor_shift_positive_definite option.
  3:  * The test matrix is the example from Kershaw's paper [J.Comp.Phys 1978]
  4:  * of a positive definite matrix for which ILU(0) will give a negative pivot.
  5:  * This means that the CG method will break down; the Manteuffel shift
  6:  * [Math. Comp. 1980] repairs this.
  7:  *
  8:  * Run the executable twice:
  9:  * 1/ without options: the iterative method diverges because of an
 10:  *    indefinite preconditioner
 11:  * 2/ with -pc_factor_shift_positive_definite option (or comment in the PCFactorSetShiftPd line below):
 12:  *    the method will now successfully converge.
 13:  *
 14:  * Contributed by Victor Eijkhout 2003.
 15:  */

 17: #include <stdlib.h>
 18:  #include petscksp.h

 22: int main(int argc,char **argv)
 23: {
 24:   KSP                solver;
 25:   PC                 prec;
 26:   Mat                A,M;
 27:   Vec                X,B,D;
 28:   MPI_Comm           comm;
 29:   PetscScalar        v;
 30:   KSPConvergedReason reason;
 31:   PetscInt           i,j,its;
 32:   PetscErrorCode     ierr;

 35:   PetscInitialize(&argc,&argv,0,0);
 36:   PetscOptionsSetValue("-options_left",PETSC_NULL);
 37:   comm = MPI_COMM_SELF;
 38: 
 39:   /*
 40:    * Construct the Kershaw matrix
 41:    * and a suitable rhs / initial guess
 42:    */
 43:   MatCreateSeqAIJ(comm,4,4,4,0,&A);
 44:   VecCreateSeq(comm,4,&B);
 45:   VecDuplicate(B,&X);
 46:   for (i=0; i<4; i++) {
 47:     v=3;
 48:     MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);
 49:     v=1;
 50:     VecSetValues(B,1,&i,&v,INSERT_VALUES);
 51:     VecSetValues(X,1,&i,&v,INSERT_VALUES);
 52:   }

 54:   i=0; v=0;
 55:   VecSetValues(X,1,&i,&v,INSERT_VALUES);

 57:   for (i=0; i<3; i++) {
 58:     v=-2; j=i+1;
 59:     MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
 60:     MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
 61:   }
 62:   i=0; j=3; v=2;
 63:   MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
 64:   MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
 65:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 67:   VecAssemblyBegin(B);
 68:   VecAssemblyEnd(B);
 69:   printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);

 71:   /*
 72:    * A Conjugate Gradient method
 73:    * with ILU(0) preconditioning
 74:    */
 75:   KSPCreate(comm,&solver);
 76:   KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN);

 78:   KSPSetType(solver,KSPCG);
 79:   KSPSetInitialGuessNonzero(solver,PETSC_TRUE);

 81:   /*
 82:    * ILU preconditioner;
 83:    * this will break down unless you add the Shift line,
 84:    * or use the -pc_factor_shift_positive_definite option */
 85:   KSPGetPC(solver,&prec);
 86:   PCSetType(prec,PCILU);
 87:   /*  PCFactorSetShiftPd(prec,PETSC_TRUE); */

 89:   KSPSetFromOptions(solver);
 90:   KSPSetUp(solver);

 92:   /*
 93:    * Now that the factorisation is done, show the pivots;
 94:    * note that the last one is negative. This in itself is not an error,
 95:    * but it will make the iterative method diverge.
 96:    */
 97:   PCGetFactoredMatrix(prec,&M);
 98:   VecDuplicate(B,&D);
 99:   MatGetDiagonal(M,D);
100:   printf("\nPivots:\n\n"); VecView(D,0);

102:   /*
103:    * Solve the system;
104:    * without the shift this will diverge with
105:    * an indefinite preconditioner
106:    */
107:   KSPSolve(solver,B,X);
108:   KSPGetConvergedReason(solver,&reason);
109:   if (reason==KSP_DIVERGED_INDEFINITE_PC) {
110:     printf("\nDivergence because of indefinite preconditioner;\n");
111:     printf("Run the executable again but with -pc_factor_shift_positive_definite option.\n");
112:   } else if (reason<0) {
113:     printf("\nOther kind of divergence: this should not happen.\n");
114:   } else {
115:     KSPGetIterationNumber(solver,&its);
116:     printf("\nConvergence in %d iterations.\n",(int)its);
117:   }
118:   printf("\n");

120:   VecDestroy(X);
121:   VecDestroy(B);
122:   VecDestroy(D);
123:   MatDestroy(A);
124:   KSPDestroy(solver);
125:   PetscFinalize();
126:   return(0);
127: }