Fortranの高性能ルックアップテーブル

Question

シミュレーションコードの一部には、特定の密度と温度で不透明度を見つける必要があります。このための分析的関係はありません。標準的な方法は、opacity(i,j)がdensity(i)とtemperature(j)の不透明度に対応する2D配列を使用し、正確な不透明度を見つけるために双線形補間を実行することです。

私のグループコードのボトルネックは、現在のところ、補間ルーチンが異なる密度と温度で約1億回と呼ばれ、ランタイムの約20％を占めています。現在のコードを以下に示します。改善のために使用できるトリックはありますか？私は私が正しく理解していれば、あなたはopc内の値を見つけるために、後で使用するインデックスを見つけるために、rhoTの全体をループをやっているオプションを使用して、インテル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の列がソートされている場合、バイナリ検索を書き込むほうが速いかもしれません（オーバーヘッドがあるため、テストする必要があります）。

また、私はridとtidを割り当ててもよくわかりません（r <= rho < rpの場合にのみridを指定するのは理にかなっているようです）。私はrhoTがどのように構築されているかについて何か不足しているかもしれません。

---

試すかもしれないトリック：rho*Mは（精度の3桁を取得するために1000で乗算）Mは、例えば2又は10の電源であってよい整数に変換します。整数に丸められた値は、要素が正しい（または最も近い）値のridである配列へのインデックスになります。正確なridが得られない場合でも、確認する範囲はずっと小さいかもしれません。スケールが線形でない場合は、rhoを最初に変換することができます。

---

別の可能なトリック：前ridとtidインデックスを格納します。呼び出しがかなり進化し続ける場合、新しいインデックスは以前のものに近い可能性があります。しかし、コードがある時点でパラレル化されなければならない場合は、呼び出し間で逐次的な依存関係が導入されるため、良い考えではありません。