KAGAKUGIJYUTUKEISAN99-1-BUTUSEIKEN2

「二次元古典イジング(Ising)模型」と呼ばれる模型を考える。この模型では、二次元の正方格子の各格子点に「1 (上向き)」、もしくは「-1 (下向き)」のどちらかの値を持つ「スピン」が存在し、それぞれの隣り合う二つのスピンがお互いに同じ向きを持つか、反対の向きを持つかによって、異なるエネルギーが与えられる。この模型は、現実の磁性体の相転移を記述する最も単純な模型の一つとして、非常に詳しく研究されてきた。この系の自由エネルギー(free energy)を(計算機の数値精度の範囲で)厳密に求めようとすると、もっとも素直な方法では、位相空間の大きさと同じ、O(2^N)の手間がかかる。(ここでNは系に含まれるスピンの個数である。 $L \times L$ の系を考える場合、計算の手間は O(2^L²)となる。) この計算の手間は、以下に述べる「転送行列 (transfer matrix)法」を用いることにより、 O(L² 2^L)にまで減らすことができることが知られている。

転送行列法では、次元が $2^L\times2^L$ の行列を考える。この行列は、二次元格子を一方向にスライスした時の各スライスの自由エネルギーへの寄与を表わしている。この転送行列をL回掛け合わせた後、トレースを取ることにより、系の分配関数(partition function)が、さらにその対数を取ることにより、自由エネルギーが求められる。

転送行列法のもう一つの利点は、半無限系( $L\times\infty$ )を扱うことができるという点にある。転送行列の無限回の積を考えると、自由エネルギーに寄与するのは、転送行列の最大固有値のみとなる。また、最大固有値と第二固有値の比から、系の相関長を求めることもできる。こうして、古典イジング模型の半無限系の自由エネルギーを求めるという問題は、「転送行列の固有値問題」へと帰着する。一般に、転送行列は $2^L\times2^L$ の正定値の実対称行列である。全ての行列要素の保存にはO(2^2L)の主記憶容量が必要であり、べき乗法、ランチョス法などの反復法を用いて最大固有値を求めるとしても、O(2^2L)の計算時間が必要となる。したがって、一台のプロセッサを用いた場合には、L=14が実質的な限界となる( $2^{2\times14}\times 8{\rm (byte)} = 2 {\rm (GB)}$ )。

しかし、系の持つ物理的特性(すなわち、相互作用が近距離の2体相互作用の和からなる)を利用することによって、転送行列はL個の「疎な」行列の積に書き直すことができる。それぞれの疎行列は、各行(各列)で対角成分と一個の非対角成分を除いては全て零となっている。この疎行列分解を用いると、反復法に必要な、ベクトルと行列の積の計算にかかる手間は $O(L\times2^L)$ にまで減少する。さらに、行列を記憶するために必要な主記憶容量も同程度にまで減少する。(零でない行列要素を「その場で」計算すれば、必要な主記憶容量をO(2^L)にまで減らすことも可能である。) この方法により、1プロセッサで最大 $L=26\sim28$ 程度の大きさの系まで扱うことができる。

２．２並列化

ここでは、疎行列の要素をあらかじめ主記憶上に記憶しておくのではなく、全てその場で計算する場合について考える。この場合、行列とベクトルの積の計算を行なうのに必要となる主記憶容量は、 $N_{\rm d}=2^L$ の長さのベクトル一本分のみである。(計算結果をもとのベクトルに記憶することにより一本で済む。) このベクトルを各プロセッサに均等に分割して配置する。ただし、プロセッサ数( $N_{\rm p}$ )は2のべき乗であると仮定する。反復法を行なうのに必要な、ベクトルのノルム計算、ベクトルの内積計算、ベクトルの定数倍等の操作についての並列化は自明であるので、以下、行列とベクトルの積の計算についてのみ考えることにする。

転送行列全体としては密であるので、最も単純に考えると、各プロセッサがそれぞれに割り当てられたベクトルについての計算を行なうためには、 $N_{\rm d}(N_{\rm p}-1)/N_{\rm p}$ の量のデータをネットワーク経由で他のプロセッサから取り寄せなければならないように思える。しかし、すでに述べた疎行列分解の方法をつかうと、以下に示す通り、データ転送量は $2N_{\rm d}(\log_2 N_{\rm p})/N_{\rm p}$ で済むことがわかる。

図1に、 $N_{\rm p}=8$ の場合の転送行列の疎行列分解を示す。もともとの転送行列は密であるが、四つの疎行列に分解した後では、それぞれ、対角部分行列とそれ以外にただ一つの零でない要素を持つ部分行列からなっていることがわかる。したがって、各疎行列の積のステップでは、各プロセッサは、それぞれペアを作って、ベクトルのトランスポーズ転送を行ない、各部分積を計算した後、逆トランスポーズ転送を行なうことになる。転送パターンは、FFTなどでおなじみのバタフライ演算と同じパターンであり、一プロセッサあたりの転送量は、合計 $2N_{\rm d}(\log_2 N_{\rm p})/N_{\rm p}$ となる。

さらに、「動的領域分割」の手法を用いることにより、転送量を(漸近的に)半分にすることができる。上記で述べた並列化では、各プロセッサの受け持つ部分ベクトルの割り当ては静的であった。計算の終わった部分ベクトルをもとのプロセッサに送り返す代わりに、各プロセッサへの割り当てそのものを動的に変化させることにより、二度目のトランスポーズ転送は不要になる。この方法では、最終的な辻褄を合わせるために一回の転送が余分な必要となるので、一プロセッサあたりの転送量は、合計で $N_{\rm d}(\log_2 N_{\rm p}+1)/N_{\rm p}$ となる。この場合の一プロセッサあたりの演算量と転送量との比は、

$\begin{displaymath}\frac{N_{\rm d}(\log_2 N_{\rm p}+1)/N_p}{(log_2 N_{\rm d}) N_{\rm d}/N_{\rm p}} \sim \frac{\log_2 N_p}{\log_2 N_{\rm d}} \end{displaymath}$

２．３ VPP500のMPI実行性能評価

厳密対角化法の並列化においては、行列のトランスポーズ転送が並列時間においてかなりの部分を占めている。この転送は、回数は少ないが一回の転送量は数十MB～数GBである。そのため、並列実効性能には、ネットワークの立ち上がり性能よりむしろ、スループット性能が重要となる。そこで、我々はまず、 VPP500におけるネットワークスループット実効性能を簡単なベンチマークプログラムにより評価した。MPIを用いた場合、プロセッサ間でのベクトル(配列)の交換を行なう方法としては、

まず気が付くのは、MPI_SENDRECVの転送速度が、片方向の転送の場合の半分程度であるという点である。すなわち、MPI_SENDRECVでは、ハードウェアの持つ全二重性はほとんど生かされておらず、片方向づつ順番に転送を行なうのと同程度の性能しか出ない。さらに、MPI_SENDRECV_REPLACEでは、ライブラリルーチン内で余分なメモリコピーが生じることによるオーバーヘッドを差し引いたとしても、非常に低い性能に留まっている。上記の4種類の転送方法の中で最も高速なのは、「(3) MPI_ISENDとMPI_RECVの組合せ」である(ただし、メッセージ長があまり長くない場合には、MPI_SENDRECV とほぼ同じ性能に留まる)。以下では、我々は、「(1) MPI_SENDRECV」と「(3) MPI_ISENDとMPI_RECVの組合せ」の二通りの転送方法を用いてプログラムを作成し、性能評価を行なった。

２．４並列性能評価

図3に、VPP500における16PEまでの並列化効率の実測値を示す。「MPI_ISENDとMPI_RECVの組合せ」を用いたプログラムでは、ピーク転送性能から見積もった並列化効率と比較しても、十分な並列化効率が得られている。一方、MPI_SENDRECVを用いたプログラムでは、前者と比較して1割以上の性能低下が見られる。