Actual source code: ex70.c


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

  4: #include <petscksp.h>

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

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

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

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

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

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

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

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

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

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

 85: /*TEST

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

 92: TEST*/