program main
  implicit none

! Includes files for PetsC
#include "include/finclude/petsc.h"
#include "include/finclude/petscvec.h"
#include "include/finclude/petscvec.h90"
#include "include/finclude/petscmat.h"
#include "include/finclude/petscmat.h90"
#include "include/finclude/petscksp.h"
#include "include/finclude/petscviewer.h"
#include "include/finclude/petscpc.h"
  

! Declarations

  integer                    :: mpi_err
  integer                    :: comm
  integer                    :: rank, size
  integer                    :: nx, ny
  integer                    :: vector_size
  integer                    :: m_l, n_l
  integer                    :: i
  integer                    :: its
  double precision           :: a0, b0, lx, ly
  double precision           :: dx, dy, dx2, dy2
  PetscErrorCode             :: petsc_err
  Vec                        :: x,b
  Mat                        :: A
  KSP                        :: solver
  PC                         :: precond
  KSPConvergedReason         :: reason


! Initializations of MPI/PetsC Context
  comm = MPI_COMM_WORLD
  PETSC_COMM_WORLD = comm
  call MPI_Init(mpi_err)
  call PetscInitialize(PETSC_NULL_CHARACTER, petsc_err)

! Getting informations about MPI
  call MPI_Comm_rank(comm, rank, mpi_err)
  call MPI_Comm_size(comm, size, mpi_err)

! Read grid file and broadcast informations for all processus
  if(rank == 0) then
!     call read_grid('grid.dat'//achar(0), a0, b0, lx, ly, nx, ny)
     nx = 499
     ny = 499
     lx = 1
     ly = 10
  end if

  call MPI_Bcast(a0, 1, MPI_DOUBLE_PRECISION, 0, comm, mpi_err)
  call MPI_Bcast(b0, 1, MPI_DOUBLE_PRECISION, 0, comm, mpi_err)
  call MPI_Bcast(lx, 1, MPI_DOUBLE_PRECISION, 0, comm, mpi_err)
  call MPI_Bcast(ly, 1, MPI_DOUBLE_PRECISION, 0, comm, mpi_err)
  call MPI_Bcast(nx, 1, MPI_INTEGER, 0, comm, mpi_err)
  call MPI_Bcast(ny, 1, MPI_INTEGER, 0, comm, mpi_err)

! Initializations
  dx = lx / dble(nx)
  dy = ly / dble(ny)

  dx2 = dx**2
  dy2 = dy**2

  vector_size = (nx+1)*(ny+1)


! PetsC Vectors Initializations
  call VecCreate(comm, x, petsc_err)
  call VecSetSizes(x, PETSC_DECIDE, vector_size, petsc_err)
  call VecSetFromOptions(x, petsc_err)
  call VecGetLocalSize(x, m_l, petsc_err)

  call VecDuplicate(x, b, petsc_err)
  call VecZeroEntries (x, petsc_err)
  call VecZeroEntries (b, petsc_err)


! PetsC Matrix Initialization
  call MatCreate(comm, A, petsc_err)
  call MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, vector_size, vector_size, petsc_err)
  call MatSetFromOptions(A, petsc_err)

! Assembling Matrix and Vectors
  if(rank == 0) then
     call MatSetValue(A, 0, 0, 1.d0, INSERT_VALUES, petsc_err)
     call MatSetValue(A, vector_size-1, vector_size-1, 1.d0, INSERT_VALUES, petsc_err)
     call VecSetValue(b, 0, 0.d0, INSERT_VALUES, petsc_err)
     call VecSetValue(b, vector_size-1, 1.d0, INSERT_VALUES, petsc_err)

     do i=1, vector_size-2
        call MatSetValue(A, i, i, 2.d0*dx2, INSERT_VALUES, petsc_err)
        call MatSetValue(A, i, i-1, -1.d0*dx2, INSERT_VALUES, petsc_err)
        call MatSetValue(A, i, i+1, -1.d0*dx2, INSERT_VALUES, petsc_err)
        call VecSetValue(b, i, 0.d0, INSERT_VALUES, petsc_err)
     end do

  end if
  
  call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, petsc_err)
  call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, petsc_err)

  call VecAssemblyBegin(b, petsc_err)
  call VecAssemblyEnd(b, petsc_err) 
  call VecAssemblyBegin(x, petsc_err)
  call VecAssemblyEnd(x, petsc_err) 

  
!   call MatView(A, PETSC_VIEWER_STDOUT_WORLD, petsc_err)  
!   call VecView(b, PETSC_VIEWER_STDOUT_WORLD, petsc_err)


! PetsC Solver and Preconditionner Initialization
  call KSPCreate(comm, solver, petsc_err)
  call KSPSetOperators(solver, A, A, SAME_NONZERO_PATTERN, petsc_err)
  call KSPGetPC(solver, precond, petsc_err)
  call KSPSetFromOptions(solver, petsc_err)

! Solving
  call KSPSolve(solver, b, x, petsc_err)

  call KSPView(solver, PETSC_VIEWER_STDOUT_WORLD, petsc_err)
  
  call KSPGetIterationNumber(solver, its, petsc_err)
  call KSPGetConvergedReason(solver, reason, petsc_err)

  if(rank==0) then
     print*,'Number of iterations : ', its
     print*,'Convergence          : ', reason
  end if

!  call VecView(x, PETSC_VIEWER_STDOUT_WORLD, petsc_err)



! Finalizing and Destroy...
  call KSPDestroy(solver, petsc_err)
  call VecDestroy(x, petsc_err)
  call VecDestroy(b, petsc_err)
  call MatDestroy(A, petsc_err)
  call PetscFinalize(petsc_err)
  call MPI_Finalize(mpi_err)

end program main