Program g_s
  implicit none
  integer, parameter :: n = 11
  integer :: i, j, k, maxit
  real*8, dimension(-n:n,-n:n) :: w
  real*8 h, error

  ! --------------------------------------------------------------------
  ! In this example we solve (iteratively) the equation
  ! 1 + Laplacian(w(x,y)) = 0 in the region -1<=x<=1, -1<=y<=1.
  ! The boundary conditions are the following:
  ! w(-1,y)=w(+1,y)=0, w(x,-1)=w(x,1)=0.
  !
  ! n      - 2n+1 = number of points in the interval [-1,1]
  ! h      - grid distance
  ! w(i,j) - function to be determined
  ! maxit  - maximal number of allowed iterations
  ! error  - averege difference between two successive iterations
  ! --------------------------------------------------------------------

  h = 1./n
  error = 1.e-4
  maxit = 1000

  ! boundary conditions and trial solution
  w = 0.
  
  call gs(w,n,h,error,maxit)
  print "(2f15.6)", (i*h, w(i,0), i=-n,n)
  
contains

  subroutine gs(w,n,h,error,maxit)
    implicit none
    real*8, dimension(-n:n,-n:n) :: w, resid
    real*8 :: h, hh, error
    integer :: n, nm, i, j, k, maxit
    ! iterations for 1 + D^2 w(x,y) = 0, -1<=x<=1, -1<=y<=1,
    ! w(0,j) = w(n,j) = 0, w(i,0) = w(i,n) = 0,
    nm = n - 1
    hh = h*h
    do k = 1, maxit
       do i = -nm, nm
          do j = -nm, nm
             ! mixed (Gauss-Seidel iteration)
             w(i,j) = .25 * (hh + &
                  w(i-1,j) + w(i,j-1) + &      ! new part of w
                  w(i+1,j) + w(i,j+1))         ! old part of w
          end do
       end do
       do i = -nm,nm
          do j = -nm,nm
             resid(i,j) = 1 + (-4.*w(i,j) + &
                  w(i-1,j) + w(i,j-1) + w(i+1,j) + w(i,j+1))/hh
          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

end Program g_s
