program main

c*********************************************************************72
c
cc helmholtz_openmp() solves a discretized Helmholtz equation.
c
c  Discussion:
c
c    The two dimensional region given is:
c
c      -1 <= X <= +1
c      -1 <= Y <= +1
c
c    The region is discretized by a set of M by N nodes:
c
c      P(I,J) = ( X(I), Y(J) )
c
c      X(I) = ( - M + 2*I - 1 ) / ( M - 1 )
c      Y(J) = ( - N + 2*J - 1 ) / ( N - 1 )
c
c    The Helmholtz equation for the scalar function U(X,Y) is
c
c      - Uxx(X,Y) -Uyy(X,Y) + ALPHA * U(X,Y) = F(X,Y)
c
c    where ALPHA is a positive constant.  We suppose that Dirichlet
c    boundary conditions are specified, that is, that the value of
c    U(X,Y) is given for all points along the boundary.
c
c    We suppose that the right hand side function F(X,Y) is specified in 
c    such a way that the exact solution is
c
c      U(X,Y) = ( 1 - X^2 ) * ( 1 - Y^2 )
c
c    Using standard finite difference techniques, the second derivatives
c    of U can be approximated by linear combinations of the values
c    of U at neighboring points.  Using this fact, the discretized
c    differential equation becomes a set of linear equations of the form:
c
c      A * U = F
c
c    These linear equations are then solved using a form of the Jacobi 
c    iterative method with a relaxation factor.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    Original Fortran77 version by Joseph Robicheaux, Sanjiv Shah.
c    This version by John Burkardt
c
      implicit none

      include 'omp_lib.h'

      double precision alpha
      integer it_max
      integer m
      integer n
      double precision omega
      double precision tol
      double precision wtime

      call timestamp ( )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'helmholtz_openmp():'
      write ( *, '(a)' ) '  Fortran77/OpenMP version'
      write ( *, '(a)' ) '  A program for the 2D Helmholtz equation.'
      write ( *, '(a)' ) '  Parallel execution under OpenMP.'
      write ( *, '(a,i8)' ) 
     &  '  The number of processors = ', omp_get_num_procs ( )
      write ( *, '(a,i8)' ) 
     &  '  The number of threads    = ', omp_get_max_threads ( )
c
c  Set the parameters.
c
      m = 64
      n = 64
      alpha = 0.5D+00
      omega = 0.9D+00
      tol = 1.0D-08
      it_max = 100

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The region is [-1,1] x [-1,1].'
      write ( *, '(a,i8)' ) 
     &  '  The number of X nodes is M = ', m
      write ( *, '(a,i8)' ) 
     &  '  The number of Y nodes is N = ', n
      write ( *, '(a,g14.6)' ) 
     &  '  The scalar coefficient is ALPHA = ', alpha
      write ( *, '(a,g14.6)' ) 
     &  '  The relaxation value is OMEGA = ', omega
      write ( *, '(a,g14.6)' ) 
     &  '  The error tolerance is TOL = ', tol
      write ( *, '(a,i8)' ) 
     &  '  The maximum number of Jacobi iterations is IT_MAX = ', it_max
c
c  Call the driver routine.
c
      wtime = omp_get_wtime ( )

      call driver ( m, n, it_max, alpha, omega, tol )

      wtime = omp_get_wtime ( ) - wtime

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  Elapsed wall clock time = ', wtime
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'helmholtz_openmp():'
      write ( *, '(a)' ) '  Normal end of execution.'

      call timestamp ( )

      stop
      end
      subroutine driver ( m, n, it_max, alpha, omega, tol )

c*********************************************************************72
c
cc driver() allocates arrays and solves the problem.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    Original Fortran77 version by Joseph Robicheaux, Sanjiv Shah.
c    This version by John Burkardt
c
c  Input:
c
c    integer M, N, the number of grid points in the X and Y directions.
c
c    integer IT_MAX, the maximum number of Jacobi iterations allowed.
c
c    double precision ALPHA, the scalar coefficient in the
c    Helmholtz equation.
c
c    double precision OMEGA, the relaxation parameter, which
c    should be strictly between 0 and 2.  For a pure Jacobi method,
c    use OMEGA = 1.
c
c    double precision TOL, an error tolerance for the linear
c    equation solver.
c
      implicit none

      integer m
      integer n

      double precision alpha
      double precision f(m,n)
      integer i
      integer it_max
      integer j
      double precision omega
      double precision tol
      double precision u(m,n)
c
c  Initialize the right hand side F.
c
      call initialize ( m, n, alpha, f )
c
c  Initialize the solution estimate U.
c
c$omp parallel 
c$omp& shared ( m, n, u )
c$omp& private ( i, j )

c$omp do
        do j = 1, n
          do i = 1, m
            u(i,j) = 0.0
          end do
        end do
c$omp end do

c$omp end parallel
c
c  Solve the Helmholtz equation.
c
      call jacobi ( m, n, alpha, omega, u, f, tol, it_max )
c
c  Determine the error.
c
      call error_check ( m, n, u )

      return
      end
      subroutine error_check ( m, n, u )

c*********************************************************************72
c
cc error_check() determines the error in the numerical solution.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    Original Fortran77 version by Joseph Robicheaux, Sanjiv Shah.
c    This version by John Burkardt
c
c  Input:
c
c    integer M, N, the number of grid points in the X and Y directions.
c
c    double precision U(M,N), the solution of the Helmholtz equation 
c    at the grid points.
c
      implicit none

      integer m
      integer n

      double precision dx
      double precision dy
      double precision error_norm
      integer i
      integer j
      double precision u(m,n)
      double precision u_exact
      double precision u_norm
      double precision u_true_norm
      double precision uxx
      double precision uyy
      double precision x
      double precision y

      dx = 2.0D+00 / dble ( m - 1 )
      dy = 2.0D+00 / dble ( n - 1 )

      u_norm = 0.0D+00
c
c  1: Compute u_norm
c
c$omp parallel 
c$omp& shared ( m, n, u ) 
c$omp& private ( i, j )

c$omp do reduction ( + : u_norm )
      do j = 1, n
        do i = 1, m
          u_norm = u_norm + u(i,j)**2
        end do
      end do
c$omp end do

c$omp end parallel

      u_norm = sqrt ( u_norm / m / n )
!
!  2: Compute U_True norm.
!
      u_true_norm = 0.0D+00

c$omp parallel 
c$omp& shared ( dx, dy, m, n, u ) 
c$omp& private ( i, j, u_exact, uxx, uyy, x, y ) 

c$omp do reduction ( + : u_true_norm )
      do j = 1, n
        do i = 1, m
          x = dble ( 2 * i - m - 1 ) / dble ( m - 1 )
          y = dble ( 2 * j - n - 1 ) / dble ( n - 1 )
          call exact ( x, y, u_exact, uxx, uyy )
          u_true_norm = u_true_norm + u_exact**2
        end do
      end do
c$omp end do

c$omp end parallel

      u_true_norm = sqrt ( u_true_norm / m / n )
c
c  3: Rescale U
c
      u = u * u_true_norm / u_norm
      u_norm = u_true_norm
!
!  4: Compute U-Utrue norm
!
      error_norm = 0.0D+00

c$omp parallel 
c$omp& shared ( dx, dy, m, n, u ) 
c$omp& private ( i, j, u_exact, uxx, uyy, x, y ) 

c$omp do reduction ( + : error_norm )
      do j = 1, n
        do i = 1, m
          x = -1.0D+00 + dx * dble ( i - 1 )
          y = -1.0D+00 + dy * dble ( j - 1 )
          call exact ( x, y, u_exact, uxx, uyy )
          error_norm = error_norm + ( u(i,j) - u_exact )**2
        end do
      end do
c$omp end do

c$omp end parallel

      error_norm = sqrt ( error_norm / m / n )

      write ( *, '(a)' ) ' '
      write ( *, '(a,g14.6)' ) '  RMS U :       ', u_norm
      write ( *, '(a,g14.6)' ) '  RMS UEXACT :  ', u_true_norm
      write ( *, '(a,g14.6)' ) '  RMS U-UEXACT: ', error_norm

      return
      end
      subroutine exact ( x, y, u, uxx, uyy )

c*********************************************************************72
c
cc exact() returns the exact value of the solution and derivatives.
c
c  Discussion:
c
c    This function can be used to compute the error, and to formulate the
c    appropriate boundary conditions and right hand side for any
c    desired solution.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision X, Y, the point at which the values are needed.
c
c  Output:
c
c    double precision U, UXX, UYY, the values of the function,
c    and its second derivatives with respect to X and Y.
c
      implicit none

      double precision u
      double precision uxx
      double precision uyy
      double precision x
      double precision y

      u = ( 1.0D+00 - x * x ) * ( 1.0D+00 - y * y )
      uxx = -2.0D+00 * ( 1.0D+00 + y ) * ( 1.0D+00 - y )
      uyy = -2.0D+00 * ( 1.0D+00 + x ) * ( 1.0D+00 - x )

      return
      end
      subroutine initialize ( m, n, alpha, f )

c*********************************************************************72
c
cc initialize() initializes the right hand side.
c
c  Discussion:
c
c    The routine assumes that the exact solution and its second
c    derivatives are given by the routine EXACT.
c
c    The appropriate Dirichlet boundary conditions are determined
c    by getting the value of U returned by EXACT.
c
c    The appropriate right hand side function is determined by
c    having EXACT return the values of U, UXX and UYY, and setting
c
c      F = -UXX - UYY + ALPHA * U
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    Original Fortran77 version by Joseph Robicheaux, Sanjiv Shah.
c    This version by John Burkardt
c
c  Input:
c
c    integer M, N, the number of grid points in the X and Y directions.
c
c    double precision ALPHA, the scalar coefficient in the
c    Helmholtz equation.  ALPHA should be positive.
c
c  Output:
c
c    double precision F(M,N): the right hand side function.
c
      implicit none

      integer m
      integer n

      double precision alpha
      double precision dx
      double precision dy
      double precision f(m,n)
      double precision f_norm
      integer i
      integer j
      double precision u_exact
      double precision uxx
      double precision uyy
      double precision x
      double precision y
c
c  Zero out F.
c
c$omp parallel
c$omp& shared ( f, m, n )
c$omp& private ( i, j )

c$omp do
      do j = 1, n
        do i = 1, m
          f(i,j) = 0.0D+00
        end do
      end do
c$omp end do

c$omp end parallel
c
c  Set boundary conditions.
c
c$omp parallel
c$omp& shared ( f, m, n )
c$omp& private ( i, j, x, y )

c$omp do
      do i = 1, m
        j = 1
        y = dble ( 2 * j - n - 1 ) / dble ( n - 1 )
        x = dble ( 2 * i - m - 1 ) / dble ( m - 1 )
        call exact ( x, y, u_exact, uxx, uyy )
        f(i,j) = u_exact
      end do
c$omp end do

c$omp do
      do i = 1, m
        j = n
        y = dble ( 2 * j - n - 1 ) / dble ( n - 1 )
        x = dble ( 2 * i - m - 1 ) / dble ( m - 1 )
        call exact ( x, y, u_exact, uxx, uyy )
        f(i,j) = u_exact
      end do
c$omp end do

c$omp do
      do j = 1, n
        i = 1
        x = dble ( 2 * i - m - 1) / dble ( m - 1 )
        y = dble ( 2 * j - n - 1) / dble ( n - 1 )
        call exact ( x, y, u_exact, uxx, uyy )
        f(i,j) = u_exact
      end do
c$omp end do

c$omp do
      do j = 1, n
        i = m
        x = dble ( 2 * i - m - 1 ) / dble ( m - 1 )
        y = dble ( 2 * j - n - 1 ) / dble ( n - 1 )
        call exact ( x, y, u_exact, uxx, uyy )
        f(i,j) = u_exact
      end do
c$omp end do

c$omp end parallel
c
c  Set the right hand side F.
c
      f_norm = 0.0D+00

c$omp parallel
c$omp& shared ( alpha, dx, dy, f, m, n ) 
c$omp& private ( i, j, u_exact, uxx, uyy, x, y )

c$omp do reduction ( + : f_norm )
      do j = 2, n - 1
        do i = 2, m - 1
          x = dble ( 2 * i - m - 1 ) / dble ( m - 1 )
          y = dble ( 2 * j - n - 1 ) / dble ( n - 1 )
          call exact ( x, y, u_exact, uxx, uyy )
          f(i,j) = uxx + uyy - alpha * u_exact
          f_norm = f_norm + f(i,j)**2
        end do
      end do
c$omp end do

c$omp end parallel

      f_norm = sqrt ( f_norm / m / n )

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

      return
      end
      subroutine jacobi ( m, n, alpha, omega, u, f, tol, it_max )

c*********************************************************************72
c
cc jacobi() applies the Jacobi iterative method to solve the linear system.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    02 November 2024
c
c  Author:
c
c    Original Fortran77 version by Joseph Robicheaux, Sanjiv Shah.
c    This version by John Burkardt
c
c  Input:
c
c    integer M, N, the number of grid points in the X and Y directions.
c
c    double precision ALPHA, the scalar coefficient in the
c    Helmholtz equation.  ALPHA should be positive.
c
c    double precision OMEGA, the relaxation parameter, which
c    should be strictly between 0 and 2.  For a pure Jacobi method,
c    use OMEGA = 1.
c
c    double precision U(M,N), the solution estimate.
c
c    double precision F(M,N), values of the right hand side function 
c    for the Helmholtz equation at the grid points.
c
c    double precision TOL, an error tolerance for the linear
c    equation solver.
c
c    integer IT_MAX, the maximum number of Jacobi iterations allowed.
c
c  Output:
c
c    double precision U(M,N), the updated solution estimate.
c
      implicit none

      integer m
      integer n

      double precision alpha
      double precision ax
      double precision ay
      double precision b
      double precision dx
      double precision dy
      double precision error_norm
      double precision f(m,n)
      integer i
      integer it
      integer it_max
      integer j
      double precision omega
      double precision resid
      double precision tol
      double precision u(m,n)
      double precision uold(m,n)
c
c  Initialize the coefficients.
c
      dx = 2.0D+00 / dble ( m - 1 )
      dy = 2.0D+00 / dble ( n - 1 )

      ax = 1.0D+00 / dx**2
      ay = 1.0D+00 / dy**2
      b  = - 2.0D+00 / dx**2 - 2.0D+00 / dy**2 - alpha

      do it = 1, it_max
c
c  Copy new solution into old.
c
c$omp parallel 
c$omp& shared ( m, n, u, uold )
c$omp& private ( i, j )

c$omp do
        do j = 1, n
          do i = 1, m
            uold(i,j) = u(i,j)
          end do
        end do
c$omp end do

c$omp end parallel
c
c  Compute stencil, residual, and update.
c
      error_norm = 0.0

c$omp parallel
c$omp& shared ( ax, ay, b, f, m, n, omega, u, uold )
c$omp& private ( i, j, resid )

c$omp do reduction ( + : error_norm )
        do j = 1, n
          do i = 1, m
c
c  Evaluate the residual.
c
            if ( i .eq. 1 .or. 
     &           i .eq. m .or. 
     &           j .eq. 1 .or. 
     &           j .eq. n ) then

              resid = uold(i,j) - f(i,j)

            else

              resid = ( ax * ( uold(i-1,j) + uold(i+1,j) ) 
     &                + ay * ( uold(i,j-1) + uold(i,j+1) ) 
     &                + b  *   uold(i,j) - f(i,j) ) / b

            end if
c
c  Update the solution.
c
            u(i,j) = uold(i,j) - omega * resid
c
c  Accumulate the residual error.
c
            error_norm = error_norm + resid * resid

          end do
        end do
c$omp end do

c$omp end parallel
c
c  Error check.
c
        error_norm = sqrt ( error_norm / m / n )

        write ( *, '(a,g14.6)' ) 
     &  '  RMS residual            ', error_norm

        if ( error_norm .le. tol .and. 1 < it ) then
          write ( *, *) '  Test = ', tol
          exit
        end if

      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a,i8)' ) 'Number of iterations ', it

      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