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;

 30:   PetscInitialize(&argc, &argv, 0, help);
 31:   comm = MPI_COMM_SELF;

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

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

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

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

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

 78:   KSPSetType(ksp, KSPCG);
 79:   KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);

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

 94:   KSPSetFromOptions(ksp);
 95:   KSPSetUp(ksp);

 97:   /*
 98:    * Now that the factorisation is done, show the pivots;
 99:    * note that the last one is negative. This in itself is not an error,
100:    * but it will make the iterative method diverge.
101:    */
102:   PCFactorGetMatrix(pc, &M);
103:   VecDuplicate(B, &D);
104:   MatGetDiagonal(M, D);
105:   PetscPrintf(PETSC_COMM_WORLD, "\nPivots:\n\n");
106:   VecView(D, 0);

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

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

135: /*TEST

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

140: TEST*/