program main

!*****************************************************************************80
!
!! MAIN is the main program for POISSON_SERIAL.
!
!  Discussion:
!
!    POISSON_SERIAL is a program for solving the Poisson problem.
!
!    This program runs serially.  Its output is used as a benchmark for
!    comparison with similar programs run in a parallel environment.
!
!    The Poisson equation
!
!      DEL^2 U(X,Y) = F(X,Y)
!
!    is solved on the unit square [0,1] x [0,1] using a grid of [0,NX+1] by
!    [0,NY+1] evenly spaced points.  
!
!    The boundary conditions and F are set so that the exact solution is
!
!      U(X,Y) = sin ( X * Y )
!
!    The Jacobi iteration is repeatedly applied until convergence is detected.
!
!  Modified:
!
!    18 September 2002
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: nx = 9
  integer, parameter :: ny = 9

  real a(0:nx+1,0:ny+1)
  real anorm
  real b(0:nx+1,0:ny+1)
  real bnorm
  logical converged
  real diff
  real dx
  real dy
  real error
  real f(0:nx+1,0:ny+1)
  integer i
  integer it
  integer, parameter :: it_max = 100
  integer j
  real, parameter :: tolerance = 0.001E+00
  real u_exact
  real x
  real y
!
!  Print a message.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POISSON_SERIAL:'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  A NON-MPI program for solving the Poisson equation.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  The number of interior X grid lines is ', nx
  write ( *, '(a,i8)' ) '  The number of interior Y grid lines is ', ny
!
!  Initialize the data.
!
  call init_serial ( nx, ny, dx, dy, f, a, b )

  write ( *, '(a,g14.6)' ) '  The X grid spacing is ', dx
  write ( *, '(a,g14.6)' ) '  The Y grid spacing is ', dy
!
!  Do the iteration.
!
  converged = .false.

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    'Step    ||U||         ||Unew||     ||Unew-U||     ||Unew-Exact||'
  write ( *, '(a)' ) ' '

  do it = 1, it_max
!
!  Perform one Jacobi sweep, computing B from A.
!
    call sweep_serial ( nx, ny, dx, dy, f, a, b )
!
!  Perform a second Jacobi sweep, computing A from B.
!
    call sweep_serial ( nx, ny, dx, dy, f, b, a )
!
!  Check for convergence.
!
    anorm = 0.0E+00
    bnorm = 0.0E+00
    diff = 0.0E+00
    do j = 1, ny
      do i = 1, nx
        anorm = max ( anorm, abs ( a(i,j) ) )
        bnorm = max ( bnorm, abs ( b(i,j) ) )
        diff = max ( diff, abs ( a(i,j) - b(i,j) ) )
      end do
    end do

    error = 0.0E+00
    do j = 1, ny
      y = real ( j ) / real ( ny + 1 )
      do i = 1, nx
        x = real ( i ) / real ( nx + 1 )
        u_exact = sin ( x * y )
        error = max ( error, abs ( u_exact - a(i,j) ) )
      end do
    end do

    write ( *, '(i3,3g14.6,f14.8)' ) it, bnorm, anorm, diff, error

    if ( diff <= tolerance ) then
      converged = .true.
      exit
    end if

  end do

  if ( converged ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POISSON_SERIAL:'
    write ( *, '(a)' ) '  The iteration has converged'
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POISSON_SERIAL:'
    write ( *, '(a)' ) '  The iteration has NOT converged'
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POISSON_SERIAL:'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine init_serial ( nx, ny, dx, dy, f, a, b )

!*****************************************************************************80
!
!! INIT_SERIAL initializes the arrays.
!
!  Modified:
!
!    19 September 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, the X and Y grid dimensions.
!
!    Output, real DX, DY, the spacing between grid points.
!
!    Output, real F(0:NX+1,0:NY+1), the right hand side data.
!
!    Output, real A(0:NX+1,0:NY+1), B(0:NX+1,0:NY+1), the initial
!    solution estimates.
!
  implicit none

  integer nx
  integer ny

  real a(0:nx+1,0:ny+1)
  real b(0:nx+1,0:ny+1)
  real dx
  real dy
  real f(0:nx+1,0:ny+1)
  real fnorm
  integer i
  integer j
  real u
  real unorm
  real x
  real y

  dx = 1.0E+00 / real ( nx + 1 )
  dy = 1.0E+00 / real ( ny + 1 )
!
!  Set the initial guesses to zero.
!
  a(0:nx+1,0:ny+1) = 0.0E+00
  b(0:nx+1,0:ny+1) = 0.0E+00
!
!  Just for debugging, compute the max norm of the exact solution 
!  on the interior nodes.
!
  unorm = 0.0E+00
  do j = 1, ny
    y = real ( j ) / real ( ny + 1 )
    do i = 1, nx
      x = real ( i ) / real ( nx + 1 )
      u = sin ( x * y )
      unorm = max ( unorm, abs ( u ) )
    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INIT_SERIAL:'
  write ( *, '(a,g14.6)' ) &
    '  Max norm of exact solution U at interior nodes = ', unorm
!
!  The "boundary" entries of F will store the boundary values of the solution.
!
  fnorm = 0.0E+00

  x = 0.0E+00
  do j = 0, ny+1
    y = real ( j ) / real ( ny + 1 )
    f(0,j) = sin ( x * y )
    a(0,j) = f(0,j)
    b(0,j) = f(0,j)
    fnorm = max ( fnorm, f(0,j) )
  end do

  x = 1.0E+00
  do j = 0, ny+1
    y = real ( j ) / real ( ny + 1 )
    f(nx+1,j) = sin ( x * y )
    a(nx+1,j) = f(nx+1,j)
    b(nx+1,j) = f(nx+1,j)
    fnorm = max ( fnorm, f(nx+1,j) )
  end do

  y = 0.0E+00
  do i = 0, nx+1
    x = real ( i ) / real ( nx + 1 )
    f(i,0) = sin ( x * y )
    a(i,0) = f(i,0)
    b(i,0) = f(i,0)
    fnorm = max ( fnorm, f(i,0) )
  end do

  y = 1.0E+00
  do i = 0, nx+1
    x = real ( i ) / real ( nx + 1 )
    f(i,ny+1) = sin ( x * y )
    a(i,ny+1) = f(i,ny+1)
    b(i,ny+1) = f(i,ny+1)
    fnorm = max ( fnorm, f(i,ny+1) )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INIT_SERIAL:'
  write ( *, '(a,g14.6)' ) '  Max norm of boundary values = ', fnorm
!
!  The "interior" entries of F will store the right hand sides 
!  of the Poisson equation.
!
  fnorm = 0.0E+00
  do j = 1, ny
    y = real ( j ) / real ( ny + 1 )
    do i = 1, nx
      x = real ( i ) / real ( nx + 1 )
      f(i,j) = - ( x**2 + y**2 ) * sin ( x * y )
      fnorm = max ( fnorm, abs ( f(i,j) ) )
    end do
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INIT_SERIAL:'
  write ( *, '(a,g14.6)' ) &
    '  Max norm of right hand side F at interior nodes = ', fnorm

  return
end
subroutine sweep_serial ( nx, ny, dx, dy, f, u, unew )

!*****************************************************************************80
!
!! SWEEP_SERIAL carries out one step of the Jacobi iteration.
!
!  Modified:
!
!    17 September 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NX, NY, the X and Y grid dimensions.
!
!    Input, real DX, DY, the spacing between grid points.
!
!    Input, real F(0:NX+1,0:NY+1), the right hand side data.
!
!    Input, real U(0:NX+1,0:NY+1), the previous solution estimate.
!
!    Output, real UNEW(0:NX+1,0:NY+1), the updated solution estimate.
!
  implicit none

  integer nx
  integer ny

  real dx
  real dy
  real f(0:nx+1,0:ny+1)
  integer i
  integer j
  real u(0:nx+1,0:ny+1)
  real unew(0:nx+1,0:ny+1)

  do j = 1, ny
    do i = 1, nx

      unew(i,j) = ( u(i-1,j) + u(i,j+1) + u(i,j-1) + u(i+1,j) ) / 4.0E+00 &
        - f(i,j) * dx * dy 

    end do
  end do

  return
end