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

 16: /*T
 17:    requires: x
 18: T*/

 20:  #include <petscksp.h>

 22: int main(int argc,char **argv)
 23: {
 24:   KSP                ksp;
 25:   PC                 pc;
 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,help);if (ierr) return ierr;
 36:   comm = MPI_COMM_SELF;

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

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

 56:   for (i=0; i<3; i++) {
 57:     v    = -2; j=i+1;
 58:     MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);
 59:     MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);
 60:   }
 61:   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:   PetscPrintf(PETSC_COMM_WORLD,"\nThe Kershaw matrix:\n\n");
 70:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

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

 79:   KSPSetType(ksp,KSPCG);
 80:   KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);

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

 95:   KSPSetFromOptions(ksp);
 96:   KSPSetUp(ksp);

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

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

127:   KSPDestroy(&ksp);
128:   MatDestroy(&A);
129:   VecDestroy(&B);
130:   VecDestroy(&X);
131:   VecDestroy(&D);
132:   PetscFinalize();
133:   return ierr;
134: }


137: /*TEST

139:    test:

141: TEST*/