Actual source code: ex1.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_type POSITIVE_DEFINITE option (or comment in the PCFactorSetShiftType() line below):
     13:  *    the method will now successfully converge.
     14:  *
     15:  * Contributed by Victor Eijkhout 2003.
     16:  */

     18: #include <petscksp.h>

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

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

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

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

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

     64:   PetscCall(MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES));
     65:   PetscCall(MatSetValues(A, 1, &j, 1, &i, &v, INSERT_VALUES));
     66:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     67:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
     68:   PetscCall(VecAssemblyBegin(B));
     69:   PetscCall(VecAssemblyEnd(B));

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

     78:   PetscCall(KSPSetType(solver, KSPCG));
     79:   PetscCall(KSPSetInitialGuessNonzero(solver, PETSC_TRUE));

     81:   /*
     82:    * ILU preconditioner;
     83:    * this will break down unless you add the Shift line,
     84:    * or use the -pc_factor_shift_positive_definite option */
     85:   PetscCall(KSPGetPC(solver, &prec));
     86:   PetscCall(PCSetType(prec, PCILU));
     87:   /* PetscCall(PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE)); */

     89:   PetscCall(KSPSetFromOptions(solver));
     90:   PetscCall(KSPSetUp(solver));

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

    101:   /*
    102:    * Solve the system;
    103:    * without the shift this will diverge with
    104:    * an indefinite preconditioner
    105:    */
    106:   PetscCall(KSPSolve(solver, B, X));
    107:   PetscCall(KSPGetConvergedReason(solver, &reason));
    108:   if (reason == KSP_DIVERGED_INDEFINITE_PC) {
    109:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nDivergence because of indefinite preconditioner;\n"));
    110:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n"));
    111:   } else if (reason < 0) {
    112:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nOther kind of divergence: this should not happen.\n"));
    113:   } else {
    114:     PetscCall(KSPGetIterationNumber(solver, &its));
    115:   }

    117:   PetscCall(VecDestroy(&X));
    118:   PetscCall(VecDestroy(&B));
    119:   PetscCall(VecDestroy(&D));
    120:   PetscCall(MatDestroy(&A));
    121:   PetscCall(KSPDestroy(&solver));
    122:   PetscCall(PetscFinalize());
    123:   return 0;
    124: }

    126: /*TEST

    128:    test:
    129:       args: -pc_factor_shift_type positive_definite

    131: TEST*/