static char help[] = "Reads in a PETSc binary matrix A, vec b for solving Ax=b.\n\
  -fmat <input_file> : matrix to load.\n\
  -fvec <input_file> : vec to load.\n\n";

#include "petsc.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  PetscErrorCode ierr;

  KSP            solver;
  PC             prec;

  Mat            A;
  Vec            b;
  Vec            x;

  char           matfile[PETSC_MAX_PATH_LEN], vecfile[PETSC_MAX_PATH_LEN]; 
  PetscViewer    matfd, vecfd;


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

#if defined(PETSC_USE_COMPLEX)
  SETERRQ(1,"This example does not work with complex numbers");
#endif
  int size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);


  ierr = PetscOptionsGetString(PETSC_NULL,"-fmat",matfile,PETSC_MAX_PATH_LEN-1,PETSC_NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetString(PETSC_NULL,"-fvec",vecfile,PETSC_MAX_PATH_LEN-1,PETSC_NULL);CHKERRQ(ierr);

  MatCreate(MPI_COMM_SELF, &A);
  VecCreate(MPI_COMM_SELF, &b);


  /* Read matrix and RHS */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,matfile,FILE_MODE_READ,&matfd);CHKERRQ(ierr);
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,vecfile,FILE_MODE_READ,&vecfd);CHKERRQ(ierr);
  ierr = MatLoad(A, matfd);CHKERRQ(ierr);
  ierr = VecLoad(b, vecfd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(matfd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(vecfd);CHKERRQ(ierr);

  ierr = VecDuplicate(b,&x); CHKERRQ(ierr);

  /* create ksp solver */

  ierr = KSPCreate(MPI_COMM_SELF, &solver); CHKERRQ(ierr);
  ierr = KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);

  ierr = KSPSetType(solver,KSPBCGSL); CHKERRQ(ierr);
  //ierr = KSPSetInitialGuessNonzero(solver,PETSC_TRUE); CHKERRQ(ierr);
  KSPSetTolerances(solver, 1e-12, 1e-15, PETSC_DEFAULT, 1000);
  ierr = KSPGetPC(solver,&prec); CHKERRQ(ierr);
  ierr = PCSetType(prec,PCILU); CHKERRQ(ierr);
  ierr = PCFactorSetColumnPivot(prec, 1.0); CHKERRQ(ierr);
  ierr = PCFactorSetShiftType(prec, MAT_SHIFT_NONZERO); CHKERRQ(ierr);

  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
  ierr = KSPSetUp(solver); CHKERRQ(ierr);

  ierr = KSPSolve(solver,b,x); CHKERRQ(ierr);

  ierr = PetscFinalize();CHKERRQ(ierr);
  return 0;
}