Actual source code: ex2.c

  1: /*
  2:  * Test file for the PCILUSetShift routine or -pc_shift 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_ilu_shift option (or comment in the PCILUSetShift line below):
 12:  *    the method will now successfully converge.
 13:  */

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

 20: int main(int argc,char **argv)
 21: {
 22:   SLES solver; KSP itmeth; PC prec; Mat A,M; Vec X,B,D;
 23:   MPI_Comm comm;
 24:   PetscScalar v; KSPConvergedReason reason;
 25:   int i,j,its,ierr;

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

 47:   i=0; v=0;
 48:   VecSetValues(X,1,&i,&v,INSERT_VALUES);

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

 64:   /*
 65:    * A Conjugate Gradient method
 66:    * with ILU(0) preconditioning
 67:    */
 68:   SLESCreate(comm,&solver);
 69:   SLESSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN);

 71:   SLESGetKSP(solver,&itmeth);
 72:   KSPSetType(itmeth,KSPCG);
 73:   KSPSetInitialGuessNonzero(itmeth,PETSC_TRUE);

 75:   /*
 76:    * ILU preconditioner;
 77:    * The iterative method will break down unless you comment in the SetShift
 78:    * line below, or use the -pc_ilu_shift option.
 79:    * Run the code twice: once as given to see the negative pivot and the
 80:    * divergence behaviour, then comment in the Shift line, or add the 
 81:    * command line option, and see that the pivots are all positive and
 82:    * the method converges.
 83:    */
 84:   SLESGetPC(solver,&prec);
 85:   PCSetType(prec,PCICC);
 86:   /*  PCICCSetShift(prec,PETSC_TRUE);  */

 88:   SLESSetFromOptions(solver);
 89:   SLESSetUp(solver,B,X);

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

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

118:   PetscFinalize();
119:   return(0);
120: }