C  Solution of Laplace's equation in 2 dimensions
C  using  a parallel multi-grid approach with
C  point-Jacobi iteration, 9-point difference op.,
C  and using two array sizes.
C  
C  File biglap.f
C
C  Parameters are in biglap.hp, which automatically
C  includes lapm9.hp,
C  These are:
C      Rg,Cg -  Rows, cols for big arrays
C      R, C  -  Rows, cols for small arrays
C  This version has Rg = 2*R - 1, Cg = 2*C - 1.
C
C  Large array is Restricted to smaller size for
C  parallel multi-grid solution, which result is
C  Prolonged back to larger size.
C
C  Uses programs in lapm9u.f for small array solution

C  Perform n cycles of point-Jacobi on U (with grid
C  spacing 1).
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max change over n iterations.
C  The local parameters were used to adapt this from
C  the multi-grid version.

      function Jacobi9(U, n, mxr)
      INCLUDE 'biglap.hp'
      real, array(Rg, Cg) :: Jacobi9, U
      integer n
      real mxr

      real, array(Rg, Cg) ::  W1, Wr
      integer k
      integer, parameter :: ra = 2, rb = Rg - 1
      integer, parameter :: ca = 2, cb = Cg - 1
      integer, parameter ::  q = 1

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
      Jacobi9 = W1
      end function Jacobi9


C  Full weighting restriction operator.  Boundary
C  entries use all the neighbors they have.

      function Restrict(U)
      INCLUDE 'biglap.hp'
      real U(Rg, Cg), Restrict(R, C)
CMPF  ONDPU U

      real, array(R, C) :: X, Y   ! Working arrays
      
      Restrict = 0.25*U(::2, ::2) ! Centers
      X = 0.0
      Y(:, 1:C-1) = U(::2, 2::2)  ! Same row
      X(:, 1:C-1) = Y(:, 1:C-1)
      X(:, 2:C)   = X(:, 2:C)   + Y(:, 1:C-1)
      Y(1:R-1, :) = U(2::2, ::2)  ! Same col
      X(1:R-1, :) = X(1:R-1, :) + Y(1:R-1, :)
      X(2:R, :)   = X(2:R, :)   + Y(1:R-1, :)
      Restrict    = Restrict + 0.125*X
      X = 0.0
      Y(1:R-1, 1:C-1) = U(2::2, 2::2)  ! Corners
      X(1:R-1, 1:C-1) = Y(1:R-1, 1:C-1)
      X(1:R-1, 2:C  ) = X(1:R-1, 2:C  ) + Y(1:R-1, 1:C-1)
      X(2:R,   1:C-1) = X(2:R,   1:C-1) + Y(1:R-1, 1:C-1)
      X(2:R,   2:C  ) = X(2:R,   2:C  ) + Y(1:R-1, 1:C-1)
      Restrict = Restrict + 0.0625*X
C  Adjust sides and corners
      Restrict(2:R-1, 1:C:C-1)
     *   = 1.3333333*Restrict(2:R-1, 1:C:C-1)
      Restrict(1:R:R-1, 2:C-1)
     *   = 1.3333333*Restrict(1:R:R-1, 2:C-1)
      Restrict(1:R:R-1, 1:C:C-1)
     *   = 1.7777777*Restrict(1:R:R-1, 1:C:C-1)
      end function Restrict


C  Prolong with bilinear interpolation.  Z gives
C  boundary values.

      function Prolong(U, Z)
      INCLUDE 'biglap.hp'
      real Prolong(Rg, Cg), U(R, C), Z(Rg, Cg)

      real X(Rg, Cg)      ! Working array

      Prolong = Z         ! Exact boundary
      Prolong(2:Rg-1:2, 2:Cg-1:2) = 0.0
      X(::2, ::2) = U     
      Prolong(::2, ::2) = X   ! Centers

C  Test Big Laplace system

      program TestBigLap
      INCLUDE 'biglap.hp'

      INTERFACE
	function Jacobi9(U, n, mxr)
        INCLUDE 'biglap.hp'
	real, array(Rg, Cg) :: Jacobim9, U
CMPF	ONDPU U
	integer n
	real mxr
      END INTERFACE
      INTERFACE
        function initlap(Z, I, J, r, c)
        integer r, c
        real, array(r, c) :: initlap, Z
        integer, array(r, c) :: I, J
CMPF    ONDPU Z, I, J
      END INTERFACE
      INTERFACE
      subroutine TestJacobim(U0)
        INCLUDE 'lapm9.hp'
        real, array(R, C) :: U0
CMPF    ONDPU U0
      END INTERFACE


      real, array(Rg, Cg) :: U0, Z, Zr
      real, array(R,  C ) :: Ux
      integer, array(Rg, Cg) :: I, J
      real eps, err, trc
      integer iter, titer, eltm, p, q
      character*50 buff
      logical more

      I = SPREAD( [1:Rg], 2, Cg)	! Coordinates
      J = SPREAD( [1:Cg], 1, Rg)

C     Z is a solution of the differential equation but
C     not of the difference equation.

      Z = -51.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)
     *   + log( (I - 129.0)**2 + (J + 1.0)**2 )
      Zr = Z                 ! Initial reference sol'n
      U0 = initlap(Z, I, J, Rg, Cg)

      print 21, Rg, Cg
 21   FORMAT (" Large grid is (1:", I3, ", 1:", I3, ")")
      more = .TRUE.
      do while ( more )
        print '(" J, R> " $)'
        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 '(" iter > " $)'
          read *, iter
          call MpTimerStart()
	  U0 = Jacobi9(U0, iter, eps)
          eltm = MpTimerElapsed()
          titer = titer + iter
        else if ( buff(1:1) == 'R' ) then
          print *, "Restricting . . . ."
          Ux = Restrict(U0)
          call TestJacobim(Ux)
          print *, "Prolonging . . . ."
          U0 = Prolong(Ux)
        else
          print *, "Inappropriate input"
          go to 88
        end if
        err = maxval(abs(U0 - Zr))
        trc = maxval(abs(U0 - Z ))
        print 22, titer, eps, err, trc, eltm
 22     FORMAT (1X, I6, " iterations,  change = ",
     *          E12.5, "  error = ", E12.5/
     &          "  truncation = ", 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 = initlap(Z, I, J, Rg, Cg)
          titer = 0
        end if
        if ( buff(1:1) == 'Q' ) more = .FALSE.
 88     continue
      end do
 99   continue    ! For exit
      end program TestBigLap
