!  Simple test of low-order polynomial interpolation using Neville's
!  algorithm with differences.  Test data are log10 values tabulated
!  at intervals designed to give roughly equal errors in all parts of
!  the table.  The program permits interpolation up to order 5.
!
!  Files:  neville.fcm  neville.h

!  General organization of the data:
!
!  Interpolation is defined for x(0) <= x <= x(Nx) .
!
!  The function values are stored redundantly:
!
!  Ys(0:Nx-1, -2:3):
!
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   3	|  y( 3)  |  y( 4)  |         |  y(i+3) |         | y(Nx+2) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   2	|  y( 2)  |  y( 3)  |         |  y(i+2) |         | y(Nx+1) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   1	|  y( 1)  |  y( 2)  |         |  y(i+1) |         | y( Nx ) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   0	|  y( 0)  |  y( 1)  |         |  y( i ) |         | y(Nx-1) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!  -1	|  y(-1)  |  y( 0)  |         |  y(i-1) |         | y(Nx-2) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!  -2	|  y(-2)  |  y(-1)  |         |  y(i-2) |         | y(Nx-3) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!	    0         1                    i                  Nx-1
!
!  as are the abscissae (x's)
!
!  Xs(0:Nx-1, -2:3):
!
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   3	|  x( 3)  |  x( 4)  |         |  x(i+3) |         | x(Nx+2) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   2	|  x( 2)  |  x( 3)  |         |  x(i+2) |         | x(Nx+1) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   1	|  x( 1)  |  x( 2)  |         |  x(i+1) |         | x( Nx ) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!   0	|  x( 0)  |  x( 1)  |         |  x( i ) |         | x(Nx-1) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!  -1	|  x(-1)  |  x( 0)  |         |  x(i-1) |         | x(Nx-2) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!  -2	|  x(-2)  |  x(-1)  |         |  x(i-2) |         | x(Nx-3) |
!	+---------+---------+  . . .  +---------+  . . .  +---------+
!	    0         1                    i                  Nx-1
!
!  Values to be interpolated are given by X(0:Nx-1) which should satisfy
!
!  		x(i) <= X(i) <= x(i+1)

!  ------ function nevilled(X, Xs, Ys, Err, n, r) --------------------------

!  Calculate and return interpolated function values for X based on
!  the table given by Xs, Ys using n-point (order n-1) interpolation.
!  If n is odd (even-order interpolation), r = 0 gives interpolates
!  with lower-numbered points, r = 1 interpolates with higher-numbered
!  points.  If n is even (odd-order interpolation) r = 0 uses symmetrically
!  located points.  Other values of r are normally inappropriate.
!  The program requires 2 <= n <= 6.
!  Err is assigned an error estimate for each interpolant.

      function nevilled(X, Xs, Ys, Err, n, r)
      INCLUDE 'neville.h'
      real, array(0:Nx-1) :: X, Err, nevilled
      real, array(0:Nx-1, -2:3) :: Xs, Ys
      integer n, r

      INTERFACE
      subroutine prmf(A, ra, rb, ca, cb, Nx, Ny)
      integer Nx, Ny
      real, array(Nx, Ny) :: A
      integer ra, rb, ca, cb
      END INTERFACE

!  Cs(u, j, k) is C[j:j+k] for interpolating X(u), similarly Ds
      real, array(0:Nx-1, -2:2, 1:5) :: Cs, Ds		! Differences
      real, array(0:Nx-1) :: Z, W			! Working storage
      integer i1, i2, j, k, u, k1, iin, iix

      i1 = -(n - r - 1)/2		! Lowest point in interpolation
      i2 = i1 + n - 2			! Highest j for C[j:j+1]
!     print *, "X(4):", X(4)
!     print *, "Xs(4, 0:1)", Xs(4, 0:1)
!     print *, "Ys(4, 0:1)", Ys(4, 0:1)
!     Cs = 0.0				! So as not to print garbage
!     Ds = 0.0				! from unused entries.
!     print *, "i1, i2:", i1, i2
!  Calculate linear C's, D's
      do, j = i1, i2
	W  = (Ys(:, j+1) - Ys(:, j))/(Xs(:, j+1) - Xs(:, j))
	Cs(:, j, 1) = (X  -  Xs(:, j))*W
	Ds(:, j, 1) = (X - Xs(:, j+1))*W
      end do ! j = i1, i1
!  Calculate higher-order C's, D's
      k = 1				! Previous entries are C,D[j:j+k]
      do while ( k+1 < n )
!       Calculate C,D[j:j+k+1]
	i2 = i2 - 1			! Highest j for C[j:j+k+1]
	k1 = k + 1
!	print *, "Level ", k1, i1, i2
	do, j = i1, i2
	  u = j + k + 1			! Calculate C[j:u]
	  W = (Cs(:, j+1, k) - Ds(:, j, k))/(Xs(:, u) - Xs(:, j))
	  Cs(:, j, k1) = (X - Xs(:, j))*W
	  Ds(:, j, k1) = (X - Xs(:, u))*W
	end do ! j = i1, i2
	k = k1				! Now we have C,D[j:j+k]'s
      end do ! while ( k+1 < n )
!  Now Cs, Ds(:, i1, k)	are C,D[i1:i1+k]
!  So  Cs, Ds(:, I1, k) are C,D[i1:i2]
!     print *, "Cs:"
!     call prmf(Cs(4,:,:), 1, 5, 1, 5, 5, 5)
!     print *, "Ds:"
!     call prmf(Ds(4,:,:), 1, 5, 1, 5, 5, 5)
      i2 = i1 + k
!  Form error estimate = smaller of final C, D:
      Err = Cs(:, i1, k)
      where ( abs(Err) > abs(Ds(:, i1, k)) ) Err = Ds(:, i1, k)
!  Sum differences back to data, alternating C/D:
      Z = 0.0				! To accumulate sum
      if ( iand(n, 1) == 0 ) then	! If n is even
	Z  = Cs(:, i1, k)		! Begin with C
!	print *, "C[", i1, ":", i2, "]"
	i2 = i2 - 1
	k  = k - 1
      end if ! n is even
!  Finish sum with D + C + . . . + D + C
!  Next term to add is D[i1:i2]
      do while ( k >= 2 )
	Z  = Z + Ds(:, i1, k)
!	print *, "D[", i1, ":", i2, "]"
	i1 = i1 + 1
	k  = k - 1
	Z  = Z + Cs(:, i1, k)
!	print *, "C[", i1, ":", i2, "]"
	i2 = i2 - 1
	k  = k - 1
      end do ! while ( k >= 2 )
      nevilled = Ys(:, i1) + Z		! Final result
      end function nevilled


!  Print real A(ra:rb, ca:cb), at most 6 columns

      subroutine prmf(A, ra, rb, ca, cb, Nx, Ny)
      integer Nx, Ny
      real, array(Nx, Ny) :: A
      integer ra, rb, ca, cb
      integer ir

      if ( cb >= ca+6 ) cb = ca+5
      do, ir = ra, rb
       print 21, A(ir, ca:cb)
 21    FORMAT ( 6(1X, E11.4:) )
      end do
      end subroutine prmf

!  ------ subroutine nreport(f, X, Y, E) -----------------------------------

!  Analyse and report results of interpolation.

      subroutine nreport(f, X, Y, E)
      include 'neville.h'
      real f				! X = X0 + f*D
      real, array(0:Nx-1) :: X, Y, E

      real, array(0:Nx-1) :: Yc, Z, R
      integer, array(1) :: in, ix	! For maxloc, minloc
      integer iin, iix

      Yc = log10(X)		! Correct logs (nearly)
      Z  = Yc - Y		! Actual errors
      R  = abs(E/Z)		! Ratios |estimate/actual|
!     print *, "X ",  X(1:5)
!     print *, "Yc", Yc(1:5)
!     print *, "Y ",  Y(1:5)
      in = minloc(Z)
      ix = maxloc(Z)
      iin = in(1) - 1		! minloc, maxloc assume subscripts
      iix = ix(1) - 1		! start from 1
      print 21, f, Z(iin), iin, Z(iix), iix, maxval(abs(Z/Yc)), &
                minval(R), maxval(R)
 21   FORMAT (X, F4.2, 2(X, E11.4, X, I4), 3(X, E11.4))
      end subroutine nreport

!  ------ program nevtest --------------------------------------------------

!  Test using logarithms (base 10) tabulated at unequal intervals.
!  These might not be optimum intervals.

      program nevtest
      INCLUDE 'neville.h'

      INTERFACE
      function nevilled(X, Xs, Ys, Err, n, r)
      INCLUDE 'neville.h'
      real, array(0:Nx-1) :: X, Err, nevilled
      real, array(0:Nx-1, -2:3) :: Xs, Ys
      integer n, r
      END INTERFACE

      INTERFACE
      subroutine nreport(f, X, Y, E)
      include 'neville.h'
      real f				! X = X0 + f*D
      real, array(0:Nx-1) :: X, Y, E
      END INTERFACE

      INTERFACE
      subroutine prmf(A, ra, rb, ca, cb, Nx, Ny)
      integer Nx, Ny
      real, array(Nx, Ny) :: A
      integer ra, rb, ca, cb
      END INTERFACE

      real,    array(0:Nx-1, -2:3) :: XX, YY	! Table
      integer, array(0:Nx-1, -2:3) :: Is, Js	! Coordinates
      real,    array(0:Nx-1) :: X, Y, D, E	! For interpolation
      integer, array(0:Nx-1) :: I		! Coordinates
      real h
      integer ix, n, r, ii

      I = (/ (ii, ii=0,Nx-1) /)
      Is = SPREAD(I, 2, 6)
      Js = SPREAD([-2:3], 1, Nx)
      print *, "Interpolation with log10, specify h <= 0 to exit."
      do while ( .TRUE.)
	print '(" h > "$)'
	read *, h
	if ( h <= 0.0 ) go to 99	! exit
!  Construct the table (the easy way):
        XX = (1.0 + h)**(Is + Js)	! Non-uniform spacing
!	XX = 1.0 + h*(Is + Js)		! Uniform spacing for testing
	print *,  "Range:",  XX(0, 0), XX(Nx-1, 1)
        YY = log10(XX)
	D  = 0.25*(XX(:, 1) - XX(:, 0))	! Differences for test values
!       print *, "XX:"
!	call prmf(XX, 1, 6, 1, 6, Nx, 6)
!       print *, "YY:"
!	call prmf(YY, 1, 6, 1, 6, Nx, 6)
!  Try various interpolations:
	print *, "Interpolations, n points with offset r."
	print *, "Specify n <= 0 to quit."
	do while ( .TRUE. )
	  print '(" n, r > "$)'
	  read *, n, r
	  if ( n <= 0 ) go to 88
	  if ( n > 6 ) then		! Bad
	    print *, "The program requires n <= 6."
	  else	! Conduct the test
	    X  = XX(:, 0)		! Initialize X
	    print 22
 22   FORMAT("   f   min error   at   max error   at  max rel err", &
             "  min ratio   max ratio")
	    do, ix = 1, 3
	      X = X + D
	      Y = nevilled(X, XX, YY, E, n, r)
	      call nreport(ix*0.25, X, Y, E)
	    end do ! ix = 1, 3
	  end if ! n > 6
	end do ! Various n, r combinations
 88	continue ! For exit
      end do ! Various values of h
 99   continue				! For exit
      end program nevtest
