# Refactoring GO TO

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?

2 Likes

I would write this as:

``````      if (u < x(k)) then
j = k
else
i = k
end if
if (j <= i+1) exit
``````

I would remove an extra space in `real(dp), intent(in) :: u` before `u`.

The other conceptual issue is that I would not `save` `i`, but rather return it to the caller to handle. Also, in practice, the interval is typically the same as last time, or the one next to it to the right. So I have used the following logic:

Which I believe is faster when I benchmarked it.

3 Likes

Building on the `seval` procedural changes by @certik, a suggestion will be to refactor as a parameterized derived type, `spline_t`, in a `container` approach that a consumer can then employ conveniently for the floating-point precision of their choice.

``````module spline_m
use kinds_m, only : RK1, RK2, .. ! supported real kinds
type :: spline_t(K,N)
integer, kind :: K
integer, len  :: N
real(kind=K) :: x(N)
real(kind=K) :: y(N)
real(kind=K) :: b(N)
real(kind=K) :: d(N)
real(kind=K) :: d(N)
integer :: idx !<-- interval index
end type
..
generic :: CalcSpline => CalcSpline_RK1, CalcSpline_RK2, ..
generic :: EvalSpline => EvalSpline_RK1, EvalSpline_RK2, ..
contains
elemental subroutine CalcSpline_RK1( spline )
! Argument List
type(spline_t(K=RK1,N=*)), intent(inout) :: spline
..
include "CalcSpline.i90"

elemental subroutine CalcSpline_RK2( spline )
! Argument List
type(spline_t(K=RK2,N=*)), intent(inout) :: spline
..
include "CalcSpline.i90"

elemental function EvalSpline_RK1( spline, u ) result(r)
! Argument List
type(spline_t(K=RK1,N=*)), intent(in) :: spline
real(kind=spline%K), intent(in)       :: u
! Function result
real(kind=spline%K) :: r
..

end module
``````

Here the interval index is effectively handled by the caller and the numerical method then does not need to use a static object for its workflow which, as you know, is not thread-safe and which hinders parallelization.

A further consideration can be the procedures `CalcSpline`, `EvalSpline`, etc. to be made `type-bound` for even more convenient usage by consumers of said numerical containers. But some folks may not like the fact the passed-object dummy argument then becomes polymorphic per the current Fortran standard. Should Fortran ever allow the option to derive an inextensible derived type, this wonâ€™t be an issue of course.

1 Like

Thanks @certik for the suggestions. Speaking of the `save` statement, I thought it solves precisely the problem of remaining in the same interval, with the binary search performed only when we move into a new interval. At least in the context of a serial program, I think itâ€™s fine.

@FortranFan, Iâ€™ve written such containers before but I donâ€™t think they are always worth the effort. In my specific context I am already encapsulating the spline data in a derived type representing an axisymmetric drying problem. Since I am using some custom spline boundary conditions, I just implemented these at the level of the `b`, `c`, `d` arrays directly.

In the context of a general fitting library, I do like the abstractions that MATLAB uses for splines and other types of piece-wise polynomial interpolants. A piece-wise polynomial is created with the call

``````% create a piece-wise polynomial
pp = spline(x,y)  % other valid interpolants are pchip and makima

% evaluate the spline at points xq
v = ppval(pp,xq)
``````

The `pp` object is just a struct with the fields:

``````pp = struct with fields:
form: 'pp'
breaks: [0 4 10 15]
coefs: [3x5 double]
pieces: 3
order: 5
dim: 1
``````

The `pp` struct is easy to integrate or differentiate due to the coefficient representation of the polynomials. Iâ€™d like to see something similar as part of `stdlib` in the future.

Some of the descriptions of the MATLAB spline fitting functions mention they are in fact just (re)implementations (or maybe even wrappers) of the spline-fitting tools in the Fortran library PPPACK from De Boor. The same goes for the Scipy spline fitting functions, which are wrappers of the FITPACK library.

If only there existed a Fortran endowment which would collect royalties each time people use spline curvesâ€¦

1 Like

Note that in this statement, the â€ś, saveâ€ť is redundant since local variables that are initialized in their declaration automatically have the SAVE attribute.

I believe itâ€™s been mentioned a few times before, that being explicit and adding the `save` attribute serves as a reminder the behavior is intentional.

3 Likes

No save > explicit save > implicit save.

My rule of thumb, but I think that most agree with it.

4 Likes

Yes, thatâ€™s sad Indeed. But every conformant compiler shall provide 2 REAL kinds at least, one toward `selected_real_kind( p=6 )` and another for `selected_real_kind( p=12 )`, so a library author can cover part of the ground with floating-point generics with their math libraries. Agree totally on the use case for better generics though, the good news is the message has been received and heard.

As a fan of parameterized derived types, Iâ€™ve been copying an eye on the compilers and based on the last I had tried out, Iâ€™m fairly certain â€śthat codeâ€ť I posted upthread will work with 3 compilers: IFORT as part of Intel oneAPI, NAG Fortran 7.0, gfortran 10.1.x (Experimental) (Iâ€™ve not been able to test Cray, I think itâ€™s a 4th supporting compiler for such code). The outstanding problem with some of the compilers is the support toward `polymorphic` objects of a declared-type that has a length-type parameter; this can come to â€śbiteâ€ť with type-bound procedures. But I didnâ€™t suggest TBPs in â€śthat codeâ€ť.

Thatâ€™s fairly decent compiler support toward PDTs that Fortranners should start considering the length-type parameter in user derived types. It can then prove a good push to iron out the remaining bugs in implementations and to get the other compilers to cover missing ground in their Fortran 2003 support.

1 Like

Fair enough: when code reviews that enforce â€śseparation of concerns (SoC)â€ť are not an issue, one has considerable flexibility with encapsulation. With SoC in some domains, including a spline-fitting scheme and its data and algorithm(s) into an encapsulation, say for a drying problem calculation/simulation, will be a no-no.