Actual source code: ex70.c

  1: static char help[] = "Solves an ill-conditioned tridiagonal linear system with KSP for testing GMRES breakdown tolerance.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b, u; /* approx solution, RHS, exact solution */
  8:   Mat         A;       /* linear system matrix */
  9:   KSP         ksp;     /* linear solver context */
 10:   PetscInt    i, n = 10, col[3];
 11:   PetscMPIInt size;
 12:   PetscScalar value[3];

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 16:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 17:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 19:   /*
 20:      Create vectors.  Note that we form 1 vector from scratch and
 21:      then duplicate as needed.
 22:   */
 23:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 24:   PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
 25:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
 26:   PetscCall(VecSetFromOptions(x));
 27:   PetscCall(VecDuplicate(x, &b));
 28:   PetscCall(VecDuplicate(x, &u));

 30:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 31:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 32:   PetscCall(MatSetFromOptions(A));
 33:   PetscCall(MatSetUp(A));

 35:   /*
 36:      Set big off-diag values to make the system ill-conditioned
 37:   */
 38:   value[0] = 10.0;
 39:   value[1] = 2.0;
 40:   value[2] = 1.0;
 41:   for (i = 1; i < n - 1; i++) {
 42:     col[0] = i - 1;
 43:     col[1] = i;
 44:     col[2] = i + 1;
 45:     PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
 46:   }
 47:   i      = n - 1;
 48:   col[0] = n - 2;
 49:   col[1] = n - 1;
 50:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 51:   i        = 0;
 52:   col[0]   = 0;
 53:   col[1]   = 1;
 54:   value[0] = 2.0;
 55:   value[1] = -1.0;
 56:   PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
 57:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 58:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 60:   PetscCall(VecSet(u, 1.0));
 61:   PetscCall(MatMult(A, u, b));

 63:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 64:   PetscCall(KSPSetOperators(ksp, A, A));
 65:   PetscCall(KSPSetFromOptions(ksp));
 66:   PetscCall(KSPSolve(ksp, b, x));

 68:   PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
 69:   PetscCall(PetscOptionsInsertString(NULL, "-ksp_type preonly -ksp_initial_guess_nonzero false"));
 70:   PetscCall(PetscOptionsClearValue(NULL, "-ksp_converged_reason"));
 71:   PetscCall(KSPSetFromOptions(ksp));
 72:   PetscCall(KSPSolve(ksp, b, x));

 74:   PetscCall(VecDestroy(&x));
 75:   PetscCall(VecDestroy(&u));
 76:   PetscCall(VecDestroy(&b));
 77:   PetscCall(MatDestroy(&A));
 78:   PetscCall(KSPDestroy(&ksp));

 80:   PetscCall(PetscFinalize());
 81:   return 0;
 82: }

 84: /*TEST

 86:    test:
 87:       requires: double !complex
 88:       args: -ksp_rtol 1e-18 -pc_type sor -ksp_converged_reason -ksp_gmres_breakdown_tolerance 1.e-10
 89:       output_file: output/ex70.out

 91: TEST*/