C mpfortran -pesize=8 -d_lines jacobi.f
C
C---ROB SALGADO (cps615g)
C   jacobi.f (3.f) 11/08/92 (extended from 11/06/92 PROBLEM 5.3)
C
C   1  4  1
C   4 -20 4
C   1  4  1
C
 
      subroutine showme(U)
      integer, parameter :: R=128,C=64
      real, array(R,C) :: U
      integer i,j
      print *,"------------------BOUNDARY POINTS + SAMPLE OF INTERIOR--"
      print 199, [1:C-1:int((C-1)/7)],C
      do, i=1, R-1, int((R-1)/11)
         print 200,i,U(i,1:C-1:int((C-1)/7)),U(i,C)
      end do
      print 200,R,U(R,1:C-1:int((C-1)/7)),U(R,C)
      print *,"-------------------------------------------------------/"
      print *," "
 199    FORMAT("    ",8I8)
 200      FORMAT(I4,8F8.4)
      end subroutine showme
 
C  Elementary solution of Laplace's equation in 2 dimensions
C  using the red/black SOR iteration.
C  Parameters R, C give array dimensions
 
 
C---AS USUAL, WE WON'T BOTH CHANGING ALL OF THE NAMES.
C
CC  Perform n cycles of SOR on U with relaxation parameter W
CC  Returns updated function, writes error estimate in mxr
CC  Error estimate is max abs change over all n iterations.
C
      function Sor(U, n, W, mxr)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Sor, U
      integer n
      real W, mxr
 
      real, array(R, C) ::  U1, U2
      logical, array(R, C) :: Red
      integer k
      real OnemW                        ! 1 - W
 
C
C      OnemW = 1.0 - W
CC  Test for coord difference i - j being even
C      Red = Iand( SPREAD([1:R], 2, C) - SPREAD([1:C], 1, R),
C     *            1) == 0
C      U1 = U
C      U2 = U
C
 
C-----JACOBI FOR  1  4  1 ---------------------------------------------
C                 4 -20 4
C                 1  4  1
C
      U1=U
      do, k = 1, n
 
      U1(2:R-1, 2:C-1)=(4*( U1(2:R-1, 1:C-2)+U1(2:R-1, 3:C  )
     *                     +U1(1:R-2, 2:C-1)+U1(3:R,   2:C-1))
     *                   +( U1(1:R-2, 1:C-2)+U1(1:R-2, 3:C  )
     *                     +U1(3:R  , 1:C-2)+U1(3:R  , 3:C  )))/20.0
 
 
C       Where ( Red(2:R-1, 2:C-1) )
C        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
C     *                           + U1(3:R  , 2:C-1)
C     *                           + U1(2:R-1, 1:C-2)
C     *                           + U1(2:R-1, 3:C  ) )
C        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
C     *                    + OnemW*U1(2:R-1, 2:C-1)
C       Elsewhere
C        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
C     *                           + U1(3:R  , 2:C-1)
C     *                           + U1(2:R-1, 1:C-2)
C     *                           + U1(2:R-1, 3:C  ) )
C        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
C     *                    + OnemW*U1(2:R-1, 2:C-1)
C       end Where
 
      end do
 
      mxr = maxval(abs(U1 - U)) ! error estimate
      Sor = U1
      end function Sor
 
 
C  Test Redblack SOR
 
      program TestRedblack
 
      INTERFACE
        function Sor(U, n, W, mxr)
        integer, parameter :: R = 128, C = 64
        real, array(R, C) :: Sor, U
CMPF    ONDPU Sor,U
        integer n
        real W, mxr
      END INTERFACE
 
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U0, Z
      integer, array(R, C) :: I, J
      real bsum, eps, err, W
      integer it
 
      I = SPREAD( [1:R], 2, C )         !Coordinates
      J = SPREAD( [1:C], 1, R )
 
C     Z is a solution of the difference equation
 
      Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)
 
      U0 = Z    ! Install boundary values
      print *, "U0=Z"
      call showme(U0)
 
      bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))
 
C     Intialize interior to average on boundary
C     The "-4" prevents counting corners twice
 
      U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)
      print *, "ORIGINAL, UNSOPHISTICATED"
      call showme(U0)
C     Obtain SOR parameter W:
      print *, " W >"
      read *, W
      print *, " it, eps, err"
      it = 0
      eps = 1.0
      do while ( eps > 0.01 )
        U0 = Sor(U0, 1000, W, eps)
        err = maxval(abs(U0 - Z))       ! Actual error
        it = it + 1000
        print *, it, eps, err
        call showme(U0)
      end do
 
      end program TestRedBlack

