在Fortran中实现类似于numpy的快速插值

Question

我有一个数值程序需要运行来解决某个包含几个嵌套循环的方程。最初我用Python编写了这个程序，并使用numba.jit来实现可接受的性能。然而对于大型系统，这种方法变得相当慢，因此我一直在重写程序，希望用Fortran实现加速。但是，我发现Fortran版本比Python第一版慢得多，约为2-3倍。

我相信瓶颈在于每个最内层循环中调用的线性插值函数。在Python实现中，我使用numpy.interp，在与numba.jit结合使用时似乎非常快。在Fortran中，我编写了自己的插值函数，其代码如下：

  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

请注意，我需要插值复杂数据。如果我的诊断正确，这个函数比 numpy.interp 慢得多，因为每次循环都必须调用插值函数，从而大大降低了整个程序的速度。

有没有人知道如何在Fortran中实现类似Numpy的插值速度？或者，如果我上面展示的插值函数效率极低，该怎么办呢？我对Fortran的编程经验不是很丰富。

谢谢！

- Jasper

---

1. 您是否对此进行了分析，以便知道哪些行花费的时间最长？
2. n_dat有多大？
3. x_dat是否已排序？
4. 请提供一个完整的程序，否则很难说出任何内容。
5. 注意自1977年以来不应使用“DABS”，请使用“ABS”。
6. 请注意，complex*16不是标准Fortran，也从未成为标准Fortran - 若要正确处理，请学习关于kind参数的知识 https://dev59.com/oHRA5IYBdhLWcg3wwwzD。

- Ian Bush

Answer 1

猜测（如果您想让我们的猜测更加准确，请参见@IanBush的评论），它指的是该行：

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

这一行代码是O(N), 其中N为x_dat的大小，耗费了所有的时间。如果x_dat是线性的，则可以将此行代码替换为

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

或者更好的方法是跳过im，并将i1 和 i2 计算为：

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

如果 x_dat 不是等间距的，但具有其他有用的属性，则应尽可能使用这些属性来计算 im。

- veryreverie

---

我的数据确实是线性间隔的，真傻我没有立刻想到这一点。非常感谢！ - Jasper