並列行列計算ライブラリ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: [email protected] 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-