计算科学 - Fortran，制作更高效的双线性插值 - 吾爱随笔录

Fortran，制作更高效的双线性插值

计算科学正则插值

2021-11-30 22:28:04

在阅读了一些食谱之后，我正在尝试编写一个有效的双线性（2D）插值，作为 Matlab 的 fortran-mex，它在太阳图像处理的长算法中被广泛使用，因此是我的主要瓶颈之一。我的第一次尝试工作但很慢。请参阅下面的代码。任何可以加快代码速度的建议都是最受欢迎的。我正在使用 gcc 版本 4.3.6。我不能使用以后的版本，因为它是唯一与我的 Matlab 版本完全兼容的版本。任何更高版本都阻止使用一些嵌入式 matlab-mex 函数，这些函数是在 matlab 提示中连接/显示错误消息所必需的。我的 mex-setup 可以编译 Fortran 90 和 95（到目前为止...）

到目前为止，我已经使用了一个 C++ mex 等价物，并希望编写良好的 Fortran 代码会更快，我使用 C++ 作为 Fortran 可以或不能提供的加速的参考。

SUBROUTINE L2DINTERPOL3(IntIm,Image,x,y,NPts,M,N)

  INTEGER*8        :: NPts,M,N,jj
  INTEGER*8        :: x1(NPts),y1(NPts),x2(NPts),y2(NPts)
  DOUBLE PRECISION :: IntIm(NPts)
  DOUBLE PRECISION :: Image(M,N)
  DOUBLE PRECISION :: f1,f2
  DOUBLE PRECISION :: x(NPts),y(NPts)
  DOUBLE PRECISION :: wx1(NPts),wx2(NPts),wy1(NPts),wy2(NPts)

c Take the nearest integers around the input. 

   x1 = FLOOR(x)
   x2 = CEILING(x)
   y1 = FLOOR(y)
   y2 = CEILING(y)

   wx1 = (x2-x)/(x2-x1)
   wx2 = (x-x1)/(x2-x1)
   wy1 = (y2-y)/(y2-y1)
   wy2 = (y-y1)/(y2-y1)

c Whenever the input are already integers, weights go to infinity,
c so set each pair to 1 and 0.

   WHERE(x1 .EQ. x2)
      wx1 = 1.
      wx2 = 0.
   END WHERE

   WHERE (y1 .EQ. y2)
      wy1 = 1.
      wy2 = 0.
   END WHERE

c Main calculation, loop over each element of the input. 

  DO 10 jj=1,NPts


     f1 = wx1(jj)*Image(y1(jj),x1(jj)) + wx2(jj)*Image(y1(jj),x2(jj))
     f2 = wx1(jj)*Image(y2(jj),x1(jj)) + wx2(jj)*Image(y2(jj),x2(jj))

     IntIm(jj) = wy1(jj)*f1 + wy2(jj)*f2


 10   CONTINUE

  RETURN
  END

[编辑] 因为我只有常规网格，x2-x1 = 1 或 0，所以我制作了另一个不需要除法的版本。这摆脱了 where 构造。此外，似乎用“FORALL”循环替换 DO 循环在 2 个 CPU 上提高了约 15%。我的理解是 FORALL 在存在时使用多线程，并且每个 jj 的每个计算在 jj+k 或 jj-k 处独立于另一个，这使得 FORALL 语句是明智的。不是吗？

SUBROUTINE L2DINTERPOL(IntIm,Image,x,y,NPts,M,N)
  implicit none

  mwSize, PARAMETER             :: dp = kind(0.d0) ! Double precision
  mwSize                        :: NPts, M,N       ! Input
  REAL(dp),DIMENSION(Npts)      :: x,y             ! Input
  REAL(dp),DIMENSION(M,N)       :: Image           ! Input
  mwSize                        :: jj
  mwSize, DIMENSION(NPts)       :: x1,y1,x2,y2      
  REAL(dp),DIMENSION(Npts)      :: wx1,wx2,wy1,wy2
  REAL(dp),DIMENSION(Npts)      :: IntIm           ! Output


  x1 = FLOOR(x)
  x2 = CEILING(x)
  y1 = FLOOR(y)
  y2 = CEILING(y)

  wx1 = x2-x
  wy1 = y2-y

  wx2 = 1 - wx1
  wy2 = 1 - wy1

  FORALL( jj=1 : NPts)

     IntIm(jj) = wy1(jj)*(wx1(jj)*Image(y1(jj),x1(jj))+wx2(jj)*
 $        Image(y1(jj),x2(jj))) + wy2(jj)*(wx1(jj)*
     $        Image(y2(jj),x1(jj))+wx2(jj)*Image(y2(jj),x2(jj)))

  END FORALL


  RETURN
  END

2个回答

您的解决方案对我来说看起来非常理想（假设它工作正常，我还没有检查过）。我不认为你可以优化很多——切换语言、编译器、优化标志可能会给你 20% 的加速，但没有什么能改变这个函数是你的瓶颈的事实。如果您需要使您的代码更快，那么您将需要更改频繁调用此函数的算法。

您是否尝试对 $\tt x$ 调用 L2DINTERPOL3 之前的数组？我的意思是只有 $\tt x$ 大批。

我正在考虑的是基于以下两个数组彼此相邻的事实 $M$ - 增量长度（即 Fortran 数组模型中的主要列），

Image( * , x1(jj) ) and Image( * , x2(jj) )

由于 ${\tt x2(jj)} = {\tt x1(jj)} +1$ ，自从 ${\tt x1(jj)} = [x_{\rm jj}]$ 和 ${\tt x2(jj)} = [x_{\rm jj}] +1$ . 排序 $\tt x$ 在调用 L2DINTERPOL3 之前可能会有所帮助 $\tt x1(jj+1)$ 靠近 $\tt x1(jj)$ 当我们继续 $\tt jj = jj+1$ 在循环。因此，它可以优化对存储的访问 $\tt Image(:,:)$ .

这个想法简单描述如下：

   call SORT( NPts, x, iorder )

   do j = 1, NPts
       x(j) = x( iorder(j) )
       y(j) = y( iorder(j) )
   enddo

   call L2DINTERPOL3(IntIm,Image,x,y,NPts,M,N)

   do j = 1, NPts
       IntIm( iorder(j) ) = IntIm(j)
   enddo

最后一个循环要恢复的地方 $\tt IntIm(:)$ 基于 $\tt iorder(:)$ . 你能实现这个想法并让我知道它的效果吗？

其它你可能感兴趣的问题

上一篇是否有更快的方法来计算矩阵的几何级数？下一篇证明自适应有限元的收敛 - min res FEM？