program gs_test
  implicit none
  integer, parameter :: n = 10
  integer :: i, j, k, maxit
  double precision, dimension(0:n,0:n) :: u
  double precision x, y, pi, h, error
  !
  ! Test gs_sor (gauss-seidel by successive-over-relaxation)
  ! Rozwiązywanie cząstkowego równania różniczkowego:
  !   Laplace(u(x,u)) = (4-pi**2)*exp(2*x)*cos(pi*y), 0 < x,y < 1
  !   u(x,0) = exp(2*x)                   dla 0 < x < 1
  !   u(1,y) = exp(2.)*cos(pi*y)          dla 0 < y < 1
  !   u(x,1) = -exp(2*x)                  dla 0 < x < 1
  !   u(0,y) = cos(pi*y)                  dla 0 < y < 1
  ! Rozwiązanie dokładne: u(x,y) = exp(2*x)*cos(pi*y)
  !
  h = 1./n
  error = 1.e-3
  maxit = 1000

  ! boundary conditions and trial solution
  pi = 4*atan(1d0)
  u = 0d0
  do i=0,n
     x = i*h
     y = x
     u(i,0) = exp(2d0*x)
     u(n,i) = exp(2d0)*cos(pi*y)
     u(i,n) = -u(i,0)
     u(0,i) = cos(pi*y)
  end do

  call gs(u,n,h,error,maxit)

  write(*, "(11f8.3)") ((u(i,j), j=0,n), i=0,n)

end Program gs_test

subroutine gs(w,n,h,error,maxit)
  implicit none
  double precision, dimension(0:n,0:n) :: w, resid
  double precision :: Laplace
  double precision :: rhs, h, hh, x, y, u, error, res
  integer :: n, nm, i, j, k, maxit
  nm = n - 1
  hh = h*h
  do k = 1, maxit
     do i = 1, nm
        do j = 1, nm
           x = i*h
           y = j*h
           w(i,j) = -.25 * (hh * rhs(x,y) - &
                (w(i-1,j)+w(i+1,j)+w(i,j-1)+w(i,j+1)))
        end do
     end do
     ! residual
     res = 0
     do i = 1,nm
        do j = 1,nm
           x = i*h
           y = j*h
           resid(i,j) = rhs(x,y) - Laplace(w,n,h,i,j)
        end do
     end do
     if ( sum(abs(resid)) .le. error ) then
        write(*,*) "# Number of iterations: ", k
        return
     end if
  end do
  write(*,*) "# In gs: Number of iterations exceeded.", maxit
  return
end subroutine gs

  
function rhs(x,y)
  double precision :: x, y, pi, rhs
  ! right hand side of the equation D^2 = f(x,y)
  data pi/3.141592/
  rhs = (4-pi*pi)*exp(2*x)*cos(pi*y)
end function rhs


function Laplace(w,n,h,i,j)
  double precision :: Laplace, h
  integer :: n, i, j
  double precision, dimension(0:n,0:n) :: w
  ! nomerical Laplacian(w(i,j))  value at i, j grid point
  ! h - step sizes in x and y directions
  Laplace = (w(i-1,j)+w(i+1,j)+w(i,j-1)+w(i,j+1)-4*w(i,j))/(h*h)
end function Laplace


