ループ階層とループ数は最小化するべし(粒子法による流体シミュ開発の中で)


粒子が循環するようにしてみた。
粒子の色の色相はその世代(生まれた時刻)を、粒子の色の明度はその粒子が感じている圧力の対数強度を表している。

今日はこの粒子間相互作用の計算方法のコツについてメモ。
まず空間を、相互作用のカットオフ距離と同じ長さの間隔をもつグリッドで分割し、グリッド単位で粒子をグループ分けする。ある粒子xと、他の粒子との相互作用は、粒子xと同じグリッド内に存在する全ての粒子と、その周囲グリッド(2次元だと縦横斜めの8グリッド)内の全ての粒子との相互作用を検討し、計算する。
と、ここまでの話は同じなのだが、以下のような分岐点がある

  • 各粒子に、自分が属しているグループの情報を持たせ、全粒子についてループをまわし、計算していく。
  • 各グリッドごとに、その中にある全ての粒子について、計算していく。

注目しているグリッドの周囲グリッドを決定する処理が、前者では各粒子ごとに必要となるのに対し、後者では各グリッドごとにしか発生しないので、後者の方が若干高速になると思うかも知れない。しかし実際の結果は逆となる。前者の方が、1.4倍くらい速く、しかも全粒子に対するでかいループが一つあるだけなので、並列化もしやすい。後者の場合はループが細切れに細分化されているので並列化も絶望的だ。計算量は少なくなるのだが、でかいループが一つだけあるというわけではなくなってしまうのである。こうなるとまず、どのループを並列化するべきか迷う。そして実際一番重そうなループの並列化をトライした結果(うちのマシンは2コアしかないが)、前者の方が総合的には2倍高速となった。
むしろ気を利かせた後者のスピードは倍近く遅くなったのである。
これは何故か。現象だけ見て考えると、前者の方はループ構造がより単純であるということだ。つまりループに階層がない。粒子単位で周囲グリッドを調べるから余計な参照やチェックは増えているが、ループ階層はなく、一つのでっかいループが回っているだけ、ということがいえると思う。
実際シングルスレッドでは前者が7秒、後者が10秒という程度の差なので大差はないが、並列化のしやすさまで考えると、その優劣は無視できない。このようにプロセッサ特性および並列化効率を考慮した最適なループの組み方があるのだ。(まあその道の一般論としては常識なのかも知れない…)
というわけで、シングルスレッドの場合でも、若干量の計算量の節約を考えるよりも、ループ構造を単純化することを真剣に考えた方がよい、並列化を考えるならなおさらだ、ということが言えそうだ。それはつまり、計算プログラムを組む場合は、ループ回数が大きくても、ループ自体の数と、ループ階層を、できるだけ少なく、単純化した方が、つまり、何万回も回すバカループが一つあるだけみたいな構造が、速いし、並列化もしやすい、ということだ。

粒子法 (計算力学レクチャーシリーズ)

粒子法 (計算力学レクチャーシリーズ)

粒子法シミュレーション―物理ベースCG入門

粒子法シミュレーション―物理ベースCG入門