program main

c*********************************************************************72
c
cc poisson_2d_test() tests poisson_2d().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    03 October 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer nx
      integer ny
      external rhs1
      double precision t_start
      double precision t_stop
      double precision, external :: u_exact1
      double precision, external :: uxxyy_exact1
      double precision wtime

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'poisson_2d_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test poisson_2d().'
c
c  Problem 1.
c
      t_start = wtime ( )
      nx = 161
      ny = 161
      call poisson_2d ( nx, ny, rhs1, u_exact1 )
      t_stop = wtime ( )
      write ( *, '(a)' ) ''
      write ( *, '(a,g14.6)' ) '  Wall clock time in seconds is ', 
     &  t_stop - t_start
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'poisson_2d_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop 0
      end
      subroutine rhs1 ( nx, ny, f )

c*********************************************************************72
c
cc rhs1() initializes the right hand side "vector" for case #1.
c
c  Discussion:
c
c    It is convenient for us to set up RHS as a 2D array.  However, each
c    entry of RHS is really the right hand side of a linear system of the
c    form
c
c      A * U = F
c
c    In cases where U(I,J) is a boundary value, then the equation is simply
c
c      U(I,J) = F(i,j)
c
c    and F(I,J) holds the boundary data.
c
c    Otherwise, the equation has the form
c
c      (1/DX^2) * ( U(I+1,J)+U(I-1,J)+U(I,J-1)+U(I,J+1)-4*U(I,J) ) = F(I,J)
c
c    where DX is the spacing and F(I,J) is the value at X(I), Y(J) of
c
c      pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 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  Output:
c
c    double precision F(NX,NY), the right hand side data.
c
      implicit none

      integer nx
      integer ny

      double precision f(nx,ny)
      double precision fnorm
      integer i
      integer j
      double precision rms_norm
      double precision u_exact1
      double precision uxxyy_exact1
      double precision x
      double precision y
c
c  The "boundary" entries of F will store the boundary values of the solution.
c
      x = 0.0D+00
      do j = 1, ny
        y = dble ( j - 1 ) / dble ( ny - 1 )
        f(1,j) = u_exact1 ( x, y )
      end do

      x = 1.0D+00
      do j = 1, ny
        y = dble ( j - 1 ) / dble ( ny - 1 )
        f(nx,j) = u_exact1 ( x, y )
      end do

      y = 0.0D+00
      do i = 1, nx
        x = dble ( i - 1 ) / dble ( nx - 1 )
        f(i,1) = u_exact1 ( x, y )
      end do

      y = 1.0D+00
      do i = 1, nx
        x = dble ( i - 1 ) / dble ( nx - 1 )
        f(i,ny) = u_exact1 ( x, y )
      end do
c
c  The "interior" entries of F store the right hand sides 
c  of the Poisson equation.
c
      do j = 2, ny - 1
        y = dble ( j - 1 ) / dble ( ny - 1 )
        do i = 2, nx - 1
          x = dble ( i - 1 ) / dble ( nx - 1 )
          f(i,j) = - uxxyy_exact1 ( x, y )
        end do
      end do

      fnorm = rms_norm ( nx, ny, f ) 

      write ( *, '(a,g14.6)' ) '  RMS of F = ', fnorm

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end
      function u_exact1 ( x, y )

c*********************************************************************72
c
cc u_exact1() evaluates the exact solution for case #1.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 October 2011
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision X, Y, the coordinates of a point.
c
c  Output:
c
c    double precision U_EXACT1, the value of the exact solution at (X,Y).
c
      implicit none

      double precision, parameter :: pi = 3.141592653589793D+00
      double precision u_exact1
      double precision x
      double precision y

      u_exact1 = sin ( pi * x * y )

      return
      end
      function uxxyy_exact1 ( x, y )

c*********************************************************************72
c
cc uxxyy_exact1() evaluates ( d/dx d/dx + d/dy d/dy ) for case #1.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 October 2011
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision X, Y, the coordinates of a point.
c
c  Output:
c
c    double precision UXXYY_EXACT1, the value of 
c    ( d/dx d/dx + d/dy d/dy ) of the exact solution at (X,Y).
c
      implicit none

      double precision, parameter :: pi = 3.141592653589793D+00
      double precision uxxyy_exact1
      double precision x
      double precision y

      uxxyy_exact1 = - pi * pi * ( x * x + y * y ) * sin ( pi * x * y )

      return
      end
      function wtime ( )

c*********************************************************************72
c
cc wtime() returns a reading of the wall clock time.
c
c  Discussion:
c
c    To get the elapsed wall clock time, call WTIME before and after a given
c    operation, and subtract the first reading from the second.
c
c    This function is meant to suggest the similar routines:
c
c      "omp_get_wtime ( )" in OpenMP,
c      "MPI_Wtime ( )" in MPI,
c      and "tic" and "toc" in MATLAB.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    27 April 2009
c
c  Author:
c
c    John Burkardt
c
c  Output:
c
c    dble WTIME, the wall clock reading, in seconds.
c
      implicit none

      integer clock_max
      integer clock_rate
      integer clock_reading
      double precision wtime

      call system_clock ( clock_reading, clock_rate, clock_max )

      wtime = dble ( clock_reading ) / dble ( clock_rate )

      return
      end