program petsc_poisson_test

#include <finclude/petscdef.h>
  use petsc
  use petscvec
  use petscmat
  use petscksp
  
  implicit none
  
  PetscInt :: nsize
  PetscErrorCode :: ierr
  Vec :: RHS, Sol
  Mat :: A_Mat
  KSP :: Solver
  double precision :: dx, X, F_X, min, max, PI
  integer :: ii, Row, Col
    
  nsize = 32
  min = 0.0d0
  max = 1.0d0
  
  call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
  
  !-- Create PETSc Vector for RHS and Sol
  call VecCreate(PETSC_COMM_WORLD, RHS, ierr)
  call VecSetSizes(RHS, PETSC_DECIDE, nsize, ierr)
  call VecSetFromOptions(RHS, ierr)
  call VecDuplicate(RHS, Sol, ierr)
  
  !-- Set the RHS value
  PI = acos(- 1.0d0)
  dx =(max - min)/(nsize - 1)
  do ii = 1, nsize
    X =(ii - 1)* dx
    F_X = - X *(X + 3)* exp(X)
    print*, X, F_X
    call VecSetValue(RHS, ii - 1, F_X, INSERT_VALUES, ierr)
  end do
  
  call VecAssemblyBegin(RHS, ierr)
  call VecAssemblyBegin(Sol, ierr)
  call VecAssemblyEnd(RHS, ierr)
  call VecAssemblyEnd(Sol, ierr)
  
  !-- View the RHS and SolVector
  call VecView(RHS, PETSC_VIEWER_STDOUT_WORLD, ierr)
  
  !-- Create PETSx Vector for the operator
  call MatCreate(PETSC_COMM_WORLD, A_Mat, ierr)
  call MatSetSizes &
        (A_Mat, PETSC_DECIDE, PETSC_DECIDE, nsize, nsize, ierr)
  call MatSetFromOptions(A_Mat, ierr)
  
  !-- Fill in the matrix value for 1-D A_Mat operator
  do ii = 1, nsize
    Row = ii - 1
    Col = Row 
    if(ii == 1)then
      call MatSetValue(A_Mat, Row, Col    , + 2.0d0, INSERT_VALUES, ierr)
      call MatSetValue(A_Mat, Row, Col + 1, - 1.0d0, INSERT_VALUES, ierr)
    else if(ii == nsize)then
      call MatSetValue(A_Mat, Row, Col    , + 2.0d0, INSERT_VALUES, ierr)
      call MatSetValue(A_Mat, Row, Col - 1, - 1.0d0, INSERT_VALUES, ierr)
    else
      call MatSetValue(A_Mat, Row, Col    , + 2.0d0, INSERT_VALUES, ierr)
      call MatSetValue(A_Mat, Row, Col + 1, - 1.0d0, INSERT_VALUES, ierr)
      call MatSetValue(A_Mat, Row, Col - 1, - 1.0d0, INSERT_VALUES, ierr)
    end if
  end do
  
  call MatAssemblyBegin(A_Mat, MAT_FINAL_ASSEMBLY, ierr)
  call MatAssemblyEnd(A_Mat, MAT_FINAL_ASSEMBLY, ierr)
  call MatView(A_Mat, PETSC_VIEWER_STDOUT_WORLD, ierr)
  
  !-- Create KSP Linear Solver Context
  call KSPCreate(PETSC_COMM_WORLD, Solver, ierr)
  call KSPSetOperators &
        (Solver, A_Mat, A_Mat, DIFFERENT_NONZERO_PATTERN, ierr)
  call KSPSetTolerances &
        (Solver, 1.0d-7, PETSC_DEFAULT_DOUBLE_PRECISION,     &
         PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
        
  call KSPSetFromOptions(Solver, ierr)
  
  !-- Solve the Linear System
  call KSPSolve(Solver, RHS, Sol, ierr)
  
  !-- View the solver Info and solution
  call KSPView(Solver, PETSC_VIEWER_STDOUT_WORLD, ierr)
  call VecView(Sol, PETSC_VIEWER_STDOUT_WORLD, ierr)
  
  call PetscFinalize(ierr)
  
end program petsc_poisson_test