subroutine poisson_2d ( nx, ny, rhs, u_exact )

c*********************************************************************72
c
cc poisson_2d() solves the Poisson equation in a 2D region (a square).
c
c  Discussion:
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 ( pi * x * y )
c
c    so that
c
c      - DEL^2 U(x,y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
c
c    The Jacobi iteration is repeatedly applied until convergence is detected.
c
c    For convenience in writing the discretized equations, we assume that NX = NY.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    04 October 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer NX, NY, the X and Y grid dimensions.
c
c    external RHS: the routine that evaluates the right hand side.
c
c    exteral U_EXACT: the function that evaluates the exact solution.
c
      implicit none

      integer nx
      integer ny

      logical converged
      double precision diff
      double precision dx
      double precision dy
      double precision error
      double precision f(nx,ny)
      integer i
      integer itnew
      integer itold
      integer j
      double precision rms_norm
      double precision tolerance
      double precision u(nx,ny)
      double precision u_exact
      double precision u_norm
      double precision udiff(nx,ny)
      double precision uexact(nx,ny)
      double precision unew(nx,ny)
      double precision unew_norm
      double precision x
      double precision y

      tolerance = 0.000001D+00
      dx = 1.0D+00 / dble ( nx - 1 )
      dy = dx
c
c  Print a message.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'poisson_2d():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  A program for solving the Poisson equation.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  -DEL^2 U = F(X,Y)'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  on the rectangle 0 <= X <= 1, 0 <= Y <= 1.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '  F(X,Y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )'
      write ( *, '(a)' ) ' '
      write ( *, '(a,i6)' ) 
     &  '  The number of interior X grid points is ', nx
      write ( *, '(a,i6)' ) 
     &  '  The number of interior Y grid points is ', ny
      write ( *, '(a,g14.6)' ) '  The X grid spacing is ', dx
      write ( *, '(a,g14.6)' ) '  The Y grid spacing is ', dy
c
c  Initialize the data.
c
      call rhs ( nx, ny, f )
c
c  Set the initial solution estimate.
c  We are "allowed" to pick up the boundary conditions exactly.
c
      do j = 1, ny
        do i = 1, nx
          if ( i == 1 .or. i == nx .or.
     &         j == 1 .or. j == ny ) then
            unew(i,j) = f(i,j)
          else
            unew(i,j) = 0.0D+00
          end if
        end do
      end do

      unew_norm = rms_norm ( nx, ny, unew )
c
c  Set up the exact solution.
c
      do j = 1, ny 
        y = dble ( j - 1 ) / dble ( ny - 1 )
        do i = 1, nx
          x = dble ( i - 1 ) / dble ( nx - 1 )
          uexact(i,j) = u_exact ( x, y )
        end do
      end do
      u_norm = rms_norm ( nx, ny, uexact )
      write ( *, '(a,g14.6)' ) '  L2 norm of exact solution = ', u_norm
c
c  Do the iteration.
c
      converged = .false.

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

      do j = 1, ny
        do i = 1, nx
          udiff(i,j) = unew(i,j) - uexact(i,j)
        end do
      end do

      error = rms_norm ( nx, ny, udiff )
      write ( *, '(2x,i4,2x,g14.6,2x,14x,2x,g14.6)' ) 
     &  0, unew_norm, error

      itnew = 0

      do

        itold = itnew
        itnew = itold + 500
c
c  Carry out 500 Jacobi steps in parallel before we come
c  back to check for convergence.
c
        call jacobi ( nx, ny, dx, dy, f, itold, itnew, u, unew )
c
c  We declare convergence if successive iterates change very little.
c
        u_norm = unew_norm
        unew_norm = rms_norm ( nx, ny, unew )
        udiff = unew(1:nx,1:ny) - u(1:nx,1:ny)
        diff = rms_norm ( nx, ny, udiff )
        udiff = unew(1:nx,1:ny) - uexact(1:nx,1:ny)
        error = rms_norm ( nx, ny, udiff )

        write ( *, '(2x,i6,2x,g14.6,2x,g14.6,2x,g14.6)' ) 
     %    itnew, unew_norm, diff, error

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

      end do

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

      return
      end
      function rms_norm ( m, n, a )

c*********************************************************************72
c
cc rms_norm() returns the RMS of a vector stored as an R8MAT.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 May 2008
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    Input, integer M, N, the number of rows and columns in A.
c
c    Input, double precision A(M,N), the data whose RMS is desired.
c
c    Output, double precision rms_norm, the RMS of A.
c
      implicit none

      integer m
      integer n

      double precision a(m,n)
      integer i
      integer j
      double precision rms_norm
      double precision value

      value = 0.0D+00
      do j = 1, n
        do i = 1, m
          value = value + a(i,j) * a(i,j)
        end do
      end do
      value = sqrt ( value / dble ( m * n ) )

      rms_norm = value

      return
      end
      subroutine jacobi ( nx, ny, dx, dy, f, itold, itnew, u, unew )

c*********************************************************************72
c
cc jacobi() carries out one step of the Jacobi iteration.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    26 October 2011
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer NX, NY, the X and Y grid dimensions.
c
c    double precision DX, DY, the spacing between grid points.
c
c    double precision F(NX,NY), the right hand side data.
c
c    integer ITOLD, the iteration index on input.
c
c    integer ITNEW, the desired iteration index on output.
c
c    double precision U(NX,NY), the previous solution estimate.
c
c  Output:
c
c    double precision UNEW(NX,NY), the updated solution estimate.
c
      implicit none

      integer nx
      integer ny

      double precision dx
      double precision dy
      double precision f(nx,ny)
      integer i
      integer it
      integer itnew
      integer itold
      integer j
      double precision u(nx,ny)
      double precision unew(nx,ny)

      do it = itold + 1, itnew
c
c  Save the current estimate.
c
        do j = 1, ny
          do i = 1, nx
            u(i,j) = unew(i,j)
          end do
        end do
c
c  Compute a new estimate.
c
        do j = 1, ny
          do i = 1, nx

            if ( i == 1 .or. i == nx .or. j == 1 .or. j == ny ) then
              unew(i,j) = f(i,j)
            else
              unew(i,j) = 0.25D+00 * ( 
     &          u(i-1,j) + u(i,j+1) + u(i,j-1) + u(i+1,j) 
     &            + f(i,j) * dx * dy )
            end if

          end do
        end do

      end do

      return
      end