Achieving numpy like fast interpolation in Fortran

Question

I have a numerical routine that I need to run to solve a certain equation, which contains a few nested four loops. I initially wrote this routine into Python, using numba.jit to achieve an acceptable performance. For large system sizes however, this method becomes quite slow, so I have been rewriting the routine into Fortran hoping to achieve a speed-up. However I have found that my Fortran version is much slower than the first version in Python, by a factor of 2-3.

I believe the bottleneck is a linear interpolation function that is called at each innermost loop. In the Python implementation I use numpy.interp, which seems to be pretty fast when combined with numba.jit. In Fortran I wrote my own interpolation function, which reads,

  complex*16 function interp(x_dat, y_dat, N_dat, x)
    implicit none
    integer, intent(in) :: N_dat
    real*8, dimension(N_dat), intent(in) :: x_dat
    complex*16, dimension(N_dat), intent(in) :: y_dat
    real*8, intent(in) :: x

    complex*16 :: y, y1, y2
    integer :: i, i1, i2, im

    if(x <= x_dat(1)) then
      y = y_dat(1)
    else if(x >= x_dat(N_dat)) then
      y = y_dat(N_dat)
    else
      im = MINLOC(DABS(x_dat - x), DIM=1)
      if(x_dat(im) >=x ) then
        i1 = im
        i2 = im - 1
      else
        i1 = im + 1
        i2 = im
      end if

      y1 = y_dat(i1)
      y2 = y_dat(i2)

      y = y1 + (x-x_dat(i1))*(y2 - y1)/(x_dat(i2) - x_dat(i1))
    end if
    interp = y
    return
  end function interp

Note that I need to interpolate complex data. If my diagnostics are correct this function is much slower than numpy.interp, which, since the interpolation has to be called in each loop, greatly reduces the speed of the whole program.

Does anyone know if there is a way to achieve Numpy like interpolation speeds in Fortran? Or if my interpolation function as shown above here is somehow horribly inefficient? I don't have much experience coding Fortran yet.

Thanks!

1) Have you profiled this so you know which line(s) are taking the most time? 2) How big is n_dat? 3) Is x_dat sorted? 4) Can you please, provide a complete program, without that it is much more difficult to say anything 5) Note `DABS` should not have been used since 1977, use `ABS` 6) Note complex*16 is not standard Fortran, and has never been standard Fortran - learn about kinds https://stackoverflow.com/questions/838310/fortran-90-kind-parameter to do it properly — Ian Bush, Feb 06 '22 at 12:57

veryreverie · Accepted Answer · 2022-02-06T15:42:33.067

At a guess (and see @IanBush's comments if you want to enable us to do better than guessing), it's the line

im = MINLOC(DABS(x_dat - x), DIM=1)

which is taking all your time, as this line is O(N) in the size of x_dat, where everything else is O(1).

If x_dat is linearly spaced then you can replace this line with

im = 1 + nint((N_dat-1)*(x-x_dat(1))/(x_dat(N_dat)-x_dat(1)))

or better yet, skip im entirely and calculate i1 and i2 as

i1 = 1 + floor((N_dat-1)*(x-x_dat(1))/(x_dat(N_dat)-x_dat(1)))
i2 = i1 + 1

If x_dat is not linearly spaced, but has other useful properties, then you want to use these properties to calculate im if possible.

My data is indeed linearly spaced, silly that I didn't think of this immediately. Thanks a bunch! — Jasper, Feb 06 '22 at 19:28

Achieving numpy like fast interpolation in Fortran

1 Answers1