1: /*
  2:  * Test file for the PCFactorSetShiftType() routine or -pc_factor_shift_type 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 PCFactorSetShiftType() line below):
 12:  *    the method will now successfully converge.
 13:  */

 15: #include <stdlib.h>
 16: #include <petscksp.h>

 20: int main(int argc,char **argv)
 21: {
 22:   KSP                ksp;
 23:   PC                 pc;
 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;

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

 52:   i=0; v=0;
 53:   VecSetValues(X,1,&i,&v,INSERT_VALUES);

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

 69:   /*
 70:    * A Conjugate Gradient method
 71:    * with ILU(0) preconditioning
 72:    */
 73:   KSPCreate(comm,&ksp);
 74:   KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);

 76:   KSPSetType(ksp,KSPCG);
 77:   KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);

 79:   /*
 80:    * ILU preconditioner;
 81:    * The iterative method will break down unless you comment in the SetShift
 82:    * line below, or use the -pc_factor_shift_positive_definite option.
 83:    * Run the code twice: once as given to see the negative pivot and the
 84:    * divergence behaviour, then comment in the Shift line, or add the 
 85:    * command line option, and see that the pivots are all positive and
 86:    * the method converges.
 87:    */
 88:   KSPGetPC(ksp,&pc);
 89:   PCSetType(pc,PCICC);
 90:   /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */

 92:   KSPSetFromOptions(ksp);
 93:   KSPSetUp(ksp);

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

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

123:   KSPDestroy(&ksp);
124:   MatDestroy(&A);
125:   VecDestroy(&B);
126:   VecDestroy(&X);
127:   VecDestroy(&D);
128:   PetscFinalize();
129:   return(0);
130: }