並列行列計算ライブラリ ScaLAPACK の VPP700/56 - 九州大学

並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
並列行列計算ライブラリ ScaLAPACKの
VPP700/56における利用法
南里豪志 *
大型計算機センターのスーパーコンピュータVPP700/56に,並列行列計算ライブラリScaLAPACK [3】を導入し
ました. ScaLAPACKは,テネシー大学のJ. Dongarra先生らによって開努されたライブラリで,以下のような特徴
を持っています.
● 豊富なライブラリ群
ScaLAPACKには,合計で約200種類のルーチンが提供されています.また,それらのほとんどで実数,複素数
のデータを扱うことが出来,精度も単精度と倍精度がサポートされています.
● 筒潔なプログラミングインタフェース
ScaLAPACKを用いることにより Fortranに関する知識と本記事で説明する並列処理に関する簡単な知識だ
けで並列プログラムを作成できます.
●高い移植性
ScaLAPACKは, PVM [1]やMPI [2]といった,一般的な通信ライブラリ上にも実装されているため様々な計
算機上で利用できます.
●高い性能
VPP700/56上のScaLAPACKは,富士通が独自にチューニングしているため,ベクトルプロセッサや高速ネッ
トワーク等のハードウェア性能を生かした並列計算を行えます.
● オープンソース
ScaLAPACKのソースコードは,以下のURLから入手できます.
http : //www.netlib. org/scalapack/
本記事では,このScaLAPACKの概要と簡単な利用法を紹介します.なおこの記事は,読者がFortranによるプロ
グラミングを行えることを前提としています.
1 ScaLAPACKとは
テネシー大学のJ. Dongarra先生率いるICL (Innovative Computing IJab.)1では,長年にわたって行
列計算ライブラリの標準化が行われてきました.その成果として公開されたライブラリのうち,逐次型計算
梯及び共有メモリ型並列計算機用のものがLAPACK [4, 5】,分散メモリ型並列計算機用に並列化されたもの
がScaLAPACKです. ScaLAPACKは,以下のようなライブラリを用いて構築されています.
メッセージパッシングライブラリ;
プロセッサ間の通信を行うルーチンを提供するライブラリ.
BLACS:
行列の転送に特化したプロセッサ間通信ライブラリ.
PBLAS:
BLACSのインタフェースを用いた,置換,転置,行列積,ベクトル行列積等の基本的な行列計算を分散
メモリ型並列計算機上で効率良く行うためのルーチンを提供するライブラリ.
これらは,実装するシステムのネットワークやプロセッサの種類,キャッシュメモリの大きさなどに合わせ
てチューニングされます.これらのライブラリが兵装環境の違いを隠蔽してくれるのでSealJAPACKはど
の計算機上でも同じインタフェースと高い性能を提供します.
'大型計算機センター・研究開発部E-mail: nanri@cc.kyushu-u.ac.jp
http://icl.cs.utk.edu/
-139-
九州大学大型計算機センター広報
Vo.32 No.3 1999
解 説
2 基本操作
本節ではScaLAPACKの基本的な利用法を説明します.まずScaLAPACKによる並列計算の基本的な
概念を説明し,次にプログラム例を用いてScaLAPACKを用いたプログラムの構成,コンパイルの方法及び
実行方法を紹介します.
2.1 ScaLAPACKを向いた並列計算
並列プログラムの起動: SPMI)モデル
ScaLAPACKを用いたプログラムは,起動時に利用する全てのプロセッサにコピーされ,それぞれ独立し
たプロセスとして実行されます.このような,同じプログラムが並列計算機の複数のプロセスとして独立し
て動作するモデルはSPMD (single program multiple data)モデルと呼ばれています. VPP700/56では,
利用するプロセッサ数と同じ数のプロセスが起動し,それぞれが同じプログラムを実行します.
処理の割り当て:一般のSPMDモデルの場合
一般にSPMDモデルの各プロセスは,番号等によって一意に区別できるようになっており,このプロセス
番号に応じて処理を分割することにより並列実行を行います
例えば,次のようなループをSPMDモデルで並列実行してみます.ただし,使用するプロセッサ数を4と
します.
例2.1 (ループ例(並列化前))
このままだと全てのプロセッサでI = 1-100のループを実行してしまいます.これではせっかく並列化
しても意味がありません.
PO
1m
P
I
P2 P3
.ll
…75
fit
100
逐次実行
SPMD並列実行
図1: SPMDモデルの例
そこで,各プロセスのプロセス番号MYPROCを用い,以下のようにループを書き換えてみます.使用するプ
ロセッサ数が4なのでMYPROCの取りうる値は0-3です.
例2.2 (ループ例(SPMDモデル))
DO IOI=MYPROC * 25 + 1, (MYPROC+1) * 25
A(I) = I
10 CONTINUE
九州大学大型計算機センター広報
Vol.32 No.3 1999
- 140
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
こうすると,プロセス0がI= 1-25,プロセス1がI=26-50というように,ループを4分割して
各プロセスがそれぞれ処理を分担することになります(図1).各プロセスは並列に実行できるので,この場
合処理時間は1/4程度になると期待できます2.
処理の割り当て: ScaLAPACKの場合
0 1 2 3
図2:プロセス格子の例(2×4)
ScaLAPACKの場合,プロセスをプロセス格子と呼ばれる2次元の格子状に並べます.そのため,例2・2
ではMYPROCのような1次元値でプロセスを区別しましたが ScaLAPACKはプロセス格子中の行(プ
ロセス行)と列(プロセス列)という2次元値で各プロセスを区別します(図2)・
この格子の形はプログラマが指定します.例えば使用するプロセス数が8の場合,プロセス格子として選
択できる形状は,1×8,2×4,4×2,1×8,の4通りです.最適なプロセス格子の形状は,計算のアルゴリズ
ムに依存します.主なScaLAPACKルーチンについて,プロセス数が4の場合のプロセス格子の最適な形
状が,富士通のマニュアル閏で紹介されています・
このプロセス格子に関する情報は,コンテキストと呼ばれる情報体に格納されます.このコンテキストは,
通常,プログラムの最初に作成されます.作成されたコンテキストの内容はメモリ上のテーブルに格納され,
各コンテキストに固有の識別番号で管理されます.その後,この識別番号を介してScaLAPACKのルーチン
内でコンテキストを参照できます.必要に応じて1つのプログラム中で複数のコンテキストを利用すること
もできます.
配列の割り当て;大域配列と配列記述子
VPP700/56のような分散メモリ型並列計算機では,配列を複数のプロセッサにまたがって配置すること
によって, 1プロセッサでは扱えないような大規模な計算を行います.例えば前述の例2.2で,配列Aを4プ
ロセッサにまたがって配置した場合を考えてみます.この場合,プログラムは以下のようになります.
例2.3 (ループ例(SPMDモデル,配列分割後))
DO 10 I =MYPRQC * 25+ 1, (MYPRQC +1) *25
A(I - (MYPROC* 25)) = I
10 CONTINUE
これにより,各プロセスがそれぞれA(1)-A(25)に対して代入を行います.例2.2では各プロセスが所有
する配列Aは100要素でしたが,例2.3では1/4の25要素になっています.このように,配列をプロセス
に割り当てることによってメモリを節約することができます.
通常のSPMDモデルでは,例2.3に示すように,プロセス毎に実行されるループ範囲の計算や,配列Aの
添字の計算が必要です.しかし,プロセス格子の形状と配列のプロセス-の割り当て方法が分かっていれば,
これらの計算方法は一意に決定できます.
2実際には,各プロセスでのループの実行開始時間のずれや並列化に伴う遅延等のため,処理時間が正確に1/4になることは殆ん
どありません・.
-141-
九州大学大型計算機センター広報
Vol.32 No.3 1999
解 説
そこでScaLAPACKでは,このような計算をルーチンの内部で行うことによって,プログラマの負担を軽
減しています.この計算に必要な,プロセス格子や配列のプロセス-の割り当て方法といった情報は,配列記
述子という情報体として,計算対象となる配列とともにScaLAPACKのルーチンに渡されます.
この配列記述子中にコンテキストの識別番号も格納されますので,その番号を介してプロセス格子を参照
することができます.また,配列のプロセス格子-の割り当て方法としては,配列をプロセス格子に割り当て
るときに用いられるブロックの大きさが格納されます.
ブロックとは長方形の部分配列です. SealノAPACKでは,このブロックを-単位として配列をプロセス格
子に割り当てます.またScaLAPACK (及びLAPACK)の計算ルーチンは,このブロック毎に計算を行う
ことにより,キャッシュメモリのような記憶階層やベクトルプロセッサの効率的な利用を図ります.
このブロックを,行方向,列方向それぞれプロセス格子の各プロセスに割り当てます.例えば25 ×25の2
次元配列を,4×4のブロック毎に2×4のプロセス格子の各プロセスに割り当てるよう指定すると,ブロッ
クの割り当ては図3のようになります.
大域配列
ローカル配列
図3:大域配列の例
このように,複数のプロセスに割り当てられる配列を大域配列と呼びます.これに対して,大域配列の実体
を格納する各プロセスの部分配列をローカル配列と呼びます.すなわちScaLAPACKの大域配列は,実体を
格納するローカル配列と,プロセス-の割り当て方法を記述した配列記述子によって定義されます.
ここで,行方向,列方向それぞれについて,大域配列の大きさがブロックの大きさで割りきれない場合や,
ブロックの大きさがプロセス数で割りきれない場合,各プロセスに割り当てられる配列の大きさは,図3に
示すように不均等になります.しかし,プロセス毎に異なるローカル配列を宣言すると配列の扱い方が複雑
になるため,通常各プロセスが同じ大きさのローカル配列を持つようにします・すなわち> Po,Oには必ず最も
大きな配列が割り当てられるので,全てのプロセスでPo,Oに割り当てられる配列と同じ大きさのローカル配
列を宣言します.
2.2 ScaLAPACKの利用法
ここでは,図4に示すプログラムを使いScaLAPACKを用いたプログラムの構成と,コンパイル,実行
方法を紹介します・このプログラムはScaLAPACKのルーチンPDSYEVを用いて,倍精度実数の対称行列
Aの固有値と各固有値に対応する固有ベクトルを計算し,それぞれ1次元配列Wと2次元配列Zに格納し
ます. Aの初期値は, Aが正定借対称行列になるようにプログラム中で設定します.計算の行列の大きさは
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
-142-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
9×9で,使用するプロセッサ数は8です.
2.2.1 プログラムの構成
ScaLAPACKを用いたプログラムの基本的な構成は,以下のようになります.
1.変数,配列の宣言
2・コンテキスト(及びプロセス格子)の作成
3.大域配列のプロセス-の割り当て
4.大域配列の初期化
5. ScaLAPACKルーチンの呼出し
6.プロセス格子の解放
図4のプログラムを用いて順に説明します.
0変数,配列の宣言
INTEGER
CONTEXT, I, INFO, J, MYCOL, MYROW, N, NB,
NPCOL, NPROW, MXLLDA, MXLOCC, LWORK
INTEGER
PARAMETER
DESCAC 9 ), DESCZ( 9 )
(M=9,
N=9,
MB=2,
NB=2,
MXLLDA = 5, MXLOCC = 3, NPROW = 2, NPCOL = 4,
令
令
LWORK = 200 )
′t
l
EXTERNAL
t﹂
Z( MXLLDA, MXLOCC
令
EiiZq Llrnu
DOUBLE PRECISION A( MXLLDA, MXLOCC
M ), WORK( LWORK ),
BLACS.EXIT, BLACS_GET, BLACS.GR工DEXIT,
S
BLACS.GRIDINFO , BLACS_GRIDINIT, BLACS_PINFO ,
令
BLACS.SETUP, DESCINIT, PDSYEV
変数,配列は,通常のFortranの配列と同様に宣言します.大域配列A, Zを格納するローカル配列の大き
さとしては,各プロセスで 0,0に割り当てられる分と同じ大きさを宣言します・例えばこの例では,5×3の
配列を宣言しています・これは, 9×9の大域配列を2×4のプロセス格子に割り当てるので,*b,Oに割り当
てられる配列の大きさが5×3であるためです(図5).
他の変数のうちcoNTEXTは2.1節で説明したコンテキストのための変数です.この変数には,プログラム
中で初期化するコンテキスト毎に付けられる識別番号が格納されます DESCA, DESCZは配列記述子を格納
する1次配列です.また woRKは作業飯域としてScaLAPACKのルーチンに渡す配列です.
0コンテキストの作成
Initialize the ScaLAPACK
CALL SL.INIT(CONTEXT, NPROW, NPCOL)
CALL BLACS_GRIDINFO( CONTEXT, NPROW, NPCOL, MYROW, MYCOL )
Bail out if this process is not a part of this context.
IF( MYROW.EQ.-1 ) GO TO 30
SL.INIT
SL_工Ⅳ工Tでは,コンテキストの作成を行います3.作成されたコンテキストには識別番号が付けられ,
内容は内部のテーブルに格納されます. SL.INITの引数には,作成されたコンテキストの識別番号を
3実はSLJENITの内部では,これらの処理をそれぞれBLACS.GET, BLACS_GRIDINITというBLACSのルーチンによって行って
います・これらのルーチンを直接呼び出すことによって,より細かい設定を行っているプログラムもあります.しかし,少なくとも
VPP700/56で利用する場合SLJNITで十分なので,本記事ではSLJNITを用いた初期設定のみを扱います.
-143-
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
解 説
PROGRAM SAMPLE_PDSYEV_CALL
INTEGER CONTEXT, I, INFO, J, MYCOL, MYROW, N, NB,
NPCOL, NPROW, MXLLDA, MXLOCC, LWORK
令
INTEGER
DESCA(
50
),
DESCZ(
50
)
PARAMETER (M=9, N=9, MB=2, NB=2,
MXLLDA = 5, MXLOCC = 3, NPROW = 2, NPCOL = 4,
i
H
U
s
'
Z( MXLLDA, MXLOCC
) )
LWORK = 200 )
DOUBLE PRECISION A( MXLLDA, MXLOCC
M ), WORK( LWORK ),
EXTERNAL BLACS_EXIT, BLACS_GET, BLACS.GRIDEXIT,
BLACS_GRIDINFO , BLACS.GRIDINIT, BLACS_PINFO ,
BLACS.SETUP, DESCINIT, PDSYEV
Initialize the ScaLAPACK
CALL SL.INIT(CONTEXT, NPROW, NPCOL)
CALL BLACS_GRIDINFO( CONTEXT, NPROW, NPCOL, MYROW, MYCOL )
Bail out if this process is not a part of this context.
IF( MYROW.EQ.-1 ) GO TO 30
basic array descriptor
CALL DESCINITC DESCA, M, N, MB, NB, 0, 0, CONTEXT, MXLLDA, INFO )
CALL DESCINITC DESCZ, M, N, MB, NB, 0, 0, CONTEXT, MXLLDA, INFO )
Build a matrix that you can create with
I a one line matlab command: hilb(n) + diag([l:-1/n:1/n])
DO 20J = 1, N
DO 10 I = 1, M
IF( I.EQ.J ) THEN
CALL PDELSETC A, I, J, DESCA,
(DBLE(N-I+l ) ) /DBLE(N )+DBLE( 1 ) /
(DBLE( I+J )-DBLE( 1 ) ) )
ELSE
CALL PDELSETC A, I, J, DESCA,
DBLE( 1 ) / (DBLE(I+J )-DBLE( 1 ) ) )
END IF
10 CONTINUI:
20 CONTINUE
Ask PDSYEV to compute the entire eigendecomposition
CALL PDSYEVC 'V, 'U¥ N, A, 1, 1, DESCA, W, Z, 1, 1,
DESCZ, WORK, LWORK, INFO )
CALL BLACS_GRIDEXIT( CONTEXT )
30 CONTINUE
CALL BLACS.EXIT( 0 )
STOP
END
図4:サンプルプログラム
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
- 144-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
格納する整数変数,及び,作成するコンテキスト内のプロセス格子の行数,列数を指定します.この例
では, NPROW- 2, NPCOL- 4なので, 2×4のプロセス格子を持つコンテキストを作成し,その識別番
号がCoNTEXTに格納されます.
BLACS_GRエDINFO
通常,コンテキストの作成後,各プロセスはBLACS_GRIDINFOを用いてプロセス格子の情報を取得し
ます.これにより,各プロセスがプロセス格子中の自分の位置(行,列)を知ることが出来ます.この例
では,それぞれMYROW, MYCOLに格納されます.またBLACS_GRIDINFOの第2引数以降が-1になる
場合はSL.INITに失敗していることを示します.この例では,このような場合処理を終了するように
指示されています.
0大域配列のプロセス-の割り当て
0 1 2 3 0
Pqo **ol ^02 03
^10 *11 ^12 *13
(b)プロセス格子上の配列Aの内容
(a)大域配列のプロセス格子への割り当て
図5:ブロックサイクリック分割の例
0
0
0
0
t一 一
一 '
I 一
B
s
' I
' _ 一
B
B
Ⅳ⋮H▲
B
s
' I
as sc
CALL DESC工NIT( DESCZ,
Ⅳ Ⅳ
CALL DESCINITC DESCA,
CONTEXT, MXLLDA,
CONTEXT, MXLLDA,
DESCINIT
大域配列のブロックの大きさと,プロセス-の割り当て方法を配列記述子に格納します.このDESC工NIT
は,密行列の割り当て方法を指定する配列記述子を作成するルーチンです.ここでは,
という割り当て方法を格納した配列記述子DESCA, DESCZを作成しています.
-145-
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
解 説
この配列記述子による大域配列Aの割り当ては図5(a)のようになります. Aはまず2×2のブロック
に分けられ,ブロック毎に2×4のプロセス格子に割り当てられます.各ブロックは,行方向,列方向
それぞれ0番目のプロセスから順に割り当て,最後のプロセスに割り当てたらまた0番目から繰り返
すという,サイクリック方式でプロセス格子に割り当てられます・その結果,各プロセス(Polo-Pl,3)
のローカル配列にはそれぞれ図5(b)に示すような要素が格納されます・
ブロックの大きさ:
分割に用いるブロックの最適な大きさは,プログラムのアルゴリズムや計算機のアーキテクチャ(キャッ
シュサイズやベクトル長)によって異なります・あまり大きすぎると, LU分解のような計算が進むに
つれて計算範囲が小さくなるようなアルゴリズムでプロセス毎の計算量に偏りが生じます.逆に小さ
すぎると,キャッシュメモリやベクトルプロセッサの性能を生かせなくなります.富士通のマニュア
ル[71によるとVPP700/56では,ブロックサイズが64 × 64の場合に殆んどのルーチンで良い性能
が得られるようです.
0配列の初期化
DO 20J = 1, N
DO 10 I = 1, M
IF( I.EQ.J ) THEN
CALL PDELSETC A, I, J, DESCA,
(DBLE(N一工+1
)
)
/DBLE(N)+DBLEC
1
)
/
(DBLE(I+J )-DBLE( 1 ) ) )
ELSE
CALL PDELSETC A, I, J, DESCA,
DBLE( 1 ) / (DBLEC I+J )-DBLE( 1 ) ) )
END IF
10 CONTINUE
20 CONTINUE
ここでは大域配列Aの各要素に値を代入しています.値は, Aが正定値対称行列になるように計算してい
ます.なお, DBLEQは与えられた整数と同じ大きさの実数を返す関数です.
PDELSET
初期化のように大域配列の各要素に個別にアクセスする場合,その要素を所有するプロセスの格子内
での位置と,その要素のローカル配列内での位置を計算する必要があります.例えば図5に示すよう
に,このプログラム中の大域配列Aの要素A(7, 7)はプロセス格子中のプロセス1,3が所有し,その
ローカル配列A内の位置はA(3, 1)です.
しかし,このような計算を大域配列の要素-のアクセス毎に行うのは面倒です.そこでこのプログラム
ではPDELSET用いて簡潔に配列の初期値を設定しています. PDELSETは,与えられた配列記述子か
ら大域配列の割り当て情報を参照し,指示された要素に値を格納するルーチンです.ここでは,配列記
述子DESCAで指示された割り当て方法でプロセス格子に割り当てられた行列Aの要素(I, J)に値を
格納しています.
大域配列の操作方法については 3.5節で説明します.
O ScaLAPACKルーチンの呼出し
CALL
PDSYEVC
'V,サU*,
N,
A,
1,
1,
DESCA,
W,
DESCZ, WORK, LWORK, INFO )
九州大学大型計算機センター広報
Vol. 32 No.3 1999
-146-
Z,
1,
1,
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
PDSYEV
このプログラムでは,対称行列の固有値求解ルーチンPDSYEVを用いています.計算の対象となる行列
は,配列記述子DESCAに格納された割り当て方法でプロセス格子に割り当てられた配列Aです.計算
した固有値はベクトルWに格納され,また,それぞれの固有値に対応した固有ベクトルは行列Zに格納
されます.このZのプロセス格子-の割り当て方法としては,配列記述子DESCZを渡しています.
このようにScaLAPACKルーチンは,引数として渡された配列記述子から処理すべき計算や通信を算
出し,ローカル配列に対して実行します.ちなみにPDSYEVの最初の引数'V'は固有ベクトルを計算す
ることを指示し, 2番目の引数'U'は対称行列の要素がAの上三角部分に格納されていることを示し
ています4.
作業領域
PDSYEVは,実行時に作業領域を必要とします.このプログラムでは,作業領域として1次元配列woRK,
また,作業領域の大きさとしてLWORKを渡しています.このようにScaLAPACKのルーチンには作
業領域を必要とするものがあります.各ルーチンについて最低限必要な作業領域の大きさは,そのルー
チン自身を用いて知ることができます.すなわち,作業領域を必要とするルーチンに,作業領域の大き
さとして-1を渡し,呼び出すと,そのルーチンの実行に最小限必要な作業簡域の大きさが,作業傍域用
の配列の先頭要素に格納されます.
例えばPDSYEVの場合,
CALL PDSYEVC 'V, >U', N, A, 1, 1, DESCA, W, Z, 1, 1,
DESCZ, WORK, -1, INFO )
とすると, WORK(l)に最小限必要な作業額域の大きさが格納されます.
すなわち,実際の計算を行う前に同じプログラムを使って予め必要な作業領域を調べ,その後LWORK
をその値に書き換えて実行すれば,必要最小限の作業領域で計算を行えます.
0プロセス格子の解放
CALL BLACS.GRIDEX工T( CONTEXT )
CALL BLACS.EXIT( 0 )
BLACS_GRIDEXIT , BLACS_EXIT
SL_エⅣエTで初期化したプロセス格子は,そのプロセス格子が格納されているコンテキストをBLACS
のルーチンBLACS.GR工DEXITに渡すことによって解放されます.その後 BLACS.EX王Tルーチンで並
列実行を終了します.
2.2.2 コンパイル
大型計算機センターのVPP700/56上でScaLAPACKを用いた並列プログラムのコンパイル例を以下に
示します. VPP700/56の利用法については国を参照してください・
例2.4 (ScaLAPACKを用いたプログラムのコンパイル触)
frt pdsyev_ex.f -o pdsyev_ex -Wl,-P,-J -lscalapack -lFblacsMPI -llapackvp -lblasvp \
-L/usr/lang/mpi/lib -lmpi -Imp -lelf -lpx
4このプログラムではAの下三角部分にも債が格納されています.
-147-
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
解 説
この例では, pdsyev.exを実行形式の出力ファイルとしてFortranのソースプログラムpdsyev.ex.fを
コンパイルしています.それ以外の-Wl.-P,-J以降の部分はVPP700/56でScaLAPACKを用いる場合
必ず必要です・なお,オプションが長いため, \で行の継続を示しています・
注意;
ScaLAPACKやPBLAS, BLACSの各ルーチンはFortran77 (及びFortran90)から利用できますが,
Fortran90/VPPの並列化機能を利用した配列は扱えません・
2.2.3 実行
ScaLAPACKを用いた並列プログラムは,他の並列プログラムと同様,バッチジョブとしてVPP700/56
の並列用バッチキューに投入することによって実行します.以下に投入する際に用いるジョブスクリプトの
例を示します.
例2.5(ジョブスクリプトの例)
港
cd EXAMPLE
pdsyev_ex
この例では,ホームディレクトリの下のEXAMPLEというディレクトリにあるpdsyev_exという名前の実
行形式ファイルを実行します.
ファイル名run.sh.に保存されたジョブスクリプトは,以下のようにジョブ投入コマンドqsub利用して
ジョブキューに投入します.
例2.6 (ジョブスクリプトの投入例)
kyu-vpp/, qsub -q s8 run.sh
qsubコマンドに指定可能なオプションのうち,よく用いられるものを以下に示します.
-q queue queueはキュー名です.省略するとVPP700/56のplキューに投入されます.
-ITtime 使用するバッチリクエストのCPU時間の上限を設定します・時間はll時:]分:】秒で指定
してください.省略値は20:00:00です.
-IP n 使用するPE数の上限を設定します.省略値はバッチキューの上限値です.
-1舛mem PEあたりの記憶債域について使用量の上限を設定します.使用量は100kb, 100mb, 1.7gb
のように指定してください.省略値は512mbです.
ここで, -qオプションで指定するキューのうち,並列ジョブを投入できるキューは表1の通りです.この
中で,特にS8キューはCPU時間の制限値が短いため,テストやデバッグに適しています.
表1:並列ジョブのバッチキュー
i:・ff:・
L馴蔓oS乱lf_;;GB
GB
GB製葦至 jPE
PE
PE
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
- 148-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
使用プロセッサ数
qsubで指定する使用プロセッサ数には,プログラム中で実際に利用するプロセス数と同じ数を指定して下
さい.プロセッサ数がプロセス数より少ない場合,プログラムを実行できません.また,プロセッサ数がプロ
セス数より多い場合,プログラムを正しく実行できないことがあります.
3 SealJAPACKにおける大域配列
前節までの内容でScaLAPACKを用いて密行列の行列演算を行う簡単なプログラムを作ることはできま
す.しかしScaLAPACKに用意されているバンド行列を扱うルーチンでは,大域配列の格納方法や扱い方
が密行列と異なります.そこで本節では,行列の種類に応じた大域配列の割り当て方法と操作方法を紹介し
ます.また,大域配列の要素をScaLAPACKのルーチンを介さず直接扱う方法についても説明します.
3.1 配列記述子のフォーマット
前述したように, ScaLAPACKの大域配列は配列記述子によって割り当て方法が記述されます.実は,こ
の配列記述子の実体は1次元配列です.この配列記述子の要素として格納される情報は,大域配列の種類に
よって大きく異なります.そのためScaLAPACKでは,配列記述子のフォーマットとして以下の4種類が
用意されています.
1.密行列
2.バンド行列
3.バンド行列を扱うScaLAPACKルーチンに渡す右辺行列
4.外部記憶装置上に配置された密行列
本記事では,これらのうち1-3について説明します.なお3.の右辺行列とは ScaLAPACKルーチンの
うち1次連立方程式Ax=bを解くルーチンに渡されるベクトルbのことです.
3.2 密行列の割り当てと操作
一般の密行列のための配列記述子のフォーマットは表2に示す通りです.
表2:密行列の配列記述子
要素番号
密行列の配列記述子はDESCエNITを用いて初期化できます.例えば図4のプログラム中の
CALL DESCINITC DESCA, M, N, MB, NB, 0, 0, CONTEXT, MXLLDA, INFO )
-149-
九州大学大型計算機センター広報
Vol.32 No. 3 1999
解 説
は,以下の代入を行うことと等価です.
DESCA(l) = 1
DESCAC2) = CONTEXT
DESCAC3) = M
DESCA(4) = N
DESCA(5) = MB
DESCAC6) = NB
DESCAC7) = 0
DESCAC8) = 0
DESCAC9) = MXLLDA
密行列の扱い方は,前節で説明した通りです.
3.3 バンド行列の割り当てと操作
(a)バンド行列
(b)バンド行列のプロセス-の割り当て
図6: ScaLAPACKにおけるバンド行列の分割
バンド行列の配列記述子のフォーマットは表3に示す通りです.このフォーマットは,一般バンド行列用
ルーチンや三重対角行列用ルーチン等,バンド行列用のScaLAPACKルーチンに渡される大域配列の配列
記述子に利用します.
バンド行列のローカル配列-の格納方法は, 2節で説明した格納方法と大きく異なります.まずScaLAPACKでは,バンド行列を割り当てる対象となるプロセス格子は1 ×Pの形でなければなりません.また, 1
プロセスに割り当てられるのは高々1ブロックです.すなわち,割り当てに用いられるブロックの列数は,行
列の列数をNとすると, 「N/P]です.例えば図6(a)に示すような非対称バンド行列は,列数9,プロセス
数4なので, 「9/41 -3でブロックの列数は3となります.そのためこの行列は図6(b)に示すように各プ
ロセスに格納されます・この例の-0,3のように,列数とプロセス数によっては配列が割り当てられないプロ
セスもあります.
また,これは非対称バンド行列の例でしたが,同じバンド行列でも,行列が対称であるか否か,三重対角で
あるか否か,及び, LU分解時にピボット選択が必要か否かでさらに格納方法が違います.
バンド行列の格納方法や扱い方についての詳細は文献[3】を参照してください.また,この文献は以下の
URIJから参照できます.
http : //www. netlib. org/scalapack/slug/scalapack_slug. html
九州大学大型計算機センター広報
Vol.32 No. 3 1999
-150-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
表3:バンド行列の配列記述子
3.4 バンド行列による1次元連立方程式計算ルーチンに渡される右辺行列
バンド行列用ルーチンのうち, 1次元連立方程式を解くルーチンでは,右辺配列として渡される配列につい
てもプロセス格子-の割り当て方法が故められています.まず,プロセス格子の形状はP x lでなければな
りません.また,バンド行列の割り当てと同様,1プロセスに割り当てられるブロック数は高々1ブロックで
す.バンド行列用ルーチンに渡す右辺配列の配列記述子のフォーマットを表4に示します.
表4:バンド行列用ルーチンに渡す右辺配列の配列記述子
3.5 大域配列の割り当てと操作
配列記述子の情報を用いると,大域配列の要素からその要素を所有するプロセスとそのローカル配列内で
の位置を計算できます.
例えば行数がNの行列を,行方向に,ブロックの大きさNBでPプロセスに割り当てた場合に, Z番目
の要素が割り当てられているプロセスp,プロセス内のブロック番号l及びブロック内の位置3:は,それぞ
れ以下の式で計算できます.
p-1(1-1)/NB¥modP, l-[(I-1)/(P*NB)」 x-((I-1)modNB)+l
列方向についても同様の計算が成り立ちます.そのため,この式を利用して,各プロセスで自分の所有する要
素の操作を行うことができます.
しかし,大域配列の操作の度にこのような計算を行うのは面倒です.そこで, ScaLAPACKには大域配列
を簡単に操作するため以下のようなルーチンが用意されています.
PxELSETC A,工A, JA, DESCA, ALPHA ),
PxELGETC SCOPE, TOP, ALPHA, A,工A, JA, DESCA )
-151
九州大学大型計算機センター広報
Vol.32 No. 3 1999
解 説
ここでPxELSET, PxELGETのⅩには,大域配列Aのデータ型に応じて, I (整数), s (単精度実数), D(倍精度
実数), c (単精度複素数), z (倍精度複素数)が入ります.
PxELSETは,大域配列Aの要素(工A, JA)に値ALPHAを代入するものです. DESCAはAの配列記述子で
す.
またPxELGETは,大域配列Aの要素(工A, JA)の値をALPHAに格納するものです. PxELGETのSCOPE
にはALPHAに値を格納するプロセスの範囲を指定します. SCOPEとして指定できるのは以下の通りですが,
通常は'A'を指定すれば良いでしょう.
SCOPE
プロセス格子の中でA(IA, JA)を所有するプロセスと同じ行にある全プロセス
プロセス格子の中でA(IA, JA)を所有するプロセスと同じ列にある全プロセス
プロセス格子中の全プロセス
A(工A, JA)を所有するプロセスのみ
またTOPには要素を他のプロセスに転送するときの手段を指示します. VPP700/56では' 'を渡します・
4 参考資料
本記事はScaLAPACKを開発したグループの作成したマニュアル[3】及び,富士通が提供する説明書[7】
を主に参考にして作成しました.どちらも大型計算機センター図書室で閲覧できます.特に文献同には,
ScaLAPACKの各ルーチンの利用方法だけでなく, BLACS, PBLASの利用方法も書かれています. ScaLAPACKの各ルーチンの説明はkyu-vppのmanコマンドで参照できるオンラインマニュアルでもある程度
調べることができますが,多少分かりにくいので,可能であれば文献l3]を参照されることをお勧めします・
以下のURLでアクセスできるScaLAPACKの開発元のホームページにも,様々な情報があります.
http : //www.netlib. org/scalapack/
ここでは,文献[3]の一部がHTMLファイルで提供されています・また,以下のURLから取得できるサ
ンプルプログラムはScaLAPACKの勉強に適しています.
http : //www. netlib. org/scalapack/examples/
さらに ScaLAPACKのソースコードもこのホームページから入手できます.ちょっと大変ですが,行列
計算の参考として読んでみるといいかも知れません5.
またSealJAPACKは,アーキテクチャや通信ライブラリ等のシステムに依存する部分をユーザから隠蔽
していますが,実際に利用する上で,システムに依存したトラブルが発生することがあります. VPP700/56
の場合,通信ライブラリにMPIを用いていますので,以下のURLから参照できる MPIの利用に関する注
意事項を参照してください.
http
:
//www.
cc.kyushu-u.
ac.
jp/library/etc/mpi一faq.html
5ただし, VPP700/56上のScaLAPACKはこのソースコードを富士通がチューニングしたものです.
九州大学大型計算棟センター広報
Vol.32 No. 3 1999
-152-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
5 ScaLAPACKのルーチン一覧
ScaLAPACKのルーチンは,一般に問題を解くためのドライバルーチン, LU分解,三重対角化など,個々
の計算を実行する計算ルーチン,及び,それらを補助する補助ルーチンから構成されます.ここでは,そのう
ち,ドライバルーチンと計算ルーチンを列挙します.利用法の詳細は,マニュアル等を参照してください.
5.1 ScaLAPACKドライバルーチン
連立1次方程式
連立1次方程式のドライバルーチンには,通常の単純ドライバに加えて,いくつかの機能を追加したエキ
スパート・ドライバがあります.エキスパートドライバは,条件数の算定,行列の特異性のチェック,解の反
復改良,誤差解析等の機能をサポートします.ただしエキスパート・ドライバを実行するためには,単純ドラ
イバの約2倍の記憶容量が必要です.
なお,表中の(未サポート)は, VPP700/56でまだ実装されていないルーチンを示します.
線形最小二乗問題
標準固有値問題,特異債分解
対称固有値問題の単純ドライバルーチンは,固有値全てと選択された固有ベクトルを計算します.エキス
パート・ドライバは指定した範囲の固有値と対応する固有ベクトルが計算できます.どちらも対称Hermite
行列を扱います.
また,特異値分解のドライバルーチンは一般の矩形行列を対象とします.
型
- 153
九州大学大型計算機センター広報
Vol. 32 No.3 1999
解 説
一般化固有値問題
AとBが対称またはHermite行列で, Bが正定値である場合の固有値問題Az - ABz, ABz- Azま
たはBAz-入Zの解を求めます.
5.2 ScaLAPACK計算ルーチン
連立1次方程式
行列の型と格納形式
一般行列
(partial pivoting)
一般バンド行列
(partial pivoting)
複素数
LU分解
分解を使って求解
条件数を推定する
解の誤差限界を計算する
分解を使って逆行列を求
める
方程式を均衡化する
LU分解
分解を使って求解(未サ
PSGETRF
PS GETRS
PSGECON
PS GERFS
PS GETRI
P CGETRF
PD GETRF
PCGETRS
PD GETRS
PZGETRS
PCGECON
PDGECON
PZGECON
PCGERFS
PCGETRI
PD GERFS
PD GETRI
PZGERFS
P ZGETRI
PSGEEQU
PS GBTRF
PSGBTRS
PCGEEQU
PCGBTRF
PDGEEQU
PD GBTRF
PZGEEQU
PZGBTRF
PCGBTRS
PD GBTRS
PZGBTRS
PSDBTRF
PSDBTRS
PSDTTRF
PSDTTRS
PSPOTRF
PCDBTRF
PDDBTRF
PZDBTRF
PCDBTRS
PDDBTRS
PZDBTRS
PCDTTRF
PDDTTRF
P ZD TTRF
P CD TTRS
PDD TTRS
PZDTTRS
P CPOTRF
PDPOTRF
PZPOTRF
PSPOTRS
PSPOCON
PSPORFS
PSP OTRI
PCPOTRS
PCPOCON
PCPORFS
PCPOTRI
PDPOTRS
PDPOCON
PDPORFS
PDPOTRI
PZPOTRS
PZPOCON
PZPORFS
PZPOTRI
PSPOEQU
PSPBTRF
PSPBTRS
PSPTTRF
PSPTTRS
PS TRTRS
PSTRCON
PS TRRF、S
P S TRTRI
PCPOEQU
PCPBTRF
PCPBTRS
P CPTTRF
PCPTTRS
PDPOEQU
PDPBTRF
PDPBTRS
PDPTTRF
PDPTTRS
P CTRTRS
PCTRCON
PD TRTRS
PDTRCON
P CTRRFS
PD TRRFS
P CTRTRI
PDTRTRI
PZPOEQU
PZPBTRF
PZPBTRS
P ZPTTRF
PZPTTRS
P ZTRTRS
PZTRCON
P ZTRRFS
P ZTRTRI
PZGETRF
ポート)
一般バンド行列
(no pivoting)
一般三重対角行列
(no pivoting)
対称/ Hermite正定値行
列
対称/ Hermite
正定借バンド行列
対称/ Hermite
正定借三重対角行列
三角行列
LU分解
分解を使って求解
LU分解
分解を使って求解
Cholesky分解
分解を使って求解
条件数を推定する
解の誤差限界を計算する
分解を使って逆行列を求
める
方程式を均等化する
Cholesky分解
分解を使って求解
LDL分解
分解を使って求解
LU分解
条件数を推定する
解の誤差限界を計算する
逆行列を求める
九州大学大型計算機センター広報
Vol. 32 No.3 1999
-154-
並列行列計算ライブラリScaLAPACKのVPP700/56における利用法
直交分解と線形最小二乗問題
一般化直交分解と線形最小二乗問題
対称固有値問題
非対称固有値問題
- 155
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
解 説
特異値分解
対称正定借一般化固有値問題
5.3 VPP700/56における未サポートルーチン
VPP700/56上のScaLAPACKでは,以下のルーチンはサポートされていません・
PxGBSV, PxGBTRF, PxGBTRS, PxGESVD
ただし, Xには,ルーチンのデータ型に応じて, S (単精度実数), D(倍精度実数). C (単精度複素数), Z(倍
精度複素数)が入ります.
6 むすび
本記事では,並列行列計算ライブラリScaLAPACKを紹介しました・大型計算機センターのVPP700/56
には,同様のライブラリとして既にSSLII/VPPが導入されています・こちらはVPP700/56のハードウェ
アを作った会社によるものなので,性能については(きちんと比較していませんが) ScaLAPACKを上回る
と思います.これに対してScaLAPACKは PVMもしくはMPIが動作すれば簡単に導入できるので,
VPP700/56以外の様々な計算機上でも同じプログラムを利用できるという利点があります・そのため,と
にかくVPP700/56で簡単に高性能を得たい場合はSSL II/VPP,他の計算機でも同じプログラムを利用し
たい場合はScaLAPACKというように,用途に応じて使い分けてください.
参考文献
[1] http : //www. epm. ornl. gov/pvm/pvm_home.html
[2] http : //www.mpi-fornm. org/
[3] L.S. Black ford, J. Choi, A. Cleary, E. D'Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet, K. Stanley, D. Walker and R.C. Whaley, ScaLAPACK Users'Guide,
SIAM, Philadelphia, ISBN 0-89871-397-8, 1997.
[4] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, A. Hammarling, A. McKenney, S. Ostrouchov and D. Sorensen,小国力訳,行列計算パッケージLAPACK
利用の手引,丸善ISBN 4-621-04076-6, 1995.
[51渡部善隆,行列計算パッケージLAPACK/VPの紹介,九州大学大型計算機センター広報, Vol. 32, No.
1, pp.ト10, (1999).
囲VPP700/56利用の手引(第2版),九州大学大型計算機センター, (1998).
閏UXP/V ScaLAPACK VIOソフトウェア説明書,富士通, (1999).
九州大学大型計算機センター広報
Vol. 32 No. 3 1999
-156-