! CSC/PATC Introduction to parallel programming using MPI and OpenMP
! (c) CSC 2013
program jacobi
  implicit none

  integer, parameter :: rk = kind(1d0)
  integer, parameter :: N = 10000

  real(kind=rk), dimension(0:N+1) :: rho  ! Source term
  real(kind=rk), dimension(0:N+1) :: phi, phi_old ! Potential

  integer :: i
  real(kind=rk) :: error = 1000.0
  real(kind=rk) :: residual
  real(kind=rk), parameter :: tol = 6.80e-7

  real(kind=rk) :: a, b  ! interval
  real(kind=rk) :: h, h2, x
  integer :: niter=1

  a = 0.0;
  b = 1.0;

  h = (b-a) / N;
  h2 = h*h;

  x = a;
  ! Initialize rho and phis
  do i = 0, N+1
     x = x + h;
     rho(i) = 2*x*x - 1
  end do

  phi(:) = 0.0
  phi_old(:) = 0.0

  ! Start your modifications from here
  ! and insert appropriate openMP pragmas
  !$omp parallel
  do while (error > tol)
     ! insert pragma which makes sure that error is not assigned
     ! to zero before the statement in while

     ! How many threads should execute the following block?
     phi_old(:) = phi(:)
     error = 0


     !$omp do private(residual) reduction(+:error)
     do i=1,N
        phi(i) = 0.5 * (phi_old(i+1) + phi_old(i-1) + h2 * rho(i))
        residual = phi_old(i+1) - 2*phi_old(i) + phi_old(i-1) + h2 * rho(i)
        error = error + residual*residual
     end do
     !$omp end do

     ! How many threads should execute the following block?
     error = sqrt(error)
     niter = niter + 1

     ! Print output every 100th iteration
     ! How many threads should execute the following block?
     if (mod(niter, 100) == 0) then
        write(*,'(A,I0,A,G0)') "Error at iteration ", niter, ": ", error
     end if

     if (niter > 5000) then
        write(*,'(A,I0,A)') "Failed to converge in ", niter, " iterations"
     end if
  end do
  !$OMP end parallel

  write(*,'(A,I0,A)') "Jacobi solver converged in ", niter, " iterations"
  write(*,'(A,G0)') "Error in solution: ", error
end program jacobi
