program main

c*********************************************************************72
c
cc MAIN is the main program for POISSON_SERIAL.
c
c  Discussion:
c
c    POISSON_SERIAL is a program for solving the Poisson problem.
c
c    This program runs serially.  Its output is used as a benchmark for
c    comparison with similar programs run in a parallel environment.
c
c    The Poisson equation
c
c      DEL^2 U(X,Y) = F(X,Y)
c
c    is solved on the unit square [0,1] x [0,1] using a grid of [0,NX+1] by
c    [0,NY+1] evenly spaced points.  
c
c    The boundary conditions and F are set so that the exact solution is
c
c      U(X,Y) = sin ( X * Y )
c
c    The Jacobi iteration is repeatedly applied until convergence is detected.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    18 September 2002
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Gropp, Ewing Lusk, Anthony Skjellum,
c    Using MPI: Portable Parallel Programming with the
c    Message-Passing Interface,
c    Second Edition,
c    MIT Press, 1999,
c    ISBN: 0262571323,
c    LC: QA76.642.G76.
c
c    Marc Snir, Steve Otto, Steven Huss-Lederman, David Walker, 
c    Jack Dongarra,
c    MPI: The Complete Reference,
c    Volume I: The MPI Core,
c    Second Edition,
c    MIT Press, 1998,
c    ISBN: 0-262-69216-3,
c     LC: QA76.642.M65.
c
      implicit none

      integer nx
      integer ny

      parameter ( nx = 9 )
      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 it_max
      integer j
      real tolerance
      real u_exact
      real x
      real y

      it_max = 100
      tolerance = 0.001E+00
c
c  Print a message.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'POISSON_SERIAL:'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  A program for solving the Poisson equation.'
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 
     &  '  The number of interior X grid lines is ', nx
      write ( *, '(a,i6)' ) 
     &  '  The number of interior Y grid lines is ', ny
c
c  Initialize the data.
c
      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
c
c  Do the iteration.
c
      converged = .false.

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

      do it = 1, it_max
c
c  Perform one Jacobi sweep, computing B from A.
c
        call sweep_serial ( nx, ny, dx, dy, f, a, b )
c
c  Perform a second Jacobi sweep, computing A from B.
c
        call sweep_serial ( nx, ny, dx, dy, f, b, a )
c
c  Check for convergence.
c
        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 .le. 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 )

c*********************************************************************72
c
cc INIT_SERIAL initializes the arrays.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    19 September 2002
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer NX, NY, the X and Y grid dimensions.
c
c    Output, real DX, DY, the spacing between grid points.
c
c    Output, real F(0:NX+1,0:NY+1), the right hand side data.
c
c    Output, real A(0:NX+1,0:NY+1), B(0:NX+1,0:NY+1), the initial
c    solution estimates.
c
      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 )
c
c  Set the initial guesses to zero.
c
      do i = 0, nx+1
        do j = 0, ny+1
          a(i,j) = 0.0E+00
          b(i,j) = 0.0E+00
        end do
      end do
c
c  Just for debugging, compute the max norm of the exact solution 
c  on the interior nodes.
c
      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
c
c  The "boundary" entries of F will store the boundary values of the solution.
c
      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
c
c  The "interior" entries of F will store the right hand sides 
c  of the Poisson equation.
c
      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 )

c*********************************************************************72
c
cc SWEEP_SERIAL carries out one step of the Jacobi iteration.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    17 September 2002
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer NX, NY, the X and Y grid dimensions.
c
c    Input, real DX, DY, the spacing between grid points.
c
c    Input, real F(0:NX+1,0:NY+1), the right hand side data.
c
c    Input, real U(0:NX+1,0:NY+1), the previous solution estimate.
c
c    Output, real UNEW(0:NX+1,0:NY+1), the updated solution estimate.
c
      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