東京大学情報基盤センター准教授塙敏博 LU 分解法（３）

(1)

LU 分解法（３）

東京大学情報基盤センター准教授塙敏博

2016

年

7

月

5

日（火）

10:25-12:10

(2)

!" qNÛàÜ

H$&5! Y¾×¿Ù <: " ;

#$%& 0 ' D ( ?ÛRÜ %$)#(*%#)%$

(3)

講義日程（工学部共通科目）

1. 4

月

19

日

(

今日

)

：ガイダンス

2. 4

月

26

日

l

並列数値処理の基本演算（座学）

3. 5

^月

10

^{日：スパコン利用開始}

l

ログイン作業、テストプログラム実行

4. 5

月

17

日

l

高性能プログラミング技法の基礎１

（階層メモリ、ループアンローリング）

5. 5

月

24

日

l

高性能プログラミング技法の基礎2

（キャッシュブロック化）

6. 5

月

31

日

l

行列-ベクトル積の並列化

7. 6

月

7

日

(8:30-10:15)

★大演習室

2

l べき乗法の並列化

8. 6

月

7

日

(10:25-12:10)

l 行列-行列積の並列化（１）

9. 6

月

14

日

(8:30-10:15)

★大演習室

2

l 行列－行列積の並列化（２）

10. 6

月

14

日

(10:25-12:10)

l ＬＵ分解法（１）

l コンテスト課題発表

11. 6

月

28

日

l ＬＵ分解法（２）

12.

７月

5

日

l ＬＵ分解法（３）

13. 7

月

12

日

l 新しいスパコンの紹介・お試し、

他

(4)

LU 分解法の演習日程

1. 今週

•

講義＆並列化の検討

2. 今週

• LU

分解法並列化実習

3. 今週

• LU

分解法並列化実習

(5)

講義の流れ

1. 並列化実習の続き

2. 並列化のヒント（その２）の説明

(6)

LU 分解並列化のヒント（２）

C 言語版

ほぼ解答が載っています

(7)

LU 分解部分 (1)

• ib = n/numprocs;

istart = myid * ib;

iend = (myid+ １ )* ib;

/* LU decomposition --- */

for (k=0; k<iend; k++) { idiagPE = k / ib;

if (idiagPE == myid) { /*

枢軸列をもつ

PE */

dtemp = 1.0 / A[k][k];

枢軸列の計算と、

buf[ ]

へ枢軸列をコピー；

for (i=myid+ １ ; i<numprocs; i++) { /*

枢軸列の転送

*/

MPI_Send(&buf[…], … , MPI_DOUBLE, i, k, MPI_COMM_WORLD);

}

istart = k+ １ ; /*

担当範囲の縮小

*/

} else { /*

枢軸列を持たない

PE */

MPI_Recv(&buf[…], …, MPI_DOUBLE, idiagPE, k, MPI_COMM_WORLD,

&istatus);

}

(8)

LU 分解部分 (2)

/*

共通消去部分

*/

for (j=k+ １ ; j<n; j++) { dtemp = buf[j];

for (i=istart; i<iend; i++) { A[j][i] = A[j][i] - A[k][i]*dtemp;

} }

} /* End of k-loop --- */

/*

前進消去にメッセージがかぶらないように同期

--- */

MPI_Barrier(MPI_COMM_WORLD);

(9)

前進代入部分 (1)

• istart = myid * ib; iend = (myid+ １ ) * ib; /*

担当範囲の初期化

*/

/* Forward substitution --- */

for (k=0; k<n; k++)

c[k] = 0.0; /* cの初期化 */

for (k=0; k<n; k+=ib) { /* 対角ブロック判定用ループ */

if (k >= istart) { /* 担当するブロックがある */

idiagPE = k / ib;

if (myid != 0)

/*

左隣りPEからデータを受け取る

*/

MPI_Recv(&c[k], ib, MPI_DOUBLE, myid-１, k, MPI_COMM_WORLD,

&istatus);

if (myid == idiagPE) { /* 対角ブロックをもつPE*/

/* 対角ブロックだけ先行計算し値を確定させる */

for (kk=0; kk<ib; kk++) {

c[k+kk] = b[k+kk] + c[k+kk];/*

途中結果が送られてくるため必要な変更点*/

for (j=istart; j<istart+kk; j++) c[k+kk] -= A[k+kk][ j ] * c[j];

}

(10)

前進代入部分 (2)

} else { /*

対角ブロックを持たない

PE */

/*

自分の所有範囲のデータのみ計算（まだ最終結果ではない）

*/

for (kk=0; kk<ib; kk++) for (j=istart; j<iend; j++)

c[k+kk] -= A[k+kk][j]*c[j];

/* 右隣のPEに、自分の担当範囲のデータを用いた演算結果を送る */

if (myid != numprocs-１)

MPI_Send(&c[k], ib, MPI_DOUBLE, myid+１, k, MPI_COMM_WORLD);

}

} /* End of if(

担当するブロックがある

) --- */

} /* End of k-loop --- */

(11)

LU 分解並列化のヒント（２）

FORTRAN 言語版

ほぼ解答が載っています

(12)

LU 分解部分 (1)

•

ib = n/numprocs istart = myid * ib + １ iend = (myid+１)* ib

c --- LU decomposition --- do k=１, iend

idiagPE = (k-１) / ib c ---

枢軸列をもつPE

if (idiagPE .eq. myid) then dtemp = 1.0 / A(k, k)

枢軸列の計算

c ---

枢軸列の転送

do i=myid+１, numprocs – １

call MPI_Send(A(k,k)), … , MPI_DOUBLE_PRECISION, i, k, MPI_COMM_WORLD, ierr )

enddo

c ---

担当範囲の縮小

istart = k + １ else

c ---

枢軸列を持たないPE

call MPI_Recv(A(k,k)), …, MPI_DOUBLE_PRECISION idiagPE, k, MPI_COMM_WORLD, istatus, ierr)

endif

(13)

LU 分解部分 (2)

c ---

共通消去部分

do j=istart, iend dtemp = A( k, j ) do i=k+ １ , n

A(i , j) = A(i , j) – A(i , k) * dtemp enddo

enddo enddo

c --- End of k-loop ---

c ---

前進消去にメッセージがかぶらないように同期

---

call MPI_Barrier(MPI_COMM_WORLD, ierr)

(14)

前進代入部分 (1)

c ---

担当範囲の初期化

istart = myid * ib + １ iend = (myid+ １ ) * ib

c --- Forward substitution --- c --- c

の初期化

do k= １ , n

c[k] = 0.0 enddo

c ---

対角ブロック判定用ループ

do k= １ , n, ib

if (k .le. istart) then idiagPE = (k- １ ) / ib

c ---

担当するブロックがある

if (myid .ne. 0) then

c ---

左隣りPEからデータを受け取る

call MPI_Recv(c(k), ib,

& MPI_DOUBLE_PRECISION,

& myid- １ , k, MPI_COMM_WORLD,

& istatus, ierr)

if (myid .eq. idiagPE) then c ---

対角ブロックをもつPE

do kk= １ , ib

c ---

途中結果が送られてくるため必要な変更点

c(k+kk- １ ) = b(k+kk- １ ) + c(k+kk- １ )

c ---対角ブロックだけ先行計算し値を確定させる

do j=istart, istart+kk-2

c(k+kk- １ ) = c(k+kk- １ ) - A(k+kk- １ , j ) * c( j ) enddo

enddo

(15)

前進代入部分 (2)

else

c ---

対角ブロックを持たないPE

do kk=１, ib

do j=istart, iend-１

c(k+kk-１) = c(k+kk-１) –A(k+kk-１, j ) * c( j ) enddo

enddo

c ---

自分の所有範囲のデータのみ計算（まだ最終結果ではない）

if (myid .ne. numprocs-１) then

c ---

右隣のPEに、自分の担当範囲のデータを用いた演算結果を送る

call MPI_Send(c(k), ib, MPI_DOUBLE_PRECISION, myid+ １ ,

& k, MPI_COMM_WORLD, ierr) endif

endif endif

c --- End of if 担当するブロックがある ---

enddo

c --- End of k-loop ---

(16)

おわり

お疲れ様でした

(17)

yi?^Û .&]ZÜ

%, - D %. ? / ? 0 â ¶µÀ×½

#, - D #& ?

! =U¡ GQaÛ1&Ü

"# $ ^D !% ^{?â½Åº×V%}

! Ö¹µ× LÂ½ÄÉÖ¹ÒÍ)m -, ( D %' ?

! 4kÉÖ¹ÒÌ×¹7N¡ \Þ Û-ÎÏÓÔÙÉ´×ÖÙÓ×

¹Ü (, ( D #- ?

! 4kÉÖ¹ÒÌ×¹7N¡ \ # Û·ÐÁ»ÑÈÖÁ¸Ü

&, ( D 1% ?

! m * Ê¸ÄÔ_¡

', & D ' ? /2)1$*%$)%(0

$Qj* #

!

¦N¡

2, & D ' ? /%$)#(*%#)%$0

!

m * m_¡ÛÞÜ

., & D %- ? /2)1$*%$)%(0

$Qj* #

!

mÝm_¡ÛßÜ

%$, & D %- ? /%$)#(*%#)%$0

!

äåqNÛÞÜ

!

º×Â½ÄxXn

%%, & D #2 ?

!

äåqNÛßÜ

%#, áD ( ?

!

äåqNÛàÜ

%1, ' D %# ?

!

>½Åº×¡cØu

(18)

!" qN¡Qj?^

%,

3 yiÚ¡Kt

#,

3 !" qN)j

1,

3 !" qN)j

(19)

yi¡O±

%, )j¡f

#, ¡Æ×ÄÛ¡ßÜ¡w@

(20)

!" q¡Æ×ÄÛßÜ

4 rvT

§¨q`|©

(21)

!" q /%0

3 56 789+9:;<=>?@A 5@BC=B 78;D5E F856A 5G9E 78/;D5EH Þ 0F856A

+F8!"8EG?>;<>@5B5>98********************** F+

I>=8/J7$A88JK5G9EA88JHH08L 5E5CMNO 78J8+856A

5I8/5E5CMNO 778;D5E088L88+F8 J{³¬ NO8F+

EBG;< 78%,$8+8PQJRQJRA

J{¡sa 6:IQ R ¥J{³ºÇÙã

I>=8/57;D5EH Þ A8885K9:;<=>?@A8885HH088L88+F8 J{¡z} F+

SNTUVG9E/W6:IQXRY88X8Y88SNTUZ["\!OY885Y88JY88SNTU4[SSU][^!Z0A 5@BC=B 78JH _8 Þ A +F8 82b¡h, F+

_8G`@G8L888+F8 J{³9 NO8F+

SNTU^G?a/W6:IQXRY88XY88SNTUZ["\!OY885E5CMNOY88JY8SNTU4[SSU][^!ZY88 W5@BCB:@0A

_

(22)

!" q /#0

+F88 P F+

I>=8/b7JH Þ A88bK9A88bHH08L EBG;< 786:IQbRA

I>=8/575@BC=BA885K5G9EA885HH08L PQbRQ5R878PQbRQ5R8* PQJRQ5RFEBG;<A _8

_

_8+F8O9E8>I8J*`>><8*************************************** F+

+F88 P ÎÁ¾Ù¼¤® F **************************** F+

SNTU\C==5G=/SNTU4[SSU][^!Z0A

(23)

/%0

3 5@BC=B 78;D5E F856A 5G9E 78/;D5EH Þ 08F856A8888+F8 82b¡F F+

+F8c>=dC=E8@:6@B5B:B5>98****************** F+

I>=8/J7$A8 JK9A88JHH08

?QJR878$,$A888+F8? ¡F F+8

I>=8/J7$A88JK9A88JH756088L8 +F8 +pÈÖÁ¸(VÔÙÉ F+

5I8/J8e785@BC=B088L888+F8 82°ÈÖÁ¸° F+

5E5CMNO 78J8+856A 5I8/;D5E f78$08

+F /¯ NO ®ÃÙ¿³° F+8

SNTU^G?a/W?QJRY856Y8SNTUZ["\!OY8;D5E* Þ Y8JY8SNTU4[SSU][^!ZY8 W5@BCB:@0A

5I8/;D5E 7785E5CMNO088L +F8 +pÈÖÁ¸³¬ NOF+

+F88 +pÈÖÁ¸msa³[(° F+8 I>=8/JJ7$A88JJK56A88JJHH08L

?QJHJJR8786QJHJJR8H8?QJHJJRA+F8 ~eI}®±°«3o#BS F+

I>=8/b75@BC=BA888bK5@BC=BHJJA888bHH0

?QJHJJR8*788PQJHJJRQ8b8R8F8?QbRA

_

(24)

/#0

_8G`@G8L +F8 +pÈÖÁ¸³9 NO8 F+

I>=8/JJ7$A88JJK56A88JJHH0 I>=8/b75@BC=BA888bK5G9EA888bHH0

?QJHJJR8*78PQJHJJRQbRF?QbRA

+F8 ¡ NO l¡82b¡ÃÙ¿³VQaeI³}° F+8 5I8/;D5E f789:;<=>?@* Þ 0

SNTUVG9E/W?QJRY856Y8SNTUZ["\!OY8;D5EH Þ Y8JY8 SNTU4[SSU][^!Z0A

_

_8 +F8O9E8>I85I/ 82°ÈÖÁ¸° 08********************************* F+

_8 +F8O9E8>I8J*`>><8*************************************************** F+

(25)

!" q¡Æ×ÄÛßÜ

c[^g^Ph rvT

§¨q`|©

(26)

!" q /%0

3

56 789+9:;<=>?@

5@BC=B 78;D5E F856 H8 Þ 5G9E 78/;D5EH Þ 0F856

?8888888* !"8EG?>;<>@5B5>98********************

E>8J7 Þ Y885G9E

5E5CMNO 78/J* Þ 08+856

?8888888888*** J{³¬ NO8

5I8/5E5CMNO ,Gi,88;D5E08BjG9 EBG;< 78%,$8+8P/JY8J0 J{¡sa

?88888888888888*** J{¡z}

E>857;D5EH Þ Y889:;<=>?@ k Þ

?C``8SNTUVG9E/P/JYJ00Y88X8Y88SNTUZ["\!OUN^O4TVT[hY8885Y88JY88SNTU4[SSU][^!ZY88 5G== 0

G9EE>

?8888888888888*** 82b¡h,

5@BC=B 78J8H8 Þ G`@G8

?88888888888*** J{³9 NO

?C``8SNTU^G?a/P/JYJ00Y88XY88SNTUZ["\!OUN^O4TVT[h885E5CMNOY88JY88SNTU4[SSU][^!ZY88 5@BCB:@Y85G==0

G9E5I

(27)

!" q /#0

?88888888888*** P

E>8b75@BC=BY885G9E EBG;< 78P/8JY8b80 E>857JH Þ Y89

P/5 Y8b0878P/5 Y8b08k P/5 Y8J08F8EBG;<

G9EE>

?8888888888*** O9E8>I8J*`>><8***************************************

?8888888888* P ÎÁ¾Ù¼¤® F ***************************

?C``8SNTU\C==5G=/SNTU4[SSU][^!ZY885G==0

(28)

/%0

?8888*** 82b¡F

5@BC=B 78;D5E F856 H8 Þ 5G9E 78/;D5EH Þ 08F856

?8888* c>=dC=E8@:6@B5B:B5>98****************

?8888*** ?8 ¡F

E>8J7 Þ Y889

?QJR878$,$ G9EE>

?888*** +pÈÖÁ¸(VÔÙÉ

E>8J7 Þ Y889Y8856

5I8/J8,`G,885@BC=B08BjG98 5E5CMNO 78/J* Þ 08+856

?88888888*** 82°ÈÖÁ¸°

5I8/;D5E ,9G,8$08BjG9

?888888888888*** /¯ NO ®ÃÙ¿³°

?C``8SNTU^G?a/?/J0Y856Y8

W888888888888SNTUZ["\!OUN^O4TVT[hY8 W888888888888;D5E* Þ Y8JY8SNTU4[SSU][^!ZY8 W8888888888885@BCB:@Y85G==0

5I8/;D5E ,Gi,885E5CMNO088BjG9

?8888888888888*** +pÈÖÁ¸³¬ NO

E>8JJ7 Þ Y8856

?888888888888888*** ~eI}®±°«3o#BS

?/JHJJ* Þ 08786/JHJJ* Þ 08H8?/JHJJ* Þ 08

?888888888888888*** +pÈÖÁ¸msa³[(°

E>8b75@BC=BY85@BC=BHJJ*#

?/JHJJ* Þ 0878?/JHJJ* Þ 08* P/JHJJ* Þ Y88b808F8?/8b80 G9EE>

G9EE>

(29)

/#0

G`@G8

?8888888888*** +pÈÖÁ¸³9 NO

E>8JJ7 Þ Y8856

E>8b75@BC=BY885G9E* Þ

?/JHJJ* Þ 08788?/JHJJ* Þ 08kP/JHJJ* Þ Y88b808F8?/8b80 G9EE>

G9EE>

5I8/;D5E ,9G,889:;<=>?@* Þ 08BjG9

?88888888888888*** ¡ NO l¡82b¡ÃÙ¿³VQaeI³}°

?C``8SNTUVG9E/?/J0Y856Y8SNTUZ["\!OUN^O4TVT[hY88;D5EH Þ Y8 W8888888888888888888JY8SNTU4[SSU][^!ZY885G==0

G9E5I G9E5I G9E5I

?88888*** O9E8>I85I882°ÈÖÁ¸° ***************************************

G9EE>

?88*** O9E8>I8J*`>><8********************************************************************

(30)

東京大学情報基盤センター准教授塙敏博 LU 分解法（３）

LU 分解法（３）