2016/7/5
スパコンプログラミング（１）、（Ⅰ）
LU分解法（３）
東京大学情報基盤センター 准教授 塙 敏博
2016年7月5日（火） 10:25-12:10
1
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"qNÛàÜ
H$&5!
Y¾×¿Ù <: " ;
#$%&0'D(?ÛRÜ %$)#(*%#)%$
!
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
講義日程（工学部共通科目）
7.
4月19日(今日)： ガイダンス
4月26日
1.
2.
l
並列数値処理の基本演算（座学）
5月10日：スパコン利用開始
3.
9.
5月17日
l
高性能プログラミング技法の基礎１
（階層メモリ、ループアンローリン
グ）
l
5月24日
5.
l
高性能プログラミング技法の基礎2
（キャッシュブロック化）
5月31日
6.
l
行列-ベクトル積の並列化
ＬＵ分解法（２）
７月5日
l
13.
ＬＵ分解法（１）
コンテスト課題発表
6月28日
l
12.
行列－行列積の並列化（２）
6月14日(10:25-12:10)
l
11.
行列-行列積の並列化（１）
6月14日(8:30-10:15)
★大演習室2
l
10.
べき乗法の並列化
6月7日(10:25-12:10)
l
l ログイン作業、テストプログラム実行
4.
6月7日(8:30-10:15)
★大演習室2
l
8.
2
レポートおよびコンテスト課題
（締切：
2016年8月8日（月）24時 厳守
ＬＵ分解法（３）
7月12日
l
新しいスパコンの紹介・お試し、
他
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
LU分解法の演習日程
今週
1.
•
講義＆並列化の検討
今週
2.
•
LU分解法並列化実習
今週
3.
•
LU分解法並列化実習
3
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
講義の流れ
1. 並列化実習の続き
2. 並列化のヒント（その２）の説明
4
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
LU分解並列化のヒント（２）
C言語版
ほぼ解答が載っています
5
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
6
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);
}
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
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);
7
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
8
前進代入部分(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];
}
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
9
前進代入部分(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 --------------------------------------------------- */
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
LU分解並列化のヒント（２）
FORTRAN言語版
ほぼ解答が載っています
10
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
11
LU分解部分(1)
•
c
c
c
c
c
ib = n/numprocs
istart = myid * ib + １
iend = (myid+１)* ib
--- LU decomposition ---------------------do k=１, iend
idiagPE = (k-１) / ib
--- 枢軸列をもつPE
if (idiagPE .eq. myid) then
dtemp = 1.0 / A(k, k)
枢軸列の計算
--- 枢軸列の転送
do i=myid+１, numprocs – １
call MPI_Send(A(k,k)), … , MPI_DOUBLE_PRECISION, i, k, MPI_COMM_WORLD,
ierr )
enddo
--- 担当範囲の縮小
istart = k + １
else
--- 枢軸列を持たないPE
call MPI_Recv(A(k,k)), …, MPI_DOUBLE_PRECISION idiagPE, k, MPI_COMM_WORLD,
istatus, ierr)
endif
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
LU分解部分(2)
c
c
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
--- End of k-loop ----------------------------------------- 前進消去にメッセージがかぶらないように同期 ----------------------------call MPI_Barrier(MPI_COMM_WORLD, ierr)
12
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
13
前進代入部分(1)
c
c
c
--- 担当範囲の初期化
istart = myid * ib + １
iend = (myid+１) * ib
--- Forward substitution -------------------- 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からデータを受け取る
if (myid .eq. idiagPE) then
c
--- 対角ブロックをもつPE
do kk=１, ib
c
--- 途中結果が送られてくるため必要な変更点
c(k+kk-１) = b(k+kk-１) + c(k+kk-１)
c
call MPI_Recv(c(k), ib,
&
MPI_DOUBLE_PRECISION,
&
myid-１, k, MPI_COMM_WORLD,
&
istatus, ierr)
---対角ブロックだけ先行計算し値を確定させる
do j=istart, istart+kk-2
c(k+kk-１) = c(k+kk-１) - A(k+kk-１, j ) * c( j )
enddo
enddo
2016/7/5
スパコンプログラミング（１）、（Ⅰ）
前進代入部分(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 --------------------------------------------------------------------
14
2016/7/5
おわり
お疲れ様でした
スパコンプログラミング（１）、（Ⅰ）
15
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
yi?^Û.&]ZÜ
',
-D%.?/?0â ¶µÀ×½
-D#&?
%,
#,
!
=
U¡
GQaÛ1&Ü
$D!%?â½Åº×V%
"#
.,
! Ö¹µ×L½ÄÉÖ¹ÒÍ)m
(D%'?
-,
!
%$,
4kÉÖ¹ÒÌ×¹7N¡ \Þ
Û-ÎÏÓÔÙÉ´×ÖÙÓ×
¹Ü
(D#-?
(,
!
4kÉÖ¹ÒÌ×¹7N¡
Û·ÐÁ»ÑÈÖÁ¸Ü
(D1%?
&,
!
m*ʸÄÔ_¡
&D'?/2)1$*%$)%(0
$Qj*#
!
2,
!
m*m_¡ÛÞÜ
!
mÝm_¡ÛßÜ
&D%-?/2)1$*%$)%(0
$Qj*#
&D%-?/%$)#(*%#)%$0
!
%#,
äåqNÛßÜ
áD(?
!
%1,
äåqNÛÞÜ
º×½ÄxXn
&D#2?
!
\#
¦N¡
&D'?/%$)#(*%#)%$0
!
%%,
&
ÕËÙÄ£º×½Äx
Ûgâ
&%!'0(D(?ÛDÜ&)A '
äåqNÛàÜ
'D%#?
!
>½Åº×¡cØu
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"qN¡Qj?^
%,
3
#,
yiÚ¡Kt
3
1,
!"qN)j
3
!"qN)j
"
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
yi¡O±
%, )j¡f
#, ¡Æ×ÄÛ¡ßÜ¡w@
)
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"q¡Æ×ÄÛßÜ
4rvT
§¨q`|©
$
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
'
!"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+F8J{³¬NO8F+
EBG;< 78%,$8+8PQJRQJRA
J{¡sa6:IQ R¥J{³ºÇÙã
I>=8/57;D5EHÞA8885K9:;<=>?@A8885HH088L88+F8J{¡z} F+
SNTUVG9E/W6:IQXRY88X8Y88SNTUZ["\!OY885Y88JY88SNTU4[SSU][^!Z0A
_8
5@BC=B 78JHÞA +F882b¡h, F+
_8G`@G8L888+F8J{³9NO8F+
SNTU^G?a/W6:IQXRY88XY88SNTUZ["\!OY885E5CMNOY88JY8SNTU4[SSU][^!ZY88
W5@BCB:@0A
_
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"q/#0
+F88P F+
I>=8/b7JHÞA88bK9A88bHH08L
EBG;< 786:IQbRA
I>=8/575@BC=BA885K5G9EA885HH08L
PQbRQ5R878PQbRQ5R8* PQJRQ5RFEBG;<A
_8
_
_8+F8O9E8>I8J*`>><8*************************************** F+
+F88P ÎÁ¾Ù¼¤® F **************************** F+
SNTU\C==5G=/SNTU4[SSU][^!Z0A
*
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
(
/%0
3
5@BC=B 78;D5E F856A 5G9E 78/;D5EHÞ08F856A8888+F882b¡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+F882°ÈÖÁ¸° 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#BSF+
I>=8/b75@BC=BA888bK5@BC=BHJJA888bHH0
?QJHJJR8*788PQJHJJRQ8b8R8F8?QbRA
_
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
+
/#0
_8G`@G8L +F8+pÈÖÁ¸³9NO8 F+
+F88l¡6Eb¡ÃÙ¿¡ªsaÛ©CdeI¢Ü 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+
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"q¡Æ×ÄÛßÜ
c[^g^PhrvT
§¨q`|©
!%
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!!
!"q/%0
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{³9NO
?C``8SNTU^G?a/P/JYJ00Y88XY88SNTUZ["\!OUN^O4TVT[h885E5CMNOY88JY88SNTU4[SSU][^!ZY88
5@BCB:@Y85G==0
G9E5I
3
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"q/#0
?88888888888*** P
E>8b75@BC=BY885G9E
EBG;< 78P/8JY8b80
E>857JHÞY89
P/5 Y8b0878P/5 Y8b08k P/5 Y8J08F8EBG;<
G9EE>
G9EE>
G9EE>
?8888888888*** O9E8>I8J*`>><8***************************************
?8888888888*** P ÎÁ¾Ù¼¤® F *****************************
?C``8SNTU\C==5G=/SNTU4[SSU][^!ZY885G==0
!&
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!"
/%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®ÃÙ¿³°
5I8/;D5E ,Gi,885E5CMNO088BjG9
?8888888888888*** +pÈÖÁ¸³¬NO
E>8JJ7ÞY8856
?888888888888888*** ~eI}®±°«3o#BS
?/JHJJ*Þ08786/JHJJ*Þ08H8?/JHJJ*Þ08
?888888888888888***+pÈÖÁ¸msa
³[(°
?C``8SNTU^G?a/?/J0Y856Y8
W888888888888SNTUZ["\!OUN^O4TVT[hY8
W888888888888;D5E*ÞY8JY8SNTU4[SSU][^!ZY8
W8888888888885@BCB:@Y85G==0
E>8b75@BC=BY85@BC=BHJJ*#
?/JHJJ*Þ0878?/JHJJ*Þ08* P/JHJJ*ÞY88b808F8?/8b80
G9EE>
G9EE>
#$%&+'+(
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
/#0
G`@G8
?8888888888*** +pÈÖÁ¸³9NO
E>8JJ7ÞY8856
E>8b75@BC=BY885G9E*Þ
?/JHJJ*Þ08788?/JHJJ*Þ08k P/JHJJ*ÞY88b808F8?/8b80
G9EE>
G9EE>
?8888888888*** l¡6Eb¡ÃÙ¿¡ªsaÛ©CdeI¢Ü
5I8/;D5E ,9G,889:;<=>?@*Þ08 BjG9
?88888888888888***
¡NO l¡82b¡ÃÙ¿³VQaeI³}°
?C``8SNTUVG9E/?/J0Y856Y8SNTUZ["\!OUN^O4TVT[hY88;D5EHÞY8
W8888888888888888888JY8SNTU4[SSU][^!ZY8 85G==0
G9E5I
G9E5I
G9E5I
?88888*** O9E8>I85I882°ÈÖÁ¸° ***************************************
G9EE>
?88*** O9E8>I8J*`>><8********************************************************************
!)
#$%&+'+(
²¯
W±M
½Åº×ÉÖ¹ÒÌ×¹ÛÞÜÛÜ
!$
© Copyright 2024 ExpyDoc
