Car-Parrinello法のためのリコンフィギャラブル三次元FFT専用プロセッサ

佐々木 徹, 別役 潔, 樋口 高年, 長嶋 雲兵

Return

1 はじめに

固体や表面の第一原理計算手法として、最も良く利用される手法に、結晶の周期性を考慮した平面波展開をもちいるCar-Parrinello法（CP法）[1]がある。CP法が登場する前の従来法、すなわち、ハミルトニアン行列を直接対角化する手法は、平面波基底の数をNとして、Nの３乗の計算負荷を要していた。CP法は固有値問題を繰り返し法による最小値問題に置き換えると同時にFFTを用いることで、固有状態の数をMとして、計算負荷をMN ln Nまで削減することに成功した。これにより、第一原理分子動力学(MD)への道が開かれた。
しかし、単体のPCやWSでは、メモリや計算時間の制約により、第一原理MDによって扱える原子数は100原子程度に限られる。したがって、最近のナノテクノロジー研究や材料開発現場においてしばしば必要とされる1,000原子を超える大規模計算への要望に応えられる状況にはない。
一方、科学技術計算分野において、汎用CPUに代わり、FPGAなどのリコンフィギャラブルデバイスが注目されてきており、信号処理や画像処理、流体解析[2]、X線結晶構造解析などにも使用され始めている。そのため、ソフトウェアによる高速化の研究[4 - 7]も盛んに行われており、リコンフィギャラブルデバイス上におけるFFTの高速化手法の研究[8]も行われている。
そこで、われわれは、CP法において計算負荷の高い箇所をFPGAへオフロードすることにより、大規模な第一原理計算を、個人で占有可能なハードウェアにおいて実現するための研究に着手した。ここに、個人で占有可能とは、低消費電力、高コストパフォーマンス、デスクサイドに収まるサイズを想定している。なお、本研究の一部は科学技術振興調整費総合研究「科学技術計算専用ロジック組込み型シミュレータに関する研究」において為されたものである。以下では、左記プロジェクトをEHPCプロジェクトと記す。

2 CP法

Figure 1に本稿が高速化の対象としているCP法の処理フローの概略を示す[2]。非局所擬ポテンシャルの演算を除くと、CP法の計算処理時間の大部分は、固有状態の三次元FFTと固有状態間のGram-Schmidt(GS)直交化によって占められている。


Figure 1. Flow diagram of self-consistent field calculation in the CP method

われわれはオブジェクト指向設計とC++言語を利用し、CP法プログラムを新規に開発した。簡単のため、CPの原論文と同様に擬ポテンシャルを局所的とした[1, 2]。
われわれのCPコードを、原子数を変えて逐次実行した際のベンチマーク結果をTable 1に示す。バルクシリコンを例題として、最急降下法の１ステップについて計測した。ここで計測した200原子以下の計算規模では、処理時間の90％弱は波動関数の3次元FFTに費やされる。したがって、CP法を専用ロジック化により高速化するには、三次元FFTの高速化が第一の課題となる。
また、一般に科学技術計算で行うFFTは多次元であるが故にデータサイズが極めて大きいという特徴を持っており、CP法の場合でも実用規模として、例えば1000原子程度の系を計算するには、2563タップ（浮動小数点倍精度複素数では256MB）の計算規模となる[9]。

Table 1. Percentage of 3D FFT in the CP calculation
Count of Atoms864216
Count of Eigenstates16128432
Count of Plane-Waves1237392517
Count of FFT Taps163323643
Wall-clock time for one SD iteration[sec]0.37732724
3D FFT Percentage60.684.888.0
Orthogonalization Percentage1.75.65.9

3 三次元FFT専用ロジック

3. 1 三次元FFT

Figure 2に三次元FFTの実行方法を示す。素直に計算するには、スキャンする方向を順次替えて一次元FFTを行えば良いが、メモリ上の連続並びと異なる方向にスキャンする場合には、隣接データ間のストライドが大きくなってしまうため、ソフトウェアで高速に実行するには、スキャン方向を替える毎に配列を転置するなどの工夫が必要となる[7]。


Figure 2. Pseudo Burst Transaction for 3D Array

本稿で提案する三次元FFTロジックは、三次元配列データを転置せずに高速にメモリアクセスを行うことができる。FFTのバタフライ演算を実行するButterflyコアはあまり工夫しても、それに見合うだけのメモリバンド幅が実現できなければあまり大きな意味を持たない。そこで我々は専用ロジックが解決すべき問題点をメモリアクセスに絞り、FFTの計算アルゴリズムについてはできるだけ平易に実装できるものを選び、Figure 3に示すCooley-Tukeyアルゴリズムを基にButterflyコアを作成した。解決すべき最大の問題としたメモリアクセスについては、Figure 4のように擬似的なバースト転送をそれぞれのスキャン方向について行うことによって、XYZ各軸に関して等方的なアクセスを行い、これをFigure 5のようにダブルバッファリングして、データ転送時間を隠蔽することにより解決した。これらをひとつのFPGAに実装すると、ロジックの全体構成はFigure 6のようになる。三次元FFTロジックについての詳細は文献[11]を参照されたい。


Figure 3. Butterfly Core


Figure 4. Scan directions of 3D FFT


Figure 5. Double buffering


Figure 6. FFT logic

3. 2 ゼロパディング領域の取り扱い

CP法に限らず逆変換時の精度保証のため、元のデータよりもかなり大きなデータ領域を用いてFFTを施すことも多い。本来データの定義されない領域はゼロパディング等により、計算結果に与える影響の少ないデータが埋められる。
CP法では実空間波動関数を波数空間（周波数空間）上に定義され、Figure 7のようにエネルギーの上限値により高周波成分がカットオフされるが、カットオフ球の直径の倍以上、ゼロパディングして解析領域を拡大する必要がある。これは、勾配ベクタ計算において逆フーリエ変換を施し、波動関数を実空間上の波動関数に変換した後、その点のポテンシャルを乗じて再度フーリエ変換して波数空間に戻すため、波数空間に戻したときの計算精度を確保するということが理由である。


Figure 7. FFT area for the CP methods

ゼロパディング領域では、メモリからバッファにフェッチした一本のスキャンラインがすべてゼロであることが多い。一本のスキャンラインのデータがすべてゼロであれば、FFTを施した結果もゼロなので、演算を行う必要がない。従って、ゼロパディング領域であることを検出してそのスキャンラインの演算を省略することにより演算時間の短縮を行った。

3. 3 CP法固有の処理

CP法では、単にIFFTやFFTを専用ロジックで実行するだけでは、システム上のメリットが小さいため、電荷密度の計算配ベクタの計算も専用ロジックで行うのが好ましい。
今回は暫定的にFPGAロジックの内部バスにEXTRAロジックとして接続してシステム評価を行っているが、本来Butterflyコアの演算器の接続を変形して実行させるか、あるいはリコンフィギャラブルなスキームを用いるべきであると考えている。これについては今後の検討課題である。

4 評価に使用したシステム

4. 1 システム構成

FPGAボードを用いたシステム全体はPCクラスタから構成され、CompactPCIバスを介して個々のPCにFPGAボードが接続されている。EHPCプロジェクトにおいて、ロジックの評価用にFPGAボードを作成した。本ボード上では、FFT以外にも、GSMAC-FEMに基づく流体解析専用ロジックなども稼動している[14]。　
コントロール・ノードとして組込み制御用のCPU（SH4）と、演算ノードとして300万ゲート相当のFPGAを4個搭載している。個々FPGAには512MBのSDRAMモジュールが接続できる[13]。

4. 2 API

「FFTエンジン」をタップ数等のパラメータを与えてオブジェクトとして生成し、FFTエンジンに対して波動関数のリストを渡せば電荷密度や勾配べクタが算出される。「FFTエンジン」は完全に抽象化されており、アプリケーションプログラムはFPGAボードの枚数や搭載されるFPGAデバイスの数などといったハードウェアの情報は一切知る必要なく、アプリケーションプログラムと接続することができる。

5 結果

5. 1 FFT単体性能

Figure 8に三次元FFTの単体性能の実測値を示す。ソフトウェアで実行する場合と異なり、タップ数によらず演算性能は500MFLOPSである。これは、搭載した5個の浮動小数点演算器が動作周波数100MHzで、ほぼ稼働率100%で稼動していることを示している。従って、EHPC-FPGAボードでの測定ではメモリバンド幅がボトルネックになっている、言い換えると懸案となっているメモリアクセス時間が完全に隠蔽された、と言うことができる。


Figure 8. Performance of 3D FFT on the FPGA Board

5. 2 FFTゼロパディング領域の演算省略

ゼロパディング領域を検出して、FFTを省略した場合の結果をFigure 9に示す。


Figure 9. Comparison of 3D FFT performance with or without zero Padding

なお、ここでのカットオフ半径はメッシュ数に換算すると

643TAP：13
1283TAP：19
2563TAP：36

である。カットオフ半径と解析空間の大きさの比率にも依存するが、各辺のタップ数をすべてN、有効範囲をMとすれば
演算削減率 ≒ ((M/N)2+(M/N)+ 1 )/ 3
と見積もられるので、演算量は1/2程度には削減される。

5. 3 汎用CPU との比較

FPGAボードの性能を評価するため，世界的に広く使われており，性能にも定評もあるFFTWをＣＰコードに実装した. 以下に単独計算，電荷密度計算，勾配計算の比較を示す. FFTWの結果はIntel Xeon 2.4GHzで測定した。FPGAボードでの計算ではゼロパディング削除機能を使用しており，FFTのデータパディング領域はそれぞれ27**3, 39**3, 73**3としている。 SCF計算部の処理時間は式（i）で評価できる。Si原子512個の系の計算における性能予測をおこなった. Si原子512個(1024バンド) ，カットオフエネルギー 16Haの系では，FFTWの場合には式(ii-a)，FPGAチップ1つの場合には(ii-b)のようにそれぞれ評価される。 Table 2に１チップ、１ボード（FPGA4個搭載）の場合の処理時間を示す。256×256×256の場合、FFTWよりも個別の処理で、ボードあたり10倍程度高速である。
(SCF計算) = （直交化計算＋その他） + (電荷密度計算)×3 + (勾配計算)×3 … (i)
2512 + 6.043*1024*3 + 9.900*1024*3 = 51489 sec (14h18m)…(ii-a)
2512 + 2.002*1024*3 + 3.836*1024*3 = 20446 sec (5h41m) …(ii-b)

Table 2. 3D FFT performance of general purpose CPU and FPGA chip in msec.
Number of Taps64×64×64128×128×128256×256×256
FFTW on Xeon @2.4GHzFFT only242914542
Charge Density605466043
Gradient Vector888489900
3D FFT Logic on FPGA ChipFFT only252051834
Charge Density282282002
Gradient Vector524323836
3D FFT Logic on FPGA BoardFFT only651459
Charge Density757500
Gradient Vector13108959

次にシステム性能として、Figure 1 中のSteepest Descent(SD)/Conjugate Gradient(CG)ループ一回の速度を見積ってみると、1チップによるスピードアップは5.34，非並列化部分は0.123であるので，アムダールの式で評価される加速率はFigure 10のようになる。したがって、FPGAボード2枚、すなわち、FPGAチップ8個を用いた場合に10倍の高速化が可能である．


Figure 10. Speed up of the CP code by the FPGA boards

6 考察

6. 1 演算性能

三次元FFTの持つメモリアクセスの特殊性を考慮し、演算器だけでなく、メモリシステムを三次元FFTに特化し、さらにゼロパディング領域の演算を省略することにより、Xeon@2.4GHzを搭載したPC上のソフトウェアで実行した場合と比較してボードあたり10倍程度の性能が達成できた。これは、専用ロジックを開発する場合に、単に実装する演算器や演算アルゴリズムだけでなく、メモリコントローラを含めて考えると有効な場合があることを示している、と考えられる。

6. 2 低消費電力

EHPC-FPGAボード全体の消費電力が約15WであるからPCの消費電力を100W程度に見積もると1/7程度の電力である。消費電力まで考慮すれば、PCに対して数10倍のコストパフォーマンス比となる。

6. 3 ダイナミック・リコンフィギャラブル・システム

CP法では計算規模が1000原子を超えたあたりから三次元FFTだけでなく、ベクトル直交化も大きな負荷となる[9]。例えばベクタ直交化としてGS直交化を用いた場合、Nを原子数として計算規模のオーダーは
3D-FFT：O(N2logN)
GS直交化: O(N3)
である。このふたつの処理は同時に実行されることが無いため、同一のFPGAデバイス上で両者を行うことができればHWの使用効率という点で有利である。従って、3D-FFTとGS直交化の専用ロジックを開発し、それをリコンフィギャラブルシステム上で、適時ロジックをデバイスにロードして実行するという方法が有力であると考えている。そのため、ボード上のコントローラCPU（SH4）のメモリ上にロジック・レジストリというオブジェクトを生成し、複数のロジックを予めレジストリに登録しておけば任意の箇所でロジックを動的に切り替えられる機構を設けている。
この機構を活用することによって、三次元FFTとGS直交化の両者を効率よく実行できるシステムを構築していくことも今後の課題である。

7 まとめ

FFT固有のバタフライ演算に伴うメモリアクセスを工夫し、さらにゼロパディング領域に対する演算の省略を行うことにより、低消費電力の高速三次元FFT専用プロセッサを開発することができた。
今後はさらにリコンフィギャラブルという特徴を活かすことにより、より高性能なCP法専用計算機の開発を目指す。

本研究の一部は科学技術振興調整費 総合研究「科学技術計算専用ロジック組込み型プラットフォーム・アーキテクチャに関する研究」（代表　村上和彰　九州大学教授）によるものである。

参考文献

[ 1] Car, R. and Parrinello, M., Phys. Rev. Lett., 55, 2471-2474 (1985).
 [ 2] 陽電子の第一原理計算では局所擬ポテンシャルで十分である。例えば、
Kawasuso, A., Chiba, T., and Higuchi, T., Phys. Rev. B, 71, 193204 (2005).
Kawasuso, A., Yoshikawa, M., Itoh, H., Chiba, T., Higuchi, T., Betsuyaku, K., Redmann, F., and Krause-Rehberg, R., Phys. Rev. B (2005), in press.
 [ 3] Yokokawa, Itakura, Uno, Ishihara and Kaneda, 16.4TFlops Direct Numerical Simulation of Turbulence by Fourier Spectral Method on Earth Simulator, SC2002 (2002).
 [ 4] http://www.fftw.org/
[ 5] http://www.ffte.jp/
[ 6] http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/index.html
[ 7] 高橋大介, 朴泰祐, 佐藤三久, Short Vector SIMD命令を用いた並列FFTの実現と評価, SACSIS2004, 277-286 (2004).
 [ 8] 黒滝俊輔, 鈴木紀章, 安生健一朗, 本村真人, 若林一敏, 天野英晴, FFTのDRPによるアクセラレーション手法, FPGA/PLD Conference and Exbition (2004).
 [ 9] 樋口高年,谷村直樹,大谷泰昭,佐々木徹,長嶋雲兵, 情報処理学会シンポジウムシリーズ ハイパーフォーマンスコンピューティングと計算科学シンポジウム(HPCS2003), 2003(4), 133 (2003).
 [10] 村上和彰, 稲垣祐一郎, 上原正光, 大谷泰昭, 小原繁, 小関史朗, 佐々木徹, 棚橋隆彦, 中馬寛, 塚田捷, 長嶋雲兵, 中野達也, 科学技術計算専用ロジック組込み型プラットフォームアーキテクチャの開発−プロジェクト全体像−, SWoPP2000 (2000).
 [11] 佐々木徹, 溝口大介, 長嶋雲兵, Car-Parrinello計算向け三次元FFTロジックの開発, 情報処理学会論文誌, 45, 313-320 (2004).
 [12] Payne, M. C., Teter, M. P., Allan, D. C., Arias, T. A. and Joannopoulos, J. D., Rev. Mod. Phys., 64, 1045 (1992).
 [13] 溝口大介, 荒木健悟, 佐々木徹, 青木すみえ, 棚橋隆彦, 科学技術計算向けFPGA基板の設計と評価, 情報処理学会研究報告ARC, NO.153-02 (2003).
 [14] 青木すみえ, 溝口大介, 荒木健悟, 石橋政一, 佐々木徹, 棚橋隆彦, FPGAによるGSMAC-FEM専用計算機の実装と評価, 日本流体力学会, 第17回数値流体力学シンポジウム, C11-5 (2003).

Return