Actual source code: ex2.c

      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: #include <petscksp.h>

     18: int main(int argc,char **argv)
     19: {
     20:   KSP                ksp;
     21:   PC                 pc;
     22:   Mat                A,M;
     23:   Vec                X,B,D;
     24:   MPI_Comm           comm;
     25:   PetscScalar        v;
     26:   KSPConvergedReason reason;
     27:   PetscInt           i,j,its;
     28:   PetscErrorCode     ierr;

     31:   PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
     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);
     65:   PetscPrintf(PETSC_COMM_WORLD,"\nThe Kershaw matrix:\n\n");
     66:   MatView(A,PETSC_VIEWER_STDOUT_WORLD);

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

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

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

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

     94:   /*
     95:    * Now that the factorisation is done, show the pivots;
     96:    * note that the last one is negative. This in itself is not an error,
     97:    * but it will make the iterative method diverge.
     98:    */
     99:   PCFactorGetMatrix(pc,&M);
    100:   VecDuplicate(B,&D);
    101:   MatGetDiagonal(M,D);
    102:   PetscPrintf(PETSC_COMM_WORLD,"\nPivots:\n\n");
    103:   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:     PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
    114:     PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with -pc_factor_shift_positive_definite option.\n");
    115:   } else if (reason<0) {
    116:     PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
    117:   } else {
    118:     KSPGetIterationNumber(ksp,&its);
    119:     PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n",(int)its);
    120:   }
    121:   PetscPrintf(PETSC_COMM_WORLD,"\n");

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

    133: /*TEST

    135:    test:
    136:      filter:  sed -e "s/in 5 iterations/in 4 iterations/g"

    139: TEST*/