1: static char help[] = "Test interface functions of KSPIDR with a nonsymmetric matrix.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b, xe; /* approx solution, RHS, exact solution */
  8:   Mat         A;        /* linear system matrix */
  9:   KSP         ksp;      /* linear solver context */
 10:   PC          pc;       /* preconditioner */
 11:   PetscInt    i, Istart, Iend, n = 40, s = 0;
 12:   PetscReal   nrm;
 13:   PetscRandom rctx;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 19:   /* Create bidiagonal matrix */
 20:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 21:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 22:   PetscCall(MatSetFromOptions(A));
 23:   PetscCall(MatSetUp(A));
 24:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 26:   for (i = Istart; i < Iend; i++) {
 27:     PetscCall(MatSetValue(A, i, i, 2.0, INSERT_VALUES));
 28:     if (i < n - 1) PetscCall(MatSetValue(A, i, i + 1, 1.0 / (i + 1), INSERT_VALUES));
 29:   }

 31:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 32:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 34:   PetscCall(MatCreateVecs(A, &x, &b));
 35:   PetscCall(VecDuplicate(x, &xe));
 36:   PetscCall(VecSet(xe, 1.0));
 37:   PetscCall(MatMult(A, xe, b));

 39:   /* Create linear solver context */
 40:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 41:   PetscCall(KSPSetOperators(ksp, A, A));
 42:   PetscCall(KSPSetType(ksp, KSPIDR));
 43:   PetscCall(KSPGetPC(ksp, &pc));
 44:   PetscCall(PCSetType(pc, PCNONE));

 46:   /* Pass specific IDR options, including user-provided PetscRandom */
 47:   PetscCall(KSPIDRSetS(ksp, 6));
 48:   PetscCall(KSPIDRSetCosine(ksp, 0.8));
 49:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
 50:   PetscCall(PetscRandomSetFromOptions(rctx));
 51:   PetscCall(KSPIDRSetRandom(ksp, rctx));
 52:   PetscCall(PetscRandomDestroy(&rctx));
 53:   PetscCall(KSPSetFromOptions(ksp));

 55:   PetscCall(KSPSolve(ksp, b, x));

 57:   PetscCall(VecAXPY(x, -1.0, xe));
 58:   PetscCall(VecNorm(x, NORM_2, &nrm));
 59:   PetscCall(KSPIDRGetS(ksp, &s));
 60:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Absolute error (s=%" PetscInt_FMT ") = %1.6f\n", s, (double)nrm));

 62:   /* Solve again with different s */
 63:   PetscCall(KSPIDRSetS(ksp, 4));
 64:   PetscCall(KSPSolve(ksp, b, x));

 66:   PetscCall(VecAXPY(x, -1.0, xe));
 67:   PetscCall(VecNorm(x, NORM_2, &nrm));
 68:   PetscCall(KSPIDRGetS(ksp, &s));
 69:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Absolute error (s=%" PetscInt_FMT ") = %1.6f\n", s, (double)nrm));

 71:   /* Free work space */
 72:   PetscCall(KSPDestroy(&ksp));
 73:   PetscCall(VecDestroy(&x));
 74:   PetscCall(VecDestroy(&xe));
 75:   PetscCall(VecDestroy(&b));
 76:   PetscCall(MatDestroy(&A));

 78:   PetscCall(PetscFinalize());
 79:   return 0;
 80: }

 82: /*TEST

 84:    test:
 85:       suffix: 1

 87:    test:
 88:       args: -pc_type jacobi -ksp_pc_side right
 89:       suffix: 2

 91: TEST*/