C  Solution of Laplace's equation in 2 dimensions
C  using the point-Jacobi iteration, 9-point difference op.
C  Parameters R, C give array dimensions

C  Perform n cycles of point-Jacobi on U
C  Returns updated function, writes error estimate in mxr
C  Error estimate is max change over n iterations.

      function Jacobi9(U, n, mxr)
      integer, parameter :: R = 128, C = 64
      real, array(R, C) :: Jacobi9, U
      integer n
      real mxr

      real, array(R, C) ::  W1, Wr
      integer k

C   Wr = * + *

      W1 = U
      Wr = U
      do, k = 1, n
        Wr(:, 2:C-1) = W1(:, 1:C-2) + W1(:, 3:C)
	W1(2:R-1, 2:C-1) =
     *	    0.20*(  Wr(2:R-1 ,2:C-1)
     *            + W1(1:R-2, 2:C-1) + W1(3:R, 2:C-1))
     *	  + 0.05*(  Wr(1:R-2, 2:C-1) + Wr(3:R, 2:C-1))
      end do
      mxr = maxval(abs(W1 - U))	! error estimate
      Jacobi9 = W1
      end function Jacobi9


C  Test Jacobi

      program TestJacobi

      INTERFACE
	function Jacobi9(U, n, mxr)
	integer, parameter :: R = 128, C = 64
	real, array(R, C) :: Jacobi9, U
CMPF	ONDPU U
	integer n
	real 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
      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

      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)

      eps = 1.0
      print *, " it, eps, err"
      do while ( eps > 0.01 )
D     print *, " U0"
D     print *, U0(1, 1:4)
D     print *, U0(2, 1:4)
D     print *, U0(3, 1:4)
	U0 = Jacobi9(U0, 1000, eps)
        err = maxval(abs(U0 -  Z))
        it = it + 1000
        print *, it, eps, err
      end do

      end program TestJacobi
