!===================================================================
! Author: Saleh Elmohamed,  saleh@npac.syr.edu
! Function: an HPF to solve a Laplace equation over a rectangle by the 
!         Red/Black SOR method.
! 1996, NPAC.
! Hardware platform: the IBM SP2.
!===================================================================

      program Laplace_SOR

      interface
         subroutine SOR_iteration(Grid,oldGrid,int_x,int_y,error,omega)
            double precision, dimension(1:int_x,1:int_y)::Grid,oldGrid
            integer int_x,int_y
            double precision error,omega
         end subroutine

         subroutine compute_val_of_U(Uval,x,y)
            double precision Uval,x,y
         end subroutine

         subroutine assign_boundary_values(U,int_x,int_y)
            double precision, dimension(1:int_x,1:int_y) :: U
            integer int_x,int_y
         end subroutine
      end interface

!!! Set the array grid sizes ....
 
      parameter(int_x = 100)
      parameter(int_y = 100)
      parameter(thresh = 1.0d-7)   !! or the level of tolerance ...

      double precision error,omega

      double precision, dimension(1:int_x,1:int_y) :: Grid, oldGrid
!HPF$ DISTRIBUTE Grid (*,CYCLIC)
!HPF$ DISTRIBUTE oldGrid (*,CYCLIC)

      integer iteration_num

!!! Set the overrelaxation parameter, omega

      omega = 1.9995d0

      write(*,*) '>>>>   The Mesh Resolution is .... '
      write(*,*) '>>>>      x-interval: ',1.0d0/dble(int_x-1)
      write(*,*) '>>>>      y-interval: ',1.0d0/dble(int_y-1)
      write(*,*) '>>>>      and the level of tolerance is  ',thresh

      call assign_boundary_values(Grid,int_x,int_y)
      error = 1.0d0
      iteration_num = 0
      do while (error .gt. thresh)
         call SOR_iteration(Grid, oldGrid, int_x, int_y, error, omega)
         iteration_num = iteration_num + 1
      enddo

      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
      write(*,*) 'Number of Iterations:                 ',iteration_num
      write(*,*) 'The overRelaxation Parameter (omega): ',omega
      write(*,*) 'The relaxation Error:                 ',error
      write(*,*) 'The tolerance value:                  ',thresh
      write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'

      end
!===================================================================
! This routine implements the Gauss-Seidel overrelaxation scheme
! using the odd-even algorithm. Computing the SOR update afterward.
!===================================================================
      subroutine SOR_iteration(Grid,oldGrid,int_x,int_y,error,omega)
         double precision, dimension(:,:) :: Grid, oldGrid
! HPF$ INHERIT Grid
! HPF$ INHERIT oldGrid
         integer int_x, int_y
         double precision error, omega

!!! local vars ...
      integer k, m

      oldGrid = Grid

!!! Compute Gauss-Seidel update for even points ...
!!! A point is the average of its adjacent north-south and east-west
!!! points ...

      Forall(k=2:int_x-2:2, m=2:int_y-2:2) Grid(k,m) = 
     c   (Grid(k+1,m) + Grid(k-1,m) + Grid(k,m+1) + Grid(k,m-1))/4.0
      Forall(k=3:int_x-1:2, m=3:int_y-1:2) Grid(k,m) = 
     c   (Grid(k+1,m) + Grid(k-1,m) + Grid(k,m+1) + Grid(k,m-1))/4.0

!!! Compute Gauss-Seidel update for odd points ...
!!! A point is the average of its adjacent north-south and east-west
!!! points ...

      Forall(k=3:int_x-1:2, m=2:int_y-2:2) Grid(k,m) = 
     c   (Grid(k+1,m) + Grid(k-1,m) + Grid(k,m+1) + Grid(k,m-1))/4.0
      Forall(k=2:int_x-2:2, m=3:int_y-1:2) Grid(k,m) = 
     c   (Grid(k+1,m) + Grid(k-1,m) + Grid(k,m+1) + Grid(k,m-1))/4.0

!!! Compute the error as the difference between the G-S update and the
!!!  old Grid ...

      error = maxval(abs(Grid - oldGrid))

!!! Now compute the new overRelaxed values then put them back in Grid ...

      Grid = oldGrid + omega * (Grid - oldGrid)

      end subroutine

!=================================================================
! This routine calculates an analytic function, and returns it
! in Uval
!=================================================================
      subroutine compute_val_of_U(Uval, x, y)
      double precision Uval,x,y

!!! local var ...
      double precision part1

        part1 = -4.7 + 0.2 * x - 0.3 * y + 0.09 * (x * x - y * y)
        Uval = part1 + log(sqrt( ((x + 0.2)**2.0) + ((y + 0.2)**2.0) ))

!!!  write(*,*) Uval,x,y
      end subroutine

!===================================================================
! This sub puts the boundary values on the boundary of U
!===================================================================
      subroutine assign_boundary_values(U,int_x,int_y)
!         interface
!            subroutine compute_val_of_U( Uval,x,y)
!                double precision Uval,x,y
!            end subroutine
!         end interface

      double precision, dimension(:,:) :: U
      integer int_x,int_y
      double precision, dimension(1:int_x) :: x_list
      double precision, dimension(1:int_y) :: y_list

! HPF$ INHERIT U

!!! Reset the grid ...

      U = 0.0d0                !! "d" here is for double precision ...

!!! Compute the x's and y's of the grid ...
      do i=1,int_x
         x_list(i) = dble(i-1)/dble(int_x-1)
      enddo
      do i=1,int_y
         y_list(i) = dble(i-1)/dble(int_y-1)
      enddo

!!! Compute those points forming the boundary region ...
!!! write(*,*) '==== Value of U ====', '==== X value ====', '=== Y value ==='

      do i=1,int_x
         call compute_val_of_U( U(i,1) ,x_list(i), 0.0d0)
         call compute_val_of_U( U(i,int_y), x_list(i), 1.0d0)
      enddo
      do i=1,int_y
         call compute_val_of_U( U(1,i), 0.0d0, y_list(i))
         call compute_val_of_U( U(int_x,i), 1.0d0, y_list(i))
      enddo
      end subroutine