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;

     29:   PetscFunctionBegin;
     30:   PetscFunctionBeginUser;
     31:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
     32:   comm = MPI_COMM_SELF;

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

     49:   i = 0;
     50:   v = 0;
     51:   PetscCall(VecSetValues(X, 1, &i, &v, INSERT_VALUES));

     53:   for (i = 0; i < 3; i++) {
     54:     v = -2;
     55:     j = i + 1;
     56:     PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
     57:     PetscCall(MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES));
     58:   }
     59:   i = 0;
     60:   j = 3;
     61:   v = 2;

     63:   PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
     64:   PetscCall(MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES));
     65:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     66:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     67:   PetscCall(VecAssemblyBegin(B));
     68:   PetscCall(VecAssemblyEnd(B));
     69:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nThe Kershaw matrix:\n\n"));
     70:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));

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

     79:   PetscCall(KSPSetType(ksp, KSPCG));
     80:   PetscCall(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:   PetscCall(KSPGetPC(ksp, &pc));
     92:   PetscCall(PCSetType(pc, PCICC));
     93:   /* PetscCall(PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE)); */

     95:   PetscCall(KSPSetFromOptions(ksp));
     96:   PetscCall(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:   PetscCall(PCFactorGetMatrix(pc, &M));
    104:   PetscCall(VecDuplicate(B, &D));
    105:   PetscCall(MatGetDiagonal(M, D));
    106:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nPivots:\n\n"));
    107:   PetscCall(VecView(D, 0));

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

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

    136: /*TEST

    138:    test:
    139:      filter:  sed -e "s/in 5 iterations/in 4 iterations/g"

    141: TEST*/