C mpfortran -pesize=8 -d_lines 2.f
C
C----ROB SALGADO (cps615g)
C    2.F 11/12/92 (extended from 11/06/92 PROBLEM 5.2) 
C
C    MORE-SOPHISTICATED INITIAL-VALES
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  Perform n cycles of SOR on U with relaxation parameter W
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max abs change over all n iterations.

      function Sor(U, n, W, offset)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Sor, U
      integer n
      real W, offset
      real maxU1, maxU0
      COMMON maxU1, maxU0

      real, array(R, C) ::  U1, U2, UU
      logical, array(R, C) :: Red
      integer k
      real OnemW                        ! 1 - W

      OnemW = 1.0 - W
C  Test for coord difference i - j being even
      Red = Iand( SPREAD([1:R], 2, C) - SPREAD([1:C], 1, R),
     *            1) == 0
      U1 = U
      U2 = U
      do, k = 1, n
         UU=U
       Where ( Red(2:R-1, 2:C-1) )
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       Elsewhere
        U2(2:R-1, 2:C-1) = 0.25*(  U1(1:R-2, 2:C-1)
     *                           + U1(3:R  , 2:C-1)
     *                           + U1(2:R-1, 1:C-2)
     *                           + U1(2:R-1, 3:C  ) )
        U1(2:R-1, 2:C-1) =      W*U2(2:R-1, 2:C-1)
     *                    + OnemW*U1(2:R-1, 2:C-1)
       end Where
      U1=U1-offset*UU
       maxU0=maxU1
       maxU1 = maxval(abs(U1))  ! maximum absolute value obtained
C  Re-normalize to the maximum value
      U1= U1/maxU1
      end do
C  Re-normalize to the maximum value
      Sor = U1
      print *,"Current:                    ",maxU1
      end function Sor


            
C  Test Redblack SOR

      program TestRedblack

      INTERFACE
        function Sor(U, n, W, offset)
        integer, parameter :: R = 128, C = 64
        real, array(R, C) :: Sor, U
CMPF    ONDPU Sor,U
        integer n
        real W, offset
      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, epsx, offset, maxU0, maxU1
      integer it, nits, notsure
      COMMON maxU1, maxU0
C--Ask for number of iterations
C
      print *, "How many NITS? "
      read *,nits

      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
C      print  *, "ORIGINAL ANALYTIC SOLUTION... ZERO BCs IMPOSED LATER"
C      call showme(U0)

       bsum = SUM(U0(1,:)) + SUM(U0(R,:)) ! Sum on boundary
     *     + SUM(U0(:,1)) + SUM(U0(:,C))
C
CC     Intialize interior to average on boundary
CC     The "-4" prevents counting corners twice
C
       U0(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)
C       print *,"INITIAL GUESS: UNSOPHISTICATED"
C       call showme(U0)
      Z(1,1:C)=0.0
      Z(R,1:C)=0.0
      Z(1:R,1)=0.0
      Z(1:R,C)=0.0
      U0(1,1:C)=0.0
      U0(R,1:C)=0.0
      U0(1:R,1)=0.0
      U0(1:R,C)=0.0

C     Obtain SOR parameter W:
      print *, "SOR param: W >"
      read *, W
      print *, "epsilon: >"
      read *, epsx
      print *, "offset: >"
      read *,offset
      print *, "        ITS, EPS"
      it = 0
      notsure=0
      do while (notsure < 2)
      eps=1.0
       do while ( (abs(eps) > epsx))
        U0 = Sor(U0, nits, W, offset)    !Number of ITerationS=nits
         eps=maxU1-maxU0
        it = it + nits
        print *, it, eps
D        call showme(U0)
       end do
      notsure=notsure+1
      end do

      end program TestRedBlack

