Refactoring GO TO

Help
1

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
2

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
3

@ivanpribec ,

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
5

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
6

Note that in this statement, the “, save” is redundant since local variables that are initialized in their declaration automatically have the SAVE attribute.

7

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
8

No save > explicit save > implicit save.

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

4 Likes
9

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
10

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.