static char help[] = "Solves a linear system in parallel with KSP.\n\
Input parameters include:\n\
  -view_exact_sol   : write exact solution vector to stdout\n\
  -m <mesh_x>       : number of mesh points in x-direction\n\
  -n <mesh_y>       : number of mesh points in y-direction\n\n";

/*
  Include "petscksp.h" so that we can use KSP solvers.
*/

#include <petscksp.h>
#include <sys/time.h>
#include <float.h>

#include "mmloader.h"

int main(int argc, char **args)
{
   Vec x, b; /* approx solution, RHS */
   Vec u;
   Mat A;             /* linear system matrix */
   KSP ksp;           /* linear solver context */
   PetscReal residue; /* residue of solution error */
   PetscInt m, n, its;
   PetscReal norm_delta, norm_b; /* norm of solution error */

   char filename[PETSC_MAX_PATH_LEN], bfilename[PETSC_MAX_PATH_LEN];

   PetscFunctionBeginUser;

   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));

   PetscCall(PetscOptionsGetString(NULL, NULL, "-MTX", filename, PETSC_MAX_PATH_LEN, NULL));
   PetscCall(PetscOptionsGetString(NULL, NULL, "-RHS", bfilename, PETSC_MAX_PATH_LEN, NULL));

   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mat name: %s\n", filename));

   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "rhs name: %s\n", bfilename));

   // PetscCall(MatCreateFromMTX(&A, filename, PETSC_TRUE));
   PetscCall(load_binary_mat(&A, filename));
   PetscCall(MatGetSize(A, &m, &n));

   // PetscCall(VecCreateFromRHS(&b, bfilename));
   PetscCall(load_binary_rhs(&b, bfilename));

   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
   PetscCall(VecSetFromOptions(x));

   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                 Create the linear solver and set various options
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

   /*
      Set operators. Here the matrix that defines the linear system
      also serves as the preconditioning matrix.
   */
   PetscCall(KSPSetOperators(ksp, A, A));

   /*
      Set linear solver defaults for this problem (optional).
      - By extracting the KSP and PC contexts from the KSP context,
        we can then directly call any KSP and PC routines to set
        various options.
      - The following two statements are optional; all of these
        parameters could alternatively be specified at runtime via
        KSPSetFromOptions().  All of these defaults can be
        overridden at runtime, as indicated below.
   */
   PetscCall(VecNorm(b, NORM_2, &norm_b));
   // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Atol: %lf\n", 1.e-10 * norm_b));
   PetscCall(KSPSetTolerances(ksp, 1e-12, DBL_MIN, PETSC_DEFAULT, PETSC_DEFAULT));

   /*
     Set runtime options, e.g.,
         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
     These options will override those specified above as long as
     KSPSetFromOptions() is called _after_ any other customization
     routines.
   */
   PetscCall(KSPSetFromOptions(ksp));

   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                       Solve the linear system
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

   struct timeval t_start, t_stop;
   gettimeofday(&t_start, NULL);

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

   gettimeofday(&t_stop, NULL);
   double total_time = (t_stop.tv_sec - t_start.tv_sec) * 1000.0 + (t_stop.tv_usec - t_start.tv_usec) / 1000.0;

   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSPSolve Time: %lf ms\n", total_time));

   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                       Check the solution and clean up
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

   PetscCall(VecDuplicate(b, &u));
   PetscCall(MatMult(A, x, u));
   PetscCall(VecAXPY(b, -1.0, u));
   PetscCall(VecNorm(b, NORM_2, &norm_delta));

   PetscCall(KSPGetResidualNorm(ksp, &residue));
   PetscCall(KSPGetIterationNumber(ksp, &its));

   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "norm_delta %lf, norm_b %lf \n", norm_delta, norm_b));
   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (norm_delta) / (norm_b), its));
   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residue of error %g iterations %" PetscInt_FMT "\n", (double)residue, its));

   /*
      Free work space.  All PETSc objects should be destroyed when they
      are no longer needed.
   */
   // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Free work space\n"));

   PetscCall(KSPDestroy(&ksp));
   PetscCall(VecDestroy(&x));
   PetscCall(VecDestroy(&b));
   PetscCall(MatDestroy(&A));

   /*
      Always call PetscFinalize() before exiting a program.  This routine
        - finalizes the PETSc libraries as well as MPI
        - provides summary and diagnostic information if certain runtime
          options are chosen (e.g., -log_view).
   */
   PetscCall(PetscFinalize());

   return 0;
}

/*TEST

   build:
      requires: !single

   test:
      suffix: chebyest_1
      args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_monitor_short

   test:
      suffix: chebyest_2
      args: -m 80 -n 80 -ksp_pc_side right -pc_type ksp -ksp_ksp_type chebyshev -ksp_ksp_max_it 5 -ksp_ksp_chebyshev_esteig 0.9,0,0,1.1 -ksp_esteig_ksp_type cg -ksp_monitor_short

   test:
      args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always

   test:
      suffix: 2
      nsize: 2
      args: -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always

   test:
      suffix: 3
      args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

   test:
      suffix: 4
      args: -pc_type eisenstat -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

   test:
      suffix: 5
      nsize: 2
      args: -ksp_monitor_short -m 5 -n 5 -mat_view draw -ksp_gmres_cgs_refinement_type refine_always -nox
      output_file: output/ex2_2.out

   test:
      suffix: bjacobi
      nsize: 4
      args: -pc_type bjacobi -pc_bjacobi_blocks 1 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres

   test:
      suffix: bjacobi_2
      nsize: 4
      args: -pc_type bjacobi -pc_bjacobi_blocks 2 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres -ksp_view

   test:
      suffix: bjacobi_3
      nsize: 4
      args: -pc_type bjacobi -pc_bjacobi_blocks 4 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres

   test:
      suffix: qmrcgs
      args: -ksp_type qmrcgs -pc_type ilu
      output_file: output/ex2_fbcgs.out

   test:
      suffix: qmrcgs_2
      nsize: 3
      args: -ksp_type qmrcgs -pc_type bjacobi
      output_file: output/ex2_fbcgs_2.out

   test:
      suffix: fbcgs
      args: -ksp_type fbcgs -pc_type ilu

   test:
      suffix: fbcgs_2
      nsize: 3
      args: -ksp_type fbcgsr -pc_type bjacobi

   test:
      suffix: groppcg
      args: -ksp_monitor_short -ksp_type groppcg -m 9 -n 9

   test:
      suffix: mkl_pardiso_cholesky
      requires: mkl_pardiso
      args: -ksp_type preonly -pc_type cholesky -mat_type sbaij -pc_factor_mat_solver_type mkl_pardiso

   test:
      suffix: mkl_pardiso_lu
      requires: mkl_pardiso
      args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mkl_pardiso

   test:
      suffix: pipebcgs
      args: -ksp_monitor_short -ksp_type pipebcgs -m 9 -n 9

   test:
      suffix: pipecg
      args: -ksp_monitor_short -ksp_type pipecg -m 9 -n 9

   test:
      suffix: pipecgrr
      args: -ksp_monitor_short -ksp_type pipecgrr -m 9 -n 9

   test:
      suffix: pipecr
      args: -ksp_monitor_short -ksp_type pipecr -m 9 -n 9

   test:
      suffix: pipelcg
      args: -ksp_monitor_short -ksp_type pipelcg -m 9 -n 9 -pc_type none -ksp_pipelcg_pipel 2 -ksp_pipelcg_lmax 2
      filter: grep -v "sqrt breakdown in iteration"

   test:
      suffix: sell
      args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -m 9 -n 9 -mat_type sell

   test:
      requires: mumps
      suffix: sell_mumps
      args: -ksp_type preonly -m 9 -n 12 -mat_type sell -pc_type lu -pc_factor_mat_solver_type mumps -pc_factor_mat_ordering_type natural

   test:
      suffix: telescope
      nsize: 4
      args: -m 100 -n 100 -ksp_converged_reason -pc_type telescope -pc_telescope_reduction_factor 4 -telescope_pc_type bjacobi

   test:
      suffix: umfpack
      requires: suitesparse
      args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type umfpack

   test:
      suffix: spqr
      requires: suitesparse
      args: -ksp_type preonly -pc_type qr -pc_factor_mat_solver_type spqr

   test:
     suffix: pc_symmetric
     args: -m 10 -n 9 -ksp_converged_reason -ksp_type gmres -ksp_pc_side symmetric -pc_type cholesky

   test:
      suffix: pipeprcg
      args: -ksp_monitor_short -ksp_type pipeprcg -m 9 -n 9

   test:
      suffix: pipeprcg_rcw
      args: -ksp_monitor_short -ksp_type pipeprcg -recompute_w false -m 9 -n 9

   test:
      suffix: pipecg2
      requires: !defined(PETSC_HAVE_THREADSAFETY)
      args: -ksp_monitor_short -ksp_type pipecg2 -m 9 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}

   test:
      suffix: pipecg2_2
      requires: !defined(PETSC_HAVE_THREADSAFETY)
      nsize: 4
      args: -ksp_monitor_short -ksp_type pipecg2 -m 15 -n 9 -ksp_norm_type {{preconditioned unpreconditioned natural}}

   test:
      suffix: hpddm
      nsize: 4
      requires: hpddm !__float128 !__fp16
      filter: sed -e "s/ iterations 9/ iterations 8/g"
      args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{single double}shared output} -ksp_pc_side {{left right}shared output}

   test:
      suffix: hpddm___float128
      output_file: output/ex2_hpddm.out
      nsize: 4
      requires: hpddm __float128
      filter: sed -e "s/ iterations 9/ iterations 8/g"
      args: -ksp_converged_reason -ksp_type hpddm -ksp_hpddm_precision {{double quadruple}shared output} -ksp_pc_side {{left right}shared output}

   test:
      suffix: symmetric_pc
      nsize: 1
      args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc -ksp_pc_side symmetric

   test:
      suffix: symmetric_pc2
      nsize: 1
      args: -ksp_monitor -ksp_type gmres -pc_type bjacobi -sub_pc_type icc -ksp_pc_side symmetric -pc_bjacobi_blocks 2

 TEST*/