C  Solution of Laplace's equation in 2 dimensions
C  using  a parallel multi-grid approach with the
C  point-Jacobi iteration, 9-point difference op.
C  Parameters R, C give array dimensions
C  This has routines for managing experiments,
C  but no main program.
C
C  File lapm9u.f
C
C  Parameters are in lapm9.hp.  These are:
C      R, C  -  Rows, cols for arrays

C  Perform n cycles of point-Jacobi on U with grid
C  spacing q, which is <= min(R, C) - 1.
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max change over n iterations.
C  Updates points no closer than q to boundaries.

      function Jacobim9(U, n, q, mxr)
      INCLUDE 'lapm9.h'
      real, array(R, C) :: Jacobim9, U
      integer n, q
      real mxr

      real, array(R, C) ::  W1, Wr
      integer k, ra, rb, ca, cb

      ra = 1 + q    ! Update region is (ra:rb, ca:cb)
      rb = R - q
      ca = 1 + q
      cb = C - q

C   Wr = * ... + ... *

      W1 = U
      Wr = U
      do, k = 1, n
        Wr(:, ca:cb) =
     *      W1(:, ca+q:cb+q) + W1(:, ca-q:cb-q)
	W1(ra:rb, ca:cb) =
     *	    0.20*(  Wr(ra:rb, ca:cb)
     *            + W1(ra+q:rb+q, ca:cb)
     *            + W1(ra-q:rb-q, ca:cb))
     *	  + 0.05*(  Wr(ra+q:rb+q, ca:cb)
     *            + Wr(ra-q:rb-q, ca:cb))
      end do
      mxr = maxval(abs(W1 - U))	! error estimate
      Jacobim9 = W1
      end function Jacobim9

C  Intialize interior by using 5-point difference
C  operator with unequal spacing from boundary.
C  Z provides boundary values, I, J are coordinate
C  arrays, all arrays are r by c (arguments, not
C  system parameters here).

      function initlap(Z, I, J, r, c)
      integer r, c
      real, array(r, c) ::  initlap, Z
      integer, array(R, c) :: I, J

      real, array(r, c) :: U, Rx, Cx

      Rx = (I-1)*(r-I)
      Cx = (J-1)*(c-J)

      U =  Cx*( (r-I)*SPREAD(Z(1,:), 1, r)
     *        + (I-1)*SPREAD(Z(r,:), 1, r))/(r-1)
     *   + Rx*( (c-J)*SPREAD(Z(:,1), 2, c)
     *        + (J-1)*SPREAD(Z(:,c), 2, c))/(c-1)
 
      where (     I == 1 .or. I == r
     *       .or. J == 1 .or. J == c )
        U = Z           ! Install boundary values
      elsewhere
        U = U/(Rx + Cx) ! Normalize
      end where
      initlap = U
      end function initlap


C  Interactive solution of R by C problems using
C  Jacobi or parallel W, V cycles

      subroutine TestJacobim(U0)
      INCLUDE 'lapm9.hp'
      real, array(R, C) :: U0
CMPF  ONDPU U0

      INTERFACE
	function Jacobim9(U, n, q, mxr)
        INCLUDE 'lapm9.hp'
	real, array(R, C) :: Jacobim9, U
CMPF	ONDPU U
	integer n, q
	real mxr
      END INTERFACE
      INTERFACE
        function V(U, n, p, lim, tn, mxr)
        INCLUDE 'lapm9.hp'
        real, array(R, C) :: U, V
CMPF    ONDPU U
        integer n, p, lim, tn
        real mxr
      END INTERFACE
      INTERFACE
        function W(U, n, p, lim, tn, mxr)
        INCLUDE 'lapm9.hp'
        real, array(R, C) :: U, W
CMPF    ONDPU U
        integer n, p, lim, tn
        real mxr
      END INTERFACE

      real, array(R, C) :: Us, Zr
      real eps, err
      integer iter, titer, eltm, p, q
      character*50 buff
      logical more

      print 21, R, C
 21   FORMAT (" Grid size is (1:", I3, ", 1:", I3, ")")
      Us = U0
      Zr = U0
      more = .TRUE.
      do while ( more )
        print '(" J, V, W > " $)'
        read *, buff
D       print *, " U0"
D       print *, U0(1, 61:64)
D       print *, U0(2, 61:64)
D       print *, U0(3, 61:64)
D       print *, U0(4, 61:64)
        if ( buff(1:1) == 'J' ) then
          print '(" q, iter > " $)'
          read *, q, iter
          call MpTimerStart()
	  U0 = Jacobim9(U0, iter, q, eps)
          eltm = MpTimerElapsed()
          titer = titer + iter
        else if ( buff(1:1) == 'V' ) then
          print '(" p, lim, iter > " $)'
          read *, p, lim, iter
          call MpTimerStart()
          U0 = V(U0, iter, p, lim, titer, eps)
          eltm = MpTimerElapsed()
        else if ( buff(1:1) == 'W' ) then
          print '(" p, lim, iter > " $)'
          read *, p, lim, iter
          call MpTimerStart()
          U0 = W(U0, iter, p, lim, titer, eps)
          eltm = MpTimerElapsed()
        else
          print *, "Inappropriate input"
          go to 88
        end if
        err = maxval(abs(U0 - Zr))
        print 22, titer, eps, err,  eltm
 22     FORMAT (1X, I6, " iterations,  change = ",
     *          E12.5, "  error = ", E12.5/
     *          I6, " msec.")
        print '(" M(ore), R(ef), S(tart), Q > " $)'
        buff(1:2) = '  '   ! Clear buff
        read *, buff
        if ( buff(1:1) == 'R' ) Zr = U0
        if ( buff(1:1) == 'S' .or. buff(2:2) == 'S' ) then
          U0 = Us
          titer = 0
        end if
        if ( buff(1:1) == 'Q' ) more = .FALSE.
 88     continue
      end do
 99   continue    ! For exit
      end subroutine TestJacobim
