並列アルゴリズム
2005年後期 火曜 ２限
高見 利也（青柳 睦）
Aoyagi@cc.kyushu-u.ac.jp
http://server-500.cc.kyushu-u.ac.jp/
12月13日（火）９. ＬＵ分解法とその並列化（講義）PCクラスタによる並列プログラミング（演習）

1

もくじ
1．序

並列計算機の現状

2．計算方式およびアーキテクチュアの分類
3. 並列計算の目的と課題
4. 数値計算における各種の並列化
5. MPIの基礎
6. 並列処理の性能評価
7. 集団通信（Collective Communication）
8. 領域分割（Domain Decomposition）
9. LU分解法とその並列化（講義）PCクラスタによる並列プログラミング（演習）2

LU分解について

最も基本的な線形方程式: 未知ベクトルを求める方程式
ｎ元１次連立方程式

Ax = bAn x n 行列(given)⎛ a11⎜⎜ a21⎜a⎜ 31⎜ a41⎜⎜⎜a⎝ n1a12a22a13a23a32a42#
an 2a33a43

a14 " a1n ⎞⎛ x1 ⎞ ⎛ b1 ⎞⎟⎜ ⎟ ⎜ ⎟a24a2 n ⎟⎜ x2 ⎟ ⎜ b2 ⎟a34a3n ⎟⎜ x3 ⎟ ⎜ b3 ⎟⎟⎜ ⎟ = ⎜ ⎟a44a4 n ⎟⎜ x4 ⎟ ⎜ b4 ⎟% # ⎟⎟⎜⎜ # ⎟⎟ ⎜⎜ # ⎟⎟an 4 " ann ⎟⎠⎜⎝ xn ⎟⎠ ⎜⎝ bn ⎟⎠x nベクトル(unknown)
anベクトル(given)b 数値計算のさまざまな場所で現れる。
 どうやって未知ベクトルを求めるか Gauss消去法(直感的な方法。手計算のやり方)
 LU分解法(同じ行列に対して繰り返すときに効率的)
 繰り返し法(行列にゼロ要素が多い場合など)
 etc.
n33

LU分解アルゴリズム （基礎知識）⎛ l11⎜⎜ l21 l22
⎜ll
A = LU = ⎜ 31 32⎜ l41 l42⎜#⎜⎜l
⎝ n1 ln 2

l43

l44

ln 3

ln 4i< jl 21 u 14 + l 22 u 24 = a 24

⇒ u 24i= jl 31 u 13 + l 32 u 23 + l 33 u 33 = a 33

⇒ l3 3i> jl 41 u 12 + l 42 u 22 = a 42

⇒ l 42l33

0 ⎞⎛ u11 u12⎟⎜" 0 ⎟⎜
u220 ⎟⎜⎟⎜
0 ⎟⎜% # ⎟⎟⎜⎜#
0" lnn ⎟⎠⎜⎝ 0

u13
u23
u33

0i −1

u14u1n ⎞⎟u24 " u2 n ⎟u34
u3 n ⎟⎟,
u44
u4 n ⎟% # ⎟⎟0 " unn ⎟⎠iaij = ∑ lik u kj " (i < j )k =1
jaij = ∑ lik ukj " (i ≥ j )
k =1

i −1

(i < j ) " aij = ∑ lik ukj + lii uij ①

uij = (aij − ∑ lik ukj ) / lii

j −1

j −1

k =1

(i ≥ j ) " aij = ∑ lik ukj + lij u jj ②
k =1

k =1

lij = (aij − ∑ lik ukj ) / u jj
k =14

LU分解アルゴリズム
L,Uの対角値
Doolittle型
 クラウト(Crout)型




Lii=1
Uii=1

３通りのアルゴリズム
right-looking法 （ガウス法）
メモリアクセスが広範囲，バンド幅を要求
 left-looking法 （ガウス法その２）
Crout法right-looking法よりもメモリアクセス範囲は狭い
5分解アルゴリズムの内容(Right-looking法）/* Right-looking method

k=0

Crout type (u_ii == 1)
Input

a: matrix, n: matrix size

Output a: matrix */for(k=0;k<n-1;k++)

k=1
k=2

⎛ a00
⎜⎜ a10
⎜ a20
⎜⎜
⎝ a30

a01

a02

a11
a21
a31

a12
a22
a32

a03 ⎞
⎟
a13 ⎟
a23 ⎟
⎟⎟
a33 ⎠{
for(j=k+1;j<n;j++)a[k][j] /= a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<n;j++)a[i][j] -= a[i][k] * a[k][j];}⎛ l00⎜⎜ l10
⎜ l20
⎜⎜
⎝ l30

u01 u02l11
l21l31u12l22
l32

u03 ⎞⎟
u13 ⎟u23 ⎟⎟⎟l33 ⎠6例えば n=4
for(k=0;k<n-1;k++){for(j=k+1;j<n;j++)a[k][j] /= a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<n;j++)a[i][j] -= a[i][k] * a[k][j];
}

uul

↓

↓↓

a12 = a12 / a11
a13 = a13 / a11

k=1

k=0ul

↓lu

↓

↓

↓lu

↓

↓l→

a22 = a22 − a21 ∗ a12u→

a23 = a23 − a21 ∗ a13

a 01 = a 01 / a 00l→

a 02 = a 02 / a 00

a21 = a21 − a20 ∗ a01

l→l→

a 03 = a 03 / a 00

a22 = a22 − a20 ∗ a02

a32 = a32 − a31 ∗ a12l→u→

a23 = a23 − a20 ∗ a03

a33 = a33 − a31 ∗ a13l→

a11 = a11 − a10 ∗ a01l→

a31 = a31 − a30 ∗ a01

u→

a12 = a12 − a10 ∗ a02l→

a32 = a32 − a30 ∗ a02u→

a13 = a13 − a10 ∗ a03l→

a33 = a33 − a30 ∗ a03

k=2

a23 = a23 / a22

a33 = a33 − a32 ∗ a237

i −1(①i < jの場合)uij = (aij − ∑ lik ukj ) / lii

例えば n=4k =00 −1(0 < 1, 2,3)u01 = (a01 − ∑ l0 k uk 1 ) / l00 = a01 / l00
k =0

u02 = a02 / l00

u03 = a03 / l00

1−1(1 < 2)u12 = (a12 − ∑ l1k uk 2 ) / l11 = (a12 − l10u02 ) / l11
k =0
1−1

(1 < 3)u13 = (a13 − ∑ l1k uk 3 ) / l11 = (a13 − l10u03 ) / l11
k =0（kのループでどの位置にある配列要素が参照されているかに注目）
2 −1(2 < 3)u23 = (a23 − ∑ l2 k uk 3 ) / l22 = (a23 − l20u03 − l21u13 ) / l22
k =08j −1(②i ≥ jの場合ujj = 1に注意)lij = (aij − ∑ lik ukj ) / u jj
k =0

0 −1

(3, 2,1, 0 ≥ 0)l00 = (a00 − ∑ l0 k uk 0 ) / u00 = a00 / u00
k =0

l10 = a10 / u00 l20 = a20 / u00 l30 = a30 / u00
1−1(3, 2,1 ≥ 1)l11 = (a11 − ∑ l1k uk 1 ) / u11 = (a11 − l10u01 ) / u11
k =0

l21 = (a21 − l20u01 ) / u11l31 = (a31 − l30u01 ) / u11

2 −1(3, 2 ≥ 2)l22 = (a22 − ∑ l2 k uk 2 ) / u22 = (a22 − l20u02 − l21u12 ) / u22
k =0

l32 = (a32 − l30u02 − l31u12 ) / u22
3−1(3 ≥ 3)l33 = (a33 − ∑ l3k uk 3 ) / u33 = (a33 − l30u03 − l31u13 − l32u23 ) / u33
k =09

アルゴリズムの違い（参考まで）/* Left-looking method(Up-lookuing method)
Crout type (u_ii == 1)
Input

a: matrix, n: matrix size

Output a: matrix */for(i=0;i<n;i++){/* Crout method
Crout type (u_ii == 1)
Input

Output a: matrix */for(k=0;k<n;k++){for(i=k;i<n;i++)for(j=0;j<k;j++)for(k=0;k<i;k++)for(j=k+1;j<n;j++)a[i][j] -= a[i][k] * a[k][j];a[i][k] -= a[i][j] * a[j][k];for(j=k+1;j<n;j++)for(i=0;i<k;i++)for(j=i+1;j<n;j++)a[i][j] /= a[i][i];}a: matrix, n: matrix size

a[k][j] -= a[i][j] * a[k][i];for(j=k+1;j<n;j++)a[k][j] /= a[k][k];}10LU分解の逐次プログラムを並列化する：/* Start up MPI */MPI_Init(&argc, &argv);

void
lu_s(int n, mat_t a){
int i,j,k;/* Find out process rank */MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);for(k=0;k<n-1;k++){for(j=k+1;j<n;j++)a[k][j] /= a[k][k];}}for(i=k+1;i<n;i++)for(j=k+1;j<n;j++)a[i][j] -= a[i][k] * a[k][j];並列化/* Find out number of processes */MPI_Comm_size(MPI_COMM_WORLD, &numprocs);：
genmat(n, a );：lu_p(n, a);：・行についてcyclic またはblock分割を試みる
・自Rankが持っていない要素の位置に注意
・処理が終わったらRank0に集約する
11

ＬＵ分解アルゴリズムの並列化例(R-L法）
void lu_p(int n, mat_t a){
int i,j,k;
int MAP[SIZE];
int blk;/* cyclic decomposition */for(i=0;i<n;i++) MAP[i]= i % numprocs;/* block decomposition
blk=n/numprocs;
for(i=0;i<numprocs;i++)for(j=0;j<blk;j++) MAP[i*blk + j]= i; */for(k=0;k<n-1;k++){if( MAP[k] == my_rank )for(j=k+1;j<n;j++)a[k][j] /= a[k][k];MPI_Bcast(&a[k][k+1],n-k-1,MPI_DOUBLE,
MAP[k],MPI_COMM_WORLD);

for(i=k+1;i<n;i++) {if( MAP[i] == my_rank )for(j=k+1;j<n;j++)a[i][j] -= a[i][k] * a[k][j];} /* end of i-loop */} /* end of k-loop *//* copy a[i][*] of my_rank -> a[i][*] of Rank0 */for(i=0;i<n;i++) {if(MAP[i] > 0) MPI_Send(&a[i][0],n,MPI_DOUBLE,0,i,MPI_COMM_WORLD);}
if(my_rank==0)
for(i=0;i<n;i++) {if(MAP[i] > 0 ) MPI_Recv(&a[i][0],n,MPI_DOUBLE,MAP[i],i,MPI_COMM_WORLD,&status);}} /* End of lu_p */(以下右上に続く)
12実行時間の計測について◆ 所要時間（Turn Around Time）
最初に実行を開始したプロセスの開始時刻から、最後に実行を終了したプロセスの終了時刻までを計測し，所要時間とする．
MPI_Barrier(MPI_COMM_WORLD);
t_start = MPI_Wtime();/*この間に所要時間を測定する処理を記述*/
MPI_Barrier(MPI_COMM_WORLD);
t_stop = MPI_Wtime();printf(“Turn around time =%.16f¥n",t_stop - t_start );MPI_Barrier「バリア同期」と呼ばれ、すべてのプロセスがこれを呼び出すまで各プロセスが待ち合わせる機能を持つ，MPI標準の関数．
13PCクラスタを利用した演習（12月13日）
利用計算機：
Xeon 3.06GHz dual processor 16 台から成るPCクラスター
Memory:4GB/node, Disk: 1.2TB(total)
Login UID:
講義中に配布端末からのLogin 方法：
Windows XP telnet or sshで omega.cc.kyushu-u.ac.jp へ接続
File 転送： Local ⇔ PC Cluster
IE(ブラウザ）からftp://UID:PWD@omega.cc.kyushu-u.ac.jp/
File 編集： LocalのWindows上
メモ帳，ワードパッド等，エディター
File 編集： PC Cluster上vi エディター

14

PCクラスタを利用した演習(cont.)
Login後の設定：
PCクラスタにはScoreという並列環境が入っており，これをactivateするためにコマンドプロンプトから scout –g all と入力すること．
Mpi program のコンパイル＆リンク:
mpicc test.c
Mpi program の実行：
mpirun –np [1-16のノード数] a.out [programの引数]
はじめの一歩：
講義で指定したLU分解の並列計算プログラムを翻訳＆実行する
実行時間の計測：
MPI_Barrier(COMM), MPI_Wtime() 関数を使って，時間を計測できるようにprogramを変更し，np 数を変えて経過時間を計る.
MPI_Wtime() 関数の型は double .
15

演習レポート(12月13日 出題)【課題】
講義で使用した行でcyclic分割したＬＵ分解プログラムを基に，以下について並列化効率の観点から考察しレポートにして，下記の要領でメールで提出してください．
(1) LU分解部分の経過時間を計測できるように書き換え，使用ノード数を変えて計測結果を表またはグラフにまとめる．
(2) 行でBlock分割するようにLU分解プログラムを変更し，上と同様に経過時間を計測し結果を(1)と比較する．
(3) 行でBlock-Cyclic分割するようにLU分解プログラムを変更し，(1)(2)と同様に経過時間を計測し結果を(1)(2)と比較する．
提出先：aoyagi@cc.kyushu-u.ac.jp
Subject: 並列アルゴリズム課題
締め切り：平成１７年１2月27日（火）形式：text, pdf, ps, Word, Excel
ファイルの添付可
16