CPUを代表とする高速プロセッサでは、「演算パイプラインをフル稼動させるか」と「キャッシュのミスヒットをいかに少なくするか」を追求することが、実行性能の向上にとって必須の条件です。本ノートでは連立一次方程式のLU分解を例にとり、高速化を図った計算アルゴリズムについて説明します。FastPageDRAM×128bit×2Wayインターリーブ)1.LU分解を行なう連立一次方程式プログラム連立一次方程式の解を求める時に、殆どの時間を占めるのは行列のLU分解の部分です。一般的なプログラムは下記になります。Cacheのシステムでは52.4秒かかります。このルーチンでは浮動小数点の計算が666,166,500回行われますので、12.7MFLOPSの性能です。ピーク性能が1GFLOPSのマシンとしては極めて不満足な結果です。2.LU分解プログラムの問題点と高速化アルゴリズムこのプログラムを眺めてみると最内側のループが行列の行方向(列番号の変化する方向)に計算を進めています。式の部分がDAXPYと呼ばれるルーチンです。K=1、I=2の時を例にして説明します。DAXPYルーチンでアクセスされる行列は下記の順になっています。・FORTRANでは行列は下記の様に列方向(行番号が変化)にメモリー上に配置されています。従って上のDAXPYでの行列のアクセスはメモリーアドレスの離れた所をアクセスしています。キャッシュのないシステムではアドレスの離れた所をアクセスしても、隣接する所をアクセスしてもアクセス時間は変りませんが、キャッシュが付いていると全く異なります。CPUが一旦データを読み込むと、その前後の1ブロックがキャッシュ内に格納されます。Alpha21164ではこのブロックのサイズが64バイトですので、必要とするデータの前後の8個の倍精度実数がキャッシュの中に取り込まれます。ところが、このプログラムでは一度アクセスが発生すると次には遠く離れた部分をアクセスするのでキャッシュに取り込まれたデータのうち1個しか使用されません。性能を向上するためには、計算のアルゴリズムをアクセスするアドレスが連続する様に変更しなければなりません。3.キャッシュのミスミットを抑えたアルゴリズム1章のプログラムではある対角項A(K,K)を含んだ行に基づいて行単位に処理を行なっていますが、下のプログラムでは列単位に処理を行なっています。対角項が決定すると同一列の項をその対角項で除算しLマトリックスを生成します。新しく生成されたLを用いて消去計算を行、列ともK+1以上の全領域に対して行なっています。この時、計算は列に沿って行ないますので、キャッシュに取り込まれた内容が有効に使用されます。このプログラムの実行時間は25.1秒(26.5MFLOPS)です。------------------ENDDORETURNEND4.DDOTを用いたアルゴリズムへの変更1プログラム1から2への変更で、計算時間は大きく改善されましたが、まだ不満足です。DAXPYルーチンを用いた両プログラムで例えばA(4,5)はどの様な計算結果が格納されるかを検討します。A(4,3)は生成されたLこれらの式を総合するとA(4,5)=A(4,5)-A(4,1)*A(1,5)-A(4,2)*A(2,5)-A(4,3)*A(3,5)になります。この計算方法を採用したものが下記のプログラム3です。-----式?ENDDOA(I,J)=A(I,J)-TENDDOENDDORETURNEND上右図の塗りつぶした部分はK=4の時に計算する部分です。同部の計算には×印で示した領域のデータが使用され、Kが大きくなるに従って計算する部分は矢印の方向に進行します。式?を含むループはDDOTと呼ばれるルーチンに相当します。この式ではDAXPYに比べて式当りの行列参照が1回少なく、配列への書込が発生しません。従ってプログラム2よりも高速化が期待できるのですが、結果は予期に反して29.9秒(22.3MFLOPS)と遅くなっています。これは式中の行列のアクセスが一方は列方向に進行してキャッシュを有効利用できるのに対し、他方は行方向に進行してキャッシュのミスヒットが多発するために生じたのです。5.DDOTを用いたアルゴリズムへの変更2プログラム3での問題は、あるKで行の連続アクセスによってキャッシュミスヒットが連発するだけでなく、K+1でも同じ所を再びアクセスしてミスヒットが連発し、しかもそのミスヒット時に遅い外部メモリーから読み込んだデータの1/8しか使用されないことです。プログラム4ではNSTEP(ここでは25)に達するとそれまでに計算された領域の結果を用いて未計算の部分の計算を全て行なってしまい、以後その領域を参照していません。Kが26から50あるいはそれ以降についても同様にしています。KがNSTEPに達するまではプログラム3と同様に行います。この計算はMUL_TVEC_MAT、MUL_MAT_VECとしてサブルーチン化しています。KがNSTEPだけ進むと行列同士の積の結果を行列から減じるMUL_MAT_MATを行なっています。一例として、A(60,70)に施される計算は下記の様になります。式?Kが例えば25になれば、A(*,1)からA(*,25)まではL行列が形成されていますので、その値を用いて部分的に計算しておくことができます。これが式?の部分です。この計算が右図の塗りつぶした領域の全てについて行われますが、それは行列Aと行列Bの積の結果を塗りつぶした行列から減じることと同じです。Aでは行単位に、Bでは列単位に配列の内容が参照されますが、どちらも参照する要素数は25ですので、一度参照すると8KBの一次キャッシュの中に格納されてしまいます。サブルーチンの中では列方向のループよりも行方向のループを内側にしていますが、行列Aの1行を参照すると、Alphaのキャッシュブロックが64Byteであることからキッシュの中にその行の上下8行分が格納されますので、次の行を参照した時にはキャッシュの中からデータを読み出すことができるからです。以上の改善を施したプログラム4の実行時間は7.3秒(91.2MFLOPS)と非常に高速化することができました。以上まではキャッシュを有効に利用するアルゴリズムでしたが、以下では浮動小数点演算パイプラインを効率良く働かせるためのプログラムの改良を行ないます。その準備としてプログラム4の中のサブルーチンMUL_MAT_MATを修正し、NSTEP(=25)単位に計算を進める様にしています。このプログラムの計算時間は5.7秒(116MFLOPS)になりました。C(I,J)-TENDDOENDDORETURNENDプログラム6はプログラム5で設けたサブルーチンMUL_MT_MTの内部をパイプラインが並列に動作する様に修正しています。この結果、実行時間は3.2秒と短縮され、208MFLOPSの性能を実現しました。C(I,J)-T1-T2-T3-T4-T5ENDDOENDDORETURNEND7.むすびAlphaの様な高速コンピュータの性能を最大限に引き出すためには、キャッシュミスヒットの回数をできるだけ抑え、パイプラインの動作の並列化を実現するアルゴリズムの適用が重要であることを連立一次方程式のLU分解のルーチンを例にとり説明しました。本レポートではオーソドックスなプログラムから最適プログラムまで16倍強の性能向上を実現していますが、この考え方は多くの数値計算分野に適用できるものと確信しています。このプログラムは25の整数倍の行列に適用できるものとなっていますので、使用される時には注意が必要です。