Fortran中的高性能查找表

Question

我的部分模拟代码要求我找到一个给定的密度和温度的不透明度。这方面没有解析关系；标准的方法是使用二维数组，其中opacity(i,j)将对应于density(i)和temperature(j)的不透明度，并运行双线性插值以求精确的不透明度。

这是我们组当前代码中的瓶颈--每个时间步骤，在不同的密度和温度下，插值例程被调用了大约1亿次，占运行时的20%。下面显示了当前的代码--我可以使用什么技巧来改进它吗？我使用Intel Fortran 16和选项-O3 -xavx -mcmodel=medium

function smoothopc(den, temp, ig, opc, rhoT)

    implicit none
    real(kind=8), intent(in) :: den,temp 
    integer,intent(in) :: ig
    real(kind=8), intent(in), dimension(1:50, 1:50, 1:52) :: opc
    real(kind=8), intent(in), dimension(1:50, 1:2) :: rhoT

    real(kind=8) :: rho, te, smoothopc, r1, r2, t1, t2, &
            interpolation, denominator, a, b, c, d, opc11, &
            opc12, opc21, opc22, t, tp, r, rp

    integer :: rid,tid,i

    rho = den * 1d - 3          !g/cc
    te = temp / (1.6d - 19)     !eV
    tid = -1
    rid = -1
    do i = 1, 49
        r = rhoT(i, 1)
        rp = rhoT(i + 1, 1)
        t = rhoT(i, 2)
        tp = rhoT(i + 1, 2)
        if (rho .ge. r) then
            rid = i
        endif

        if (te .ge. t) then
            tid = i
        endif
    enddo

    r1 = rhoT(rid, 1)
    r2 = rhoT(rid + 1, 1)
    t1 = rhoT(tid, 2)
    t2 = rhoT(tid + 1, 2)
    opc11 = opc(rid, tid, ig + 4)
    opc12 = opc(rid, tid + 1, ig + 4)
    opc21 = opc(rid + 1, tid, ig + 4)
    opc22 = opc(rid + 1, tid + 1, ig + 4)

    denominator = (r2 - r1) * (t2 - t1)
    a = r2 - rho
    b = rho - r1
    c = t2 - te
    d = te - t1

    interpolation = a * (c * opc11 + d * opc12) + b * &
            (c * opc21 + d * opc22)

    smoothopc = interpolation / denominator

    return

end function smoothopc 

Answer 1

如果我正确理解，您将对整个rhoT进行循环，以查找随后用于在opc中查找值的索引。

如果对rhoT的列进行了排序，则编写二进制搜索可能会更快(存在开销，因此必须进行测试)。

另外，我不太理解您分配rid和tid的条件(只在r <= rho < rp的情况下分配rid似乎是合乎逻辑的)。关于rhoT是如何构建的，我可能遗漏了一些东西。

您可能尝试的一个技巧:将rho*M转换为整数，其中M可能是2或10的幂(乘以1000获得3位精度)。舍入整数的值将是数组的索引，数组的元素是rid的正确(或最近的)值。即使您没有得到正确的rid，您也可能需要检查的范围要小得多。如果刻度不是线性的，您可以先转换rho。

另一个可能的技巧是:存储以前的rid和tid索引。如果投资者的呼声是持续不断的，那么新的指数很可能接近之前的指数。但是，如果您的代码必须在某个时候被并行化，这不是一个好主意，因为这会在调用之间引入顺序依赖关系。