! Ax=b is solved using PETSc but now, instead of each proc (ideally) writing its own
!       submatrix, proc 0 writes a (PETSc) binary file containing both A and b which is
!       then read and solved in parallel. 
! Usage: mpirun -np 2 ./a.out -ksp_type cg -pc_type bjacobi -ksp_rtol 1.0e-15 -ksp_max_it 50000


#include "include/finclude/petsc.h"
#include "include/finclude/petscmat.h"
#include "include/finclude/petscvec.h"
#include "include/finclude/petscpc.h"
#include "include/finclude/petscksp.h"
#include "include/finclude/petscviewer.h"

!      implicit real*8(a-h,o-z)
      REAL(8), POINTER :: hs(:), cmatrx(:), rld(:)
      INTEGER, POINTER :: nlrn(:),lrn(:)

      integer :: rank, ierr
      integer :: n, nnz
      real*8 fm1, fnorm, time1, time2
      Mat     :: A
      Vec     :: x, b, xtrue
      KSP     :: ksp
      integer :: its
      PetscViewer :: view,uview, bview

      call PetscInitialize(Petsc_Null_Character, ierr)
      call MPI_Comm_Rank(Petsc_Comm_World, rank, ierr)


! First proc reads/writes a given matrix (in std CSR format) alongwith the rhs to a PETSc binary file
      if (rank == 0) then
         open(2,file='input/out1.cmatrx',form='unformatted',status='old')
         read(2) nnp
         allocate(hs(nnp))
         allocate(rld(nnp))
         allocate(nlrn(nnp+1))
c
c read matrix, right hand side, and initial guess
         print *,'reading nlrn ...'
         read(2) (nlrn(j),j=1,nnp+1)
         print *,'reading lrn ...'
         allocate(lrn(nlrn(nnp+1)))
         allocate(cmatrx(nlrn(nnp+1)))
         read(2)((lrn(i),i=nlrn(j)+1,nlrn(j+1)),j=1,nnp)
         print *,'reading cmatrix and rld ...'
         read(2) ((cmatrx(i),i=nlrn(j)+1,nlrn(j+1)),rld(j),hs(j),j=1,nnp)
         close(2)


c
c out0.solve stores the result from CG with JACOBI preconditioner 
         open(12,file='input/out0.xsolve',form='unformatted',status='old')
         read(12)(hs(i),i=1,nnp)
         close(12)

         call MatCreateSeqAIJWithArrays(Petsc_Comm_Self, nnp, nnp, nlrn, lrn-1,cmatrx , A, ierr)
         call VecCreateSeqWithArray(Petsc_Comm_Self, nnp, rld, b, ierr)
         call VecCreateSeqWithArray(Petsc_Comm_Self, nnp, hs, xtrue, ierr)


! write A and rhs to a file:
         call PetscViewerBinaryOpen(Petsc_Comm_Self, 'matrhs.dat', File_Mode_Write, view, ierr)
         call PetscViewerSetFormat(view, Petsc_Viewer_Binary_Native, ierr)
         call MatView(A, view, ierr)
         call VecView(b, view, ierr)
         call PetscViewerDestroy(view, ierr)
         call MatDestroy(A, ierr)
         call VecDestroy(b, ierr)

! write xtrue to a binary file so we can compute the norm

         call PetscViewerBinaryOpen(Petsc_Comm_Self, 'xtrue.dat', File_Mode_Write, view, ierr)
         call PetscViewerSetFormat(view, Petsc_Viewer_Binary_Native, ierr)
         call VecView(xtrue, view, ierr)
         call PetscViewerDestroy(view, ierr)
         call VecDestroy(xtrue, ierr)

! deallocate array.

         deallocate(cmatrx, nlrn, lrn, rld, hs)

      endif
! Make sure that the procs do not try opening the file before it is created
      call MPI_Barrier(Petsc_Comm_World, ierr)

c
c     compute r = Ax-b, L2 norm, L infinity norm

! Now the matrix is read from the file and solved in parallel
      call PetscViewerBinaryOpen(Petsc_Comm_World, 'matrhs.dat', File_Mode_Read, view, ierr)

      call MatLoad(view, MATMPIAIJ, A, ierr)
      call VecLoad(view, VECMPI, b, ierr)
      call VecDuplicate(b, x, ierr)

! Noe read the xtrue:

      call PetscViewerBinaryOpen(Petsc_Comm_World, 'xtrue.dat', File_Mode_Read, view, ierr)
      call VecLoad(view, VECMPI, xtrue, ierr)

!---
      time1 = MPI_WTIME ()
      call KSPCreate(Petsc_Comm_World, ksp, ierr)
      call KSPSetOperators(ksp, A, A, Different_Nonzero_Pattern, ierr)
      call KSPSetFromOptions(ksp, ierr)
      call KSPSolve(ksp, b, x, ierr)
      time2 = MPI_WTIME ()

      call KSPGetIterationNumber(ksp, its, ierr)
      if (rank .eq. 0) then
         print*, "Number of iterations and Time in PETSc solver =", its, time2 - time1, " secs"
      end if


!         Print the alternate solution error information.
      fm1 = -1.0
      call VecAXPY (xtrue, fm1, x, ierr)
      call VecNorm (xtrue, NORM_2, fnorm, ierr)
    
      if (rank .eq. 0) then
        print*, '2 norm of the error from the provided & computed solution =',
     &    fnorm
      end if

      call VecNorm (xtrue, NORM_INFINITY, fnorm, ierr)
      if (rank .eq. 0) then
        print*, 'Maximum error from the provided & computed solution',
     &    ' (infinity norm) =', fnorm
      end if

      call VecNorm (xtrue, NORM_1, fnorm, ierr)
      if (rank .eq. 0) then
        print*, '1 norm of the error from the provided & computed solution',
     &   fnorm
      end if



      call KSPDestroy(ksp, ierr)
      call MatDestroy(A, ierr)
      call VecDestroy(b, ierr)
      call VecDestroy(x, ierr)
      call VecDestroy(xtrue, ierr)
      call PetscViewerDestroy(view, ierr)

      call PetscFinalize(ierr)
c
      stop
      end
c