#include "petscksp.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Vec            x, b, u;      /* approx solution, RHS, exact solution */
  Mat            A, C, CT, B;            /* linear system matrix */
  KSP            ksp;         /* linear solver context */
  PC             pc;           /* preconditioner context */
  PetscReal      norm, bnorm,
      tol=1.e-14; /* norm of solution error */
  PetscErrorCode ierr;
  PetscInt       i,nopt = 10,col[3],its;
  PetscMPIInt    size;
  PetscScalar    neg_one = -1.0,one = 1.0,value[3];
  PetscViewer    viewer1, viewer2;
  
  PetscInitialize(&argc,&args,(char *)0,NULL);
  MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  PetscOptionsGetInt(PETSC_NULL,"-n",&nopt,PETSC_NULL);
  
  PetscViewerBinaryOpen(PETSC_COMM_WORLD,
                        "Ablock.bin",FILE_MODE_READ,
                        &viewer1);
  /* load the matrix A */
  MatCreate(PETSC_COMM_WORLD,&A);
  MatSetType(A,MATSEQAIJ);
  MatLoad(A, viewer1);
  /* ----- */
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"b.bin",
                               FILE_MODE_READ,&viewer2);
  VecCreate(PETSC_COMM_WORLD,&b);
  VecLoad(b,viewer2);
  /* ---- */
  MatView(A,PETSC_VIEWER_STDOUT_WORLD);
  VecView(b,PETSC_VIEWER_STDOUT_WORLD);
  /* Create a index sets */
  IS is_row1; /* 0 and 1 */
  IS is_row2; /* 2 and 3 */
  
  PetscInt bs = 2, n1=1, n2=1, inputindices[]={0,1}; /* block size is 2 */
    
  ISCreateBlock(PETSC_COMM_SELF,bs,n1,inputindices,PETSC_COPY_VALUES,&is_row1);
  ISCreateBlock(PETSC_COMM_SELF,bs,n2,inputindices,PETSC_COPY_VALUES,&is_row2);

  ISView(is_row1,PETSC_VIEWER_STDOUT_SELF);
  ISView(is_row2,PETSC_VIEWER_STDOUT_SELF);

  KSPCreate(PETSC_COMM_WORLD,&ksp);
  KSPGetPC(ksp,&pc);
  KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);

  PCFieldSplitSetIS(pc, "1", is_row1);
  PCFieldSplitSetIS(pc, "2", is_row2);
  
  KSPSetFromOptions(ksp);
  
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
    Solve linear system
  */
  VecDuplicate(b,&x);
  ierr = KSPSolve(ksp,b,x);

  ierr = VecView(x,
                 PETSC_VIEWER_STDOUT_WORLD );
  
   KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);

  VecNorm(b,NORM_2,&bnorm);
  KSPGetResidualNorm(ksp,&norm);
  KSPGetIterationNumber(ksp,&its);
  if (norm > tol) {
    PetscPrintf(PETSC_COMM_WORLD,
                "Norm of error %G,Iterations %D\n",norm/bnorm,its);
  }
  VecDestroy(&x);
  VecDestroy(&b);
  MatDestroy(&A);
  KSPDestroy(&ksp);
  
  MatDestroy(&A);
  PetscFinalize();
  return 0;
}