I would like to refactor the following elegant piece of code from Forsyth, Malcolm, & Moler (1977), Computer Methods for Mathematical Computation, Prentice-Hall, Inc. (available as fmm on Netlib):
double precision function seval(n, u, x, y, b, c, d) integer n double precision u, x(n), y(n), b(n), c(n), d(n) c c this subroutine evaluates the cubic spline function c c seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3 c c where x(i) .lt. u .lt. x(i+1), using horner's rule c c if u .lt. x(1) then i = 1 is used. c if u .ge. x(n) then i = n is used. c c input.. c c n = the number of data points c u = the abscissa at which the spline is to be evaluated c x,y = the arrays of data abscissas and ordinates c b,c,d = arrays of spline coefficients computed by spline c c if u is not in the same interval as the previous call, then a c binary search is performed to determine the proper interval. c integer i, j, k double precision dx data i/1/ if ( i .ge. n ) i = 1 if ( u .lt. x(i) ) go to 10 if ( u .le. x(i+1) ) go to 30 c c binary search c 10 i = 1 j = n+1 20 k = (i+j)/2 if ( u .lt. x(k) ) j = k if ( u .ge. x(k) ) i = k if ( j .gt. i+1 ) go to 20 c c evaluate spline c 30 dx = u - x(i) seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) return end
The nice thing about the function is it keeps the number of comparisons and evaluations needed to find the right interval of the spline to a minimum.
The translation I have come up with is
real(dp) function seval(n, u, x, y, b, c, d) integer, intent(in) :: n real(dp), intent(in) :: u real(dp), intent(in), dimension(n) :: x, y, b, c, d integer, save :: i = 1 real(dp) :: dx if ( i >= n ) i = 1 if ( u < x(i) ) then call binary_search(i) else if ( u > x(i + 1) ) then call binary_search(i) end if dx = u - x(i) seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) contains subroutine binary_search(i) integer, intent(out) :: i integer :: j, k i = 1 j = n+1 do k = (i+j)/2 if ( u .lt. x(k) ) j = k if ( u .ge. x(k) ) i = k if ( j .le. i+1) exit end do end subroutine end function
I was tempted to introduce an
.or. to remove one of the statements calling the binary search, but I believe the
else if is closer to the original logic.
Would you do anything differently?