我正在尝试优化这段特定 fortran 代码的执行时间。

详细信息：i_gc 是一个 (ngpts, 3) 数组，其中包含每个网格点的 (i,j,k) 索引。这是完整的 NX x NY x NZ 网格的子集。子集在采样密度变化的意义上是连接的，但点不会随机分散在整个域中。

本质上，我在网格上有一个非本地算子，它是 3 个一维算子的张量积。如果我有完整的 3d 网格上的数据，那么在没有内部循环中的所有 if 语句的情况下编写这段代码会很容易。另外，我想避免为每个网格点存储一组 {i},{j},{k} 索引，但也许这是要走的路？

我想知道是否有比这更好的方法？如果需要，我确实有一些自由来重新排列这些点。也许我可以以某种方式利用它？我想我可以尝试像希尔伯特曲线这样的东西，但我不确定这对我的情况有多好，因为我有一个形状奇特的规则网格子集。

编辑：数学上的操作是这样的

$Y_{ijk} = \sum^{N^x}_{q=1} M^x_{qi} X_{qjk} + \sum^{N^y}_{q=1} M^y_{qj} X_{iqk} + \sum^{N^z}_{q=1} M^z_{qk} X_{ijq}$

其中总和仅超过子集中的点 (i_gc)

subroutine apply_lf4(ngpts, i_gc, mx,my,mz,x,y)
  implicit none
  integer, intent(in) :: ngpts
  integer, intent(in), dimension(:,:) :: i_gc
  real(8), intent(in), dimension(:,:) :: mx,my,mz
  real(8), intent(in), dimension(ngpts) :: x
  real(8), intent(out), dimension(ngpts) :: y

  integer :: ig,jg,i1,i2,i3,j1,j2,j3

  y=0.0d0

  do ig=1, ngpts
     i1=i_gc(ig,1); i2=i_gc(ig,2); i3=i_gc(ig,3)
    inner: do jg=1, ngpts
        j1=i_gc(jg,1)
        if (j1==i1) then
           j2=i_gc(jg,2)
           if (j2==i2) then
              j3 = i_gc(jg,3)
              y(ig) = y(ig) + mz(j3,i3)*x(jg)
              cycle inner
           end if
           j3=i_gc(jg,3)
           if (j3==i3) then
              y(ig) = y(ig) + my(j2,i2)*x(jg)
              cycle inner
           end if
        else
           j2 = i_gc(jg,2)
           if (j2 == i2) then
              j3 = i_gc(jg,3)
              if (j3 == i3) then
                 y(ig) = y(ig) + mx(j1,i1)*x(jg)
                 cycle inner
              end if
           end if
        end if                        
     end do inner
  end do             

end subroutine apply_lf4