Alpha CPUの性能を引き出す計算アルゴリズム(1997/6/2)

← 戻る

Alpha CPUを代表とする高速プロセッサでは、「演算パイプラインをフル稼動させるか」と「キャッシュのミスヒットをいかに少なくするか」を追求することが、実行性能の向上にとって必須の条件です。本ノートでは連立一次方程式のLU分解を例にとり、高速化を図った計算アルゴリズムについて説明します。
使用しているシステムはAlpha21164-500MHz、一次キャッシュ8KB、二次キャッシュ96KB(以上オンチップ)、三次キャッシュ2MB、メモリ(70ns FastPageDRAM×128bit×2Wayインターリーブ)

１．LU分解を行なう連立一次方程式プログラム

連立一次方程式の解を求める時に、殆どの時間を占めるのは行列のLU分解の部分です。一般的なプログラムは下記になります。

＜プログラム１＞
SUBROUTINE LU_FACTORIZE(A,NDIM)
DOUBLE PRECISION A(NDIM,NDIM),T,U
DO K=1,NDIM-1
　U=A(K,K)
　DO I=K+1,NDIM
　　T=A(I,K)/U
　　A(I,K)=T　－－－－－－－－－－－－－－－－－Lマトリックスの生成
　　DO J=K+1,NDIM　－－－－－－－－－－－－－－
　　　A(I,J)=A(I,J)-T*A(K,J)　----式?　　　　｜消去計算
　　ENDDO　　－－－－－－－－－－－－－－－－－－
　ENDDO
ENDDO
RETURN
END

行列を1000元にして、このルーチンをDEC FORTRANでコンパイルして実行すると、Alpha500MHz + 2MB Cacheのシステムでは52.4秒かかります。このルーチンでは浮動小数点の計算が666,166,500回行われますので、12.7MFLOPSの性能です。ピーク性能が1GFLOPSのマシンとしては極めて不満足な結果です。

2．LU分解プログラムの問題点と高速化アルゴリズム

このプログラムを眺めてみると最内側のループが行列の行方向（列番号の変化する方向）に計算を進めています。式の部分がDAXPYと呼ばれるルーチンです。K=1、I=2の時を例にして説明します。DAXPYルーチンでアクセスされる行列は下記の順になっています。

A(K,J): A(1,1),A(1,2),A(1,3) ･ ･ ･ ･ ･
A(I,J): A(2,1),A(2,2),A(2,3) ･ ･ ･ ･ ･

FORTRANでは行列は下記の様に列方向（行番号が変化）にメモリー上に配置されています。
A(1,1),A(2,1),A(3,1), ･ ･ ･ ･,A(NDIM,1),A(1,2),A(2,2),A(3,2), ･ ･ ･ ･ ･

従って上のDAXPYでの行列のアクセスはメモリーアドレスの離れた所をアクセスしています。
キャッシュのないシステムではアドレスの離れた所をアクセスしても、隣接する所をアクセスしてもアクセス時間は変りませんが、キャッシュが付いていると全く異なります。CPUが一旦データを読み込むと、その前後の１ブロックがキャッシュ内に格納されます。Alpha21164ではこのブロックのサイズが64バイトですので、必要とするデータの前後の８個の倍精度実数がキャッシュの中に取り込まれます。ところが、このプログラムでは一度アクセスが発生すると次には遠く離れた部分をアクセスするのでキャッシュに取り込まれたデータのうち１個しか使用されません。性能を向上するためには、計算のアルゴリズムをアクセスするアドレスが連続する様に変更しなければなりません。

3．キャッシュのミスミットを抑えたアルゴリズム

1章のプログラムではある対角項A(K,K)を含んだ行に基づいて行単位に処理を行なっていますが、下のプログラムでは列単位に処理を行なっています。まず対角項が決定すると同一列の項をその対角項で除算しLマトリックスを生成します。次いで新しく生成されたLを用いて消去計算を行、列ともK+1以上の全領域に対して行なっています。この時、計算は列に沿って行ないますので、キャッシュに取り込まれた内容が有効に使用されます。このプログラムの実行時間は25.1秒(26.5MFLOPS)です。

＜プログラム２＞
SUBROUTINE LU_FACTORIZE(A,NDIM)
DOUBLE PRECISION A(NDIM,NDIM),T,U
DO K=1,NDIM-1
U=A(K,K)　　　－－－－－－－－－－－－－－－
DO I=K+1,NDIM　　　　　　　　　　　　　　｜Lマトリックスの生成
A(I,K)=A(I,K)/U　　　　　　　　　　　　｜
ENDDO　　－－－－－－－－－－－－－－－－－
DO J=K+1,NDIM　－－－－－－－－－－－－－－
DO I=K+1,NDIM　　　　　　　　　　　　　｜
A(I,J)=A(I,J)-A(I,K)*A(K,J)　　　　　　　｜消去計算
ENDDO　　　　　　　　　　　　　　　　　　｜
ENDDO　　－－－－－－－－－－－－－－－－－－
ENDDO
RETURN
END

４．DDOTを用いたアルゴリズムへの変更１

プログラム１から２への変更で、計算時間は大きく改善されましたが、まだ不満足です。DAXPYルーチンを用いた両プログラムで例えばA(4,5)はどの様な計算結果が格納されるかを検討します。
　　　　K=1の時　A(4,5)=A(4,5)-A(4,1)*A(1,5) A(4,1)は生成されたL
　　　　K=2の時　A(4,5)=A(4,5)-A(4,2)*A(2,5) A(4,2)は生成されたL
　　　　K=3の時　A(4,5)=A(4,5)-A(4,3)*A(3,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)
になります。この計算方法を採用したものが下記のプログラム３です。

＜プログラム３＞
　　　　SUBROUTINE LU_FACTORIZE(A,NDIM)
　　　　DOUBLE PRECISION A(NDIM,NDIM),T,U
　　　　DO K=1,NDIM-1
　　　　　U=A(K,K)
　　　　　DO I=K+1,NDIM
　　　　　　A(I,K)=A(I,K)/U
　　　　　ENDDO
　　　　　I= K+1
　　　　　DO J=K+1,NDIM
　　　　　　T=0.0
　　　　　　DO L=1,K
　　　　　　　T= T + A(I,L)*A(L,J)
　　　　　　ENDDO
　　　　　　A(I,J)=A(I,J)-T
　　　　　ENDDO
　　　　　J= K+1
　　　　　DO I=K+2,NDIM
　　　　　　T=0.0
　　　　　　DO L=1,K
　　　　　　　T= T + A(I,L)*A(L,J)　-----式?
　　　　　　ENDDO
　　　　　　A(I,J)=A(I,J)-T
　　　　　ENDDO
　　　　ENDDO
　　　　RETURN
　　　　END

上右図の塗りつぶした部分はK=4の時に計算する部分です。同部の計算には×印で示した領域のデータが使用され、Kが大きくなるに従って計算する部分は矢印の方向に進行します。式?を含むループはDDOTと呼ばれるルーチンに相当します。この式ではDAXPYに比べて式当りの行列参照が１回少なく、配列への書込が発生しません。従ってプログラム２よりも高速化が期待できるのですが、結果は予期に反して29.9秒(22.3MFLOPS)と遅くなっています。これは式中の行列のアクセスが一方は列方向に進行してキャッシュを有効利用できるのに対し、他方は行方向に進行してキャッシュのミスヒットが多発するために生じたのです。

５．DDOTを用いたアルゴリズムへの変更２

プログラム３での問題は、あるKで行の連続アクセスによってキャッシュミスヒットが連発するだけでなく、K+1でも同じ所を再びアクセスしてミスヒットが連発し、しかもそのミスヒット時に遅い外部メモリーから読み込んだデータの１／８しか使用されないことです。プログラム４ではNSTEP（ここでは25）に達するとそれまでに計算された領域の結果を用いて未計算の部分の計算を全て行なってしまい、以後その領域を参照していません。Kが26から50あるいはそれ以降についても同様にしています。
KがNSTEPに達するまではプログラム３と同様に行います。この計算はMUL_TVEC_MAT、MUL_MAT_VECとしてサブルーチン化しています。

KがNSTEPだけ進むと行列同士の積の結果を行列から減じるMUL_MAT_MATを行なっています。一例として、A(60,70)に施される計算は下記の様になります。

A(60,70)=A(60,70)-A(60,1)*A(1,70)-A(60,2)*A(2,70) ･ ･ ･-A(60,25)*A(25,70) 式?

A(60,70)=A(60,70)-A(60,26)*A(26,70)-A(60,27)*A(27,70) ･ ･ ･-A(60,50)*A(50,70) 式?

A(60,70)=A(60,70)-A(60,51)*A(51,70) ･ ･ ･ -A(60,59)*A(59,70) 式?

Kが例えば25になれば、A(*,1)からA(*,25)まではL行列が形成されていますので、その値を用いて部分的に計算しておくことができます。これが式?の部分です。この計算が右図の塗りつぶした領域の全てについて行われますが、それは行列Ａと行列Ｂの積の結果を塗りつぶした行列から減じることと同じです。Ａでは行単位に、Ｂでは列単位に配列の内容が参照されますが、どちらも参照する要素数は25ですので、一度参照すると8KBの一次キャッシュの中に格納されてしまいます。サブルーチンの中では列方向のループよりも行方向のループを内側にしていますが、行列Ａの１行を参照すると、Alphaのキャッシュブロックが64Byteであることからキッシュの中にその行の上下８行分が格納されますので、次の行を参照した時にはキャッシュの中からデータを読み出すことができるからです。

以上の改善を施したプログラム４の実行時間は7.3秒(91.2MFLOPS)と非常に高速化することができました。

＜プログラム４＞
　　　　SUBROUTINE LU_FACTORIZE(A,NDIM)
　　　　DOUBLE PRECISION A(NDIM,NDIM),T,U
　　　　NSTEP=25
　　　　KX=1
　　　　DO K=1,NDIM-1
　　　　　U=A(K,K)
　　　　　DO I=K+1,NDIM
　　　　　　A(I,K)=A(I,K)/U
　　　　　ENDDO
　　　　　KK=K-(K/NSTEP)*NSTEP
　　　　　KJ= K-KX+1
　　　　　IF(KK.NE.0) THEN
　　　　　　CALL MUL_TVEC_MAT(A(K+1,KX),A(KX,K+1),A(K+1,K+1),NDIM,KJ,KX)
　　　　　　IF(K.NE.NDIM-1) THEN
　　　　　　　CALL MUL_MAT_VEC(A(K+2,KX),A(KX,K+1),A(K+2,K+1),NDIM,KJ,KX)
　　　　　　ENDIF
　　　　　ELSE
　　　　　　CALL MUL_MAT_MAT(A(K+1,KX),A(KX,K+1),A(K+1,K+1),NDIM,NSTEP,KX)
　　　　　　KX=KX+NSTEP
　　　　　ENDIF
　　　　ENDDO
　　　　RETURN
　　　　END

　　　　SUBROUTINE MUL_TVEC_MAT(A,B,C,NDIM,K,KX)
　　　　DOUBLE PRECISION A(NDIM,*), B(NDIM,*), C(NDIM,*), T
　　　　JEND= NDIM-K-(KX-1)
　　　　DO J=1,JEND
　　　　　T=0.0
　　　　DO L=1,K
　　　　　T= T + A(1,L)*B(L,J)
　　　　　ENDDO
　　　　　C(1,J)=C(1,J)-T
　　　　ENDDO
　　　　RETURN
　　　　END

　　　　SUBROUTINE MUL_MAT_VEC(A,B,C,NDIM,K,KX)
　　　　DOUBLE PRECISION A(NDIM,*),B(NDIM),C(NDIM),T
　　　　IEND= NDIM-K-1-(KX-1)
　　　　DO I=1,IEND
　　　　　T=0.0
　　　　　DO L=1,K
　　　　　　T= T + A(I,L)*B(L)
　　　　　ENDDO
　　　　　C(I)=C(I)-T
　　　　ENDDO
　　　　RETURN
　　　　END
　　　　SUBROUTINE MUL_MAT_MAT(A,B,C,NDIM,NSTEP,KX)
　　　　DOUBLE PRECISION A(NDIM,*),B(NDIM,*),C(NDIM,*),T
　　　　IEND= NDIM-NSTEP-(KX-1)
　　　　DO J=1,IEND
　　　　　DO I=1,IEND
　　　　　　T=0.0
　　　　　　DO L=1,NSTEP
　　　　　　　T= T+A(I,L)*B(L,J)
　　　　　　ENDDO
　　　　　　C(I,J)= C(I,J)-T
　　　　　ENDDO
　　　　ENDDO
　　　　RETURN
　　　　END

６．更なる高速化への挑戦

以上まではキャッシュを有効に利用するアルゴリズムでしたが、以下では浮動小数点演算パイプラインを効率良く働かせるためのプログラムの改良を行ないます。
その準備としてプログラム４の中のサブルーチンMUL_MAT_MATを修正し、NSTEP(=25)単位に計算を進める様にしています。このプログラムの計算時間は5.7秒(116MFLOPS)になりました。

＜プログラム５＞
　　　　SUBROUTINE MUL_MAT_MAT(A,B,C,NDIM,NSTEP,KX)
　　　　DOUBLE PRECISION A(NDIM,*),B(NDIM,*),C(NDIM,*),T
　　　　IEND= NDIM-NSTEP-(KX-1)
　　　　DO J=1,IEND,NSTEP
　　　　　DO I=1,IEND,NSTEP
　　　　　　CALL MUL_MT_MT(A(I,1),B(1,J),C(I,J),NDIM,NSTEP)
　　　　　ENDDO
　　　　ENDDO
　　　　RETURN
　　　　END
　　　　SUBROUTINE MUL_MT_MT(A,B,C,NDIM,NSTEP)
　　　　DOUBLE PRECISION A(NDIM,*),B(NDIM,*),C(NDIM,*),T
　　　　DO J=1,NSTEP
　　　　　DO I=1,NSTEP
　　　　　　T=0.0
　　　　　　DO L=1,NSTEP
　　　　　　　T= T+A(I,L)*B(L,J)
　　　　　　ENDDO
　　　　　　C(I,J)= C(I,J)-T
　　　　　ENDDO
　　　　ENDDO
　　　　RETURN
　　　　END

プログラム６はプログラム５で設けたサブルーチンMUL_MT_MTの内部をパイプラインが並列に動作する様に修正しています。この結果、実行時間は3.2秒と短縮され、208MFLOPSの性能を実現しました。

＜プログラム６＞
　　　　SUBROUTINE MUL_MT_MT(A,B,C,NDIM,NSTEP)
　　　　DOUBLE PRECISION A(NDIM,*),B(NDIM,*),C(NDIM,*),T1,T2,T3,T4,T5
　　　　DO J=1,NSTEP
　　　　　DO I=1,NSTEP
　　　　　　T1= A(I,1)*B(1,J) +A(I,2)*B(2,J) +A(I,3)*B(3,J)
　　　@　　　　　+A(I,4)*B(4,J) +A(I,5)*B(5,J)
　　　　　　T2= A(I,6)*B(6,J) +A(I,7)*B(7,J) +A(I,8)*B(8,J)
　　　@　　　　　+A(I,9)*B(9,J) +A(I,10)*B(10,J)
　　　　　　T3= A(I,11)*B(11,J)+A(I,12)*B(12,J)+A(I,13)*B(13,J)
　　　@　　　　　+A(I,14)*B(14,J)+A(I,15)*B(15,J)
　　　　　　T4= A(I,16)*B(16,J)+A(I,17)*B(17,J)+A(I,18)*B(18,J)
　　　@　　　　　+A(I,19)*B(19,J)+A(I,20)*B(20,J)
　　　　　　T5= A(I,21)*B(21,J)+A(I,22)*B(22,J)+A(I,23)*B(23,J)
　　　@　　　　　+A(I,24)*B(24,J)+A(I,25)*B(25,J)
　　　　　　C(I,J)= C(I,J)-T1-T2-T3-T4-T5
　　　　　ENDDO
　　　　ENDDO
　　　　RETURN
　　　　END

７．むすび

Alphaの様な高速コンピュータの性能を最大限に引き出すためには、キャッシュミスヒットの回数をできるだけ抑え、パイプラインの動作の並列化を実現するアルゴリズムの適用が重要であることを連立一次方程式のLU分解のルーチンを例にとり説明しました。本レポートではオーソドックスなプログラムから最適プログラムまで16倍強の性能向上を実現していますが、この考え方は多くの数値計算分野に適用できるものと確信しています。このプログラムは25の整数倍の行列に適用できるものとなっていますので、使用される時には注意が必要です。