Previous ToC Next

139. 半精度による相互作用計算(2017/12/30)

最近ぼちぼち色々な計算機で FP16 が使えるようになってきたので、 それで多体問題、少なくとも重力多体問題ができないかという雑談を 何人かでしていた時にでてきた方法について。

仮数長は 10ビットもあるので、GRAPE-1, 1A, 3, 5, 8 あたりに比べて長く、 相互作用計算にはいってからの精度には特に問題はありませんが、問題は

ということです。座標表現は、実効的に固定小数点にして、その上位と下位を 別々のワードに格納し、別々に引き算してから加算することで精度を確保する ことはできなくはありません。但し、2語使っても20ビットなので、3語くらい 使う必要がありそうです。

とはいえ、自分は原点近くにいて、遠くのほうは相対精度が落ちてよいとすれば、遠くは「桁を合わせない」という方法でなんかできる気がします。

積算は複数語のアキュムレータを使う方法が古くから知られています。

問題は途中でのオーバーフローです。距離のダイナミックレンジ自体は30ビッ ト程度はほしいので、途中で必要な距離の3乗の逆数を計算すると90ビットに なり、 FP16 で持つ5ビットの指数では表現できなくなります。

ここで、ツリー法を使い、質量の空間分布が2ないし3次元でほぼ一様な場合だ けをまず考えると、相互作用を計算する超粒子の質量は距離の 2 ないし3乗に 比例します。まず、色々上手くいきそうな 2次元の場合を考えると、 距離と質量の平方根が比例するので、距離を質量の平方根でわったものは ほぼ一定になります。従って、普通なら(dx 等は3次元ベクトルで内積があるとして)

   dx = xj-xi
   r2= dx*dx
   r=sqrt(r2)
   mr3inv=mj/(r*r*r)
   fi += mr3inv*dx
と書くところを

   q=1/sqrt(mj)
   dxhq = q*(xjh-xih)
   dxlq = q*(xjh-xil)
   dxq= dxhq+dxlq
   r2q2= dxq*dxq
   rq=sqrt(r2q2)
   mr3invqinv=1/(rq*rq*rq)
   fi += mr3invqinv*dxq
としてやると、 dxhq や dxlq はダイナミックレンジが小さく、基本的に1の オーダーになっているので、指数のオーバーフローやアンダーフローが ほぼ起きなくなることが期待できます。最初の、質量の逆数の平方根をとると ころは、ツリー構築を作った時に1度だけすればいいので、計算量の増加は無視できます。

3次元の場合には、一様分布の場合遠くの超粒子からの力が近くの粒子からの 力より大きくなりますが、これは距離に比例する程度なのでやはり上のような 修正を行うことでダイナミックレンジを大幅に小さくできます。

3次元で TreePM とか使うなら距離にはカットオフがあるのでこれで問題ない でしょう。開放境界の3次元問題の時には、結局、トータルの力が決まるあた りの特徴的な距離スケールがあるので、そこでオーバーフローが起きないよう にスケールを決めることができるはずです。
Previous ToC Next