C  Eigenvalues of the SOR operator for Laplace's equation
C  in 2 dimensions using the red/black ordering.
C  Parameters R, C give array dimensions

C  Perform n cycles of SOR on approximate eigenvector U
C  with relaxation parameter W, shift S, calculating
C  updated eigenvalue  mxm.
C  Updates eigenvector U, writes eigenvalue in mxm

      subroutine Sor(U, n, W, S, mxm)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U
      integer n
      real W, S, mxm, mxn

      real, array(R, C) ::  U1, U2
      integer, array(R, C) :: I, J
      COMMON /EIGENPE/ I, J
CMPF  ONDPU I, J
D     integer, array(R, C) :: I1, J1

      logical, array(R, C) :: Red
      integer k
D     integer ix, jx
      real OnemW			! 1 - W

      OnemW = 1.0 - W
C  Test for coord difference i - j being even
      Red = Iand(I - J, 1) == 0
      do, k = 1, n
D      print *, "U:"
D      call prU(U)
       U1 = U  ! Initialize to prevent garbage calculations
       U2 = U
       Where ( Red(2:R-1, 2:C-1) )
	U2(2:R-1, 2:C-1) = 0.25*(  U(1:R-2, 2:C-1)
     *				 + U(3:R  , 2:C-1)
     *				 + U(2:R-1, 1:C-2)
     *				 + U(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 - S*U           ! Unnormalized vector
D      print *, "U1:"
D      call prU(U1)
       mxm = maxval(U1)        ! Maximum entry
       mxn = minval(U1)        ! Minimum entry
D      I1 = I
D      J1 = J
D      where ( mxn < U1 .and. U1 < mx1 )
D        I1 = R+100
D        J1 = C+100
D      end where
D      ix = minval(I1)
D      where ( I1 > ix ) J1 = C+100
D      jx = minval(J1) 
D      print *, ix, jx, mxn, mxm
       if ( mxm < -mxn ) mxm = mxn  ! Entry of max abs val
       U = U1/mxm              ! Normalized eigenvector
      end do
D     print *, "U:"
D     call prU(U)
      end subroutine Sor

C  Diagnostic printing
      subroutine prU(U)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: U
CMPF  ONDPU U
      print 21, U(1,1:5)
      print 21, U(2,1:5)
      print 21, U(3,1:5)
      print 21, U(4,1:5)
 21   format (5(1X, E12.5))
      end subroutine prU


C  Look for Redblack SOR eigenvalues

      program EigenRedblack

      INTERFACE
	subroutine Sor(U, n, W, S, mxm)
	integer, parameter :: R = 128, C = 64
	real, array(R, C) :: U
CMPF    ONDPU U
	integer n
	real W, S, mxm
      END INTERFACE

      integer, parameter :: R = 128, C = 64,
     *                      RN = R*C
      real, parameter :: Two31 = 2147483648.0
      real, array(R, C) :: U0
      integer, array(R, C) :: I, J
      COMMON /EIGENPE/ I, J
CMPF  ONDPU I, J
      integer, array(RN) :: RR          ! For random nos.
      real W, S, eps, err, mxm, mxold
      integer it, nt, ilim
      character*80 ibff                 ! Terminal input

CMPF  mpl MPF_RANDOM, MPF_SRANDOM

      I = SPREAD( [1:R], 2, C )		!Coordinates
      J = SPREAD( [1:C], 1, R )
      RR = 1492 + [1:RN]
      call mpf_srandom(RR)

      print *, "Specify nt <= 0 to quit."
      print '(" nt, W, S, eps > " $)'   ! Obtain param's
      read *, nt, W, S, eps
      do while ( nt > 0 )
        call mpf_random(RR)	! Random start
        U0 = RESHAPE(shape=[R, C], source=RR)/Two31
        U0 = 2.0*U0 - 1.0
        where (     I == 1 .or. I == R
     *         .or. J == 1 .or. J == C)
     *    U0 = 0.0	! Boundary
        mxm = maxval(abs(U0))
D       print *, "max(abs(U0)) = ", mxm
        U0 = U0/maxval(abs(U0))     ! normalize
        mxm = 0.0
        it = 0
        ilim = 10*nt
        err = 1.0 + 2.0*eps
 91     print *, "        it       mxm      change"
        do while ( err > eps .and. it < ilim )
          mxold = mxm
	  call Sor(U0, nt, W, S, mxm)
	  err = abs(mxm - mxold)
          it = it + nt
          print *, it, mxm, err
        end do
        print '(" More iterations (y or n) ? " $)'
        read *, ibff
        if ( ibff(1:1) == 'y' ) then
          print '(" nt, S, eps > " $)'
          read *, nt, S, eps
          ilim = ilim + 10*nt
          go to 91
        end if
        print '(" nt, W, S, eps > " $)'   ! Obtain param's
        read *, nt, W, S, eps
      end do

      end program EigenRedBlack
