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;

 32:   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:   MatCreateSeqAIJ(comm, 4, 4, 4, 0, &A);
 40:   VecCreateSeq(comm, 4, &B);
 41:   VecDuplicate(B, &X);
 42:   for (i = 0; i < 4; i++) {
 43:     v = 3;
 44:     MatSetValues(A, 1, &i, 1, &i, &v, INSERT_VALUES);
 45:     v = 1;
 46:     VecSetValues(B, 1, &i, &v, INSERT_VALUES);
 47:     VecSetValues(X, 1, &i, &v, INSERT_VALUES);
 48:   }

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

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

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

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

 78:   KSPSetType(solver, KSPCG);
 79:   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:   KSPGetPC(solver, &prec);
 86:   PCSetType(prec, PCILU);
 87:   /* PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE); */

 89:   KSPSetFromOptions(solver);
 90:   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:   PCFactorGetMatrix(prec, &M);
 98:   VecDuplicate(B, &D);
 99:   MatGetDiagonal(M, D);

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

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

126: /*TEST

128:    test:
129:       args: -pc_factor_shift_type positive_definite

131: TEST*/