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

C  General organization of the data:
C
C  Interpolation is defined for x(0) <= x <= x(Nx), though the program
C  may be used to extrapolate, and for non-standard interpolation.
C
C  The function values are stored redundantly:
C
C  Ys(0:Nx-1, -2:3):
C
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   3	|  y( 3)  |  y( 4)  |         |  y(i+3) |         | y(Nx+2) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   2	|  y( 2)  |  y( 3)  |         |  y(i+2) |         | y(Nx+1) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   1	|  y( 1)  |  y( 2)  |         |  y(i+1) |         | y( Nx ) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   0	|  y( 0)  |  y( 1)  |         |  y( i ) |         | y(Nx-1) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C  -1	|  y(-1)  |  y( 0)  |         |  y(i-1) |         | y(Nx-2) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C  -2	|  y(-2)  |  y(-1)  |         |  y(i-2) |         | y(Nx-3) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C	    0         1                    i                  Nx-1
C
C  as are the abscissae (x's)
C
C  Xs(0:Nx-1, -2:3):
C
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   3	|  x( 3)  |  x( 4)  |         |  x(i+3) |         | x(Nx+2) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   2	|  x( 2)  |  x( 3)  |         |  x(i+2) |         | x(Nx+1) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   1	|  x( 1)  |  x( 2)  |         |  x(i+1) |         | x( Nx ) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C   0	|  x( 0)  |  x( 1)  |         |  x( i ) |         | x(Nx-1) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C  -1	|  x(-1)  |  x( 0)  |         |  x(i-1) |         | x(Nx-2) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C  -2	|  x(-2)  |  x(-1)  |         |  x(i-2) |         | x(Nx-3) |
C	+---------+---------+  . . .  +---------+  . . .  +---------+
C	    0         1                    i                  Nx-1
C
C  Values to be interpolated are given by X(0:Nx-1) which normally satisfy
C
C  		x(i) <= X(i) <= x(i+1)
C
C  More particularly, the interpolant for X(i) is obtained using the table
C  defined by Xs(:, i), Ys(:, i) .

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

C  Calculate and return interpolated function values for X based on
C  the table given by Xs, Ys using n-point (order n-1)interpolation.
C  If n is odd (even-order interpolation), r = 0 gives interpolates
C  with lower-numbered points, r = 1 interpolates with higher-numbered
C  points.  If n is even (odd-order interpolation) r = 0 uses symmetrically
C  located points.  Other values of r are normally inappropriate.
C  The program requires 2 <= n <= 6.
C  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
CMF$  LAYOUT X(), Xs(,), Ys(,), Err(), nevilled()
      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
CMF$  LAYOUT A(,)
      END INTERFACE

C  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

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


C -----------------------------------------------------------------------------

C  Set up tables for Neville, using the data given to create a table of
C  Nx subtables, treating X(a) as x(0).  The result is the largest number
C  of points which may be used for interpolation.  If this is 0, there
C  wasn't enough data to define a linear interpolant at both ends.
C  If a >= 2 and M >= Nx + 2  the table is defined for interpolation up
C  to 6 points (quintic), which is the maximum allowed.

      integer function Nevtabl(X, Y, Xs, Ys, M, Nx, a)
      integer M, Nx, a				! Array dimensions, offset
      real, array(0:M-1) :: X, Y		! The data
      real, array(0:Nx-1, -2:3) :: Xs, Ys	! Nevilled tables (result)

      integer d, r

      d = min(a, M - Nx) + 1		! Number of points at ends
      if ( d <= 0 ) then
	d = 0				! Too little to interpolate at ends
	go to 99			! Exit
      end if
C   Fill in tables
      for r = 0, d-1
	Xs(:,  r ) = X(  a-r  : Nx-r-1)
	Ys(:,  r ) = Y(  a-r  : Nx-r-1)
	Xs(:, r+1) = X( a+r+1 :  Nx+r )
	Ys(:, r+1) = Y( a+r+1 :  Nx+r )
      end for
      d = 2*d				! Max number of points allowed
 99   continue				! For exit
      Nevtabl = d
      end function Nevtabl
