LU 分解法（１）
東京大学情報基盤センター 准教授 塙 敏博
2020年6月23日（火）10:25-12:10
Site Computer/Year Vendor Cores Rmax
(TFLOPS)
Rpeak
(TFLOPS)
Power (kW) 1 Fugaku, 2020, Japan
R-CCS, RIKEN
Fujitsu PRIMEHPC FX1000, Fujitsu A64FX 48C
2.2GHz, Tofu-D 7,299,072 415,530
(= 415.5 PF) 513,854.7 28,335 2 Summit, 2018, USA
DOE/SC/Oak Ridge National Laboratory
IBM Power System AC922, IBM POWER9 22C 3.07GHz, NVIDIA Volta GV100, Dual-rail Mellanox EDR Infiniband
2,414,592 148,600 200,795 10,096
3 Sierra, 2018, USA DOE/NNSA/LLNL
IBM Power System S922LC, IBM POWER9 22C 3.1GHz, NVIDIA Volta GV100, Dual-rail Mellanox EDR Infiniband
1,572,480 94,640 125,712 7,438
4 Sunway TaihuLight, 2016, China National Supercomputing Center in Wuxi
Sunway MPP, Sunway SW26010 260C 1.45GHz,
Sunway 10,649,600 93,015 125,436 15,371
5 Tianhe-2A, 2018, China
National Super Computer Center in Guangzhou
TH-IVB-FEP Cluster, Intel Xeon E5-2692v2 12C
2.2GHz, TH Express-2, Matrix-2000 4,981,760 61,445 100,679 18,482 6 HPC5, 2020, Italy
Eni S.p.A.
Dell C4140, Xeon Gold 6252 24c 2.1GHz, NVIDIA
Volta GV100, Mellanox Infiniband HDR 669,760 35,450 51,720 2,252 7 Selene, 2020, USA
NVIDIA
NVIDIA DGX A100 SuperPOD, AMD EPYC 7742 64C 2.25GHz, NVIDIA GA100, Mellanox Infiniband HDR
277,760 27,580 34,568 1,344 8 Frontera, 2019, USA
Texas Advanced Computing Center
Dell C6420, Xeon Platinum 8280 28c 2.7GHz,
Mellanox Infiniband HDR 448,448 23,516 38,746
9 Marconi-100, 2020, Italy Cineca
IBM Power System AC922, IBM POWER9 22C 3.07GHz, NVIDIA Volta GV100, Dual-rail Mellanox EDR Infiniband
347,776 21,640 29,354 1,476
10 Piz Daint, 2017, Switzerland Swiss National Supercomputing Centre (CSCS)
Cray XC50, Xeon E5-2690v3 12C 2.6GHz, Aries
interconnect , NVIDIA Tesla P100 387,872 21,230 27,154 2,384
19 Oakforest-PACS, 2016, Japan Joint Center for Advanced HPC
PRIMERGY CX1640 M1, Intel Xeon Phi 7250 68C
1.4GHz, Intel Omni-Path 556,104 13,556 24,913 2,719
2020/6/23
http://www.hpcg-benchmark.org/
Computer Cores HPL Rmax
(Pflop/s)
TOP500 Rank
HPCG (Pflop/s)
1 Fugaku 7,299,072 415.530 1 13.400
2 Summit 2,414,592 148.600 2 2.926
3 Sierra 1,572,480 94.640 3 1.796
4 HPC5 669,760 35,450 6 0.860
5 Trinity 979,072 20.159 11 0.546
6 Selene 277,760 27.580 7 0.509
7 ABCI 391,680 19.880 12 0.509
8 Piz Daint 387,872 21.230 10 0.497
9 Sunway TaihuLight 10,649,600 93.015 4 0.481
10 Nurion (KISTI, Korea) 570,020 13.929 18 0.391
11 Oakforest-PACS 556,104 13.555 19 0.385
2020/6/23
TOP 500 Rank
System Accelerator Cores HPL Rmax
(Pflop/s)
Power (kW)
GFLOPS/
W 1 394 MN-3, Preferred Networks, Japan MN-Core 2,080 1.621 77 21.108 2 7 Selena, NVIDIA, USA NVIDIA A100 277,760 27.580 1,344 20.518
3 469 NA-1, PEZY, Japan PEZY-SC2 1,271,040 1.303 80 *18.433
4 205 A64FX Prototype, Fujitsu, Japan 36,864 1.999 118 16.876
5 27 AiMOS, USA NVIDIA V100 130,000 8.339 512 16.285
6 6 HPC5, Italy NVIDIA V100 669,760 35.450 2,252 15.740
7 422 Satori, USA NVIDIA V100 34,040 1.464 94 15.574
8 2 Summit, USA NVIDIA V100 2,414,592 148.6 10,096 14.719
9 1 Fugaku, Fujitsu, Japan 7,299,072 415.53 28,335 14.665
10 9 Marconi-100, Italy NVIDIA V100 347,776 21.64 1,476 14.671
(13)
Nov.’17
Reedbush-L, U.Tokyo, Japan NVIDIA P100 16,640 806 79 10.167
(19) Reedbush-H, U.Tokyo, Japan NVIDIA P100 17,760 802 94 8.576
2020/6/23
講義日程（工学部共通科目 ）
1. 4月14日 ガイダンス
2. 4月21日
l 並列数値処理の基本演算（座学）
3. 4月28日：スパコン利用開始
l ログイン作業、テストプログラム実行 4. 5月12日
l 高性能プログラミング技法の基礎１
（階層メモリ、ループアンローリン グ）
5. 5月19日
l 高性能プログラミング技法の基礎2
（キャッシュブロック化）
6. 5月26日
l 行列-ベクトル積の並列化
7. 6月2日
l べき乗法の並列化
8. 6月9日
l 行列-行列積の並列化（１）
9. 6月16日
l 行列－行列積の並列化（２）
10. 6月23日
l ＬＵ分解法（１）
l コンテスト課題発表
11. 6月30日
l ＬＵ分解法（２）
12. 7月7日
l ＬＵ分解法（３）、非同期通信
13. 7月14日
l RB-Hお試し、研究紹介他
LU 分解法（中級レベル以上 ）の演習日程
並列化が難しいので、３週間確保してあります。
1. 今週
• 講義（知識、アルゴリズムの理解）
• 並列化の検討
2. 来週
• LU分解法の逐次アルゴリズムの説明
• ＬＵ分解法の並列化実習（１）
3. 再来週
• LU分解法の並列化実習（２）
講義の流れ
1. LU分解法
• ガウス・ジョルダン法
• ガウス消去法
• 枢軸選択
• LU分解法
• 外積形式、内積形式、クラウト法、ブロック形式ガウス法、縦ブロックガウ ス法、前進・後退代入
2. サンプルプログラムの実行
3. 並列化のヒント
4. 実習課題
5. レポート課題
LU 分解法の概略
いろいろな変種があります
密行列に対する連立一次方程式
•
以下の式
ここで は実数の密行列 は 実数のベクトルとすると、解ベクトル を 求めること。
•
解ベクトルを求める方法は、以下の二種類が 知られている
1. 直接解法
行列操作により厳密解を求める方法
2. 反復解法
近似解を反復計算で解に収束させ求める方法
b Ax =
A x, b
x
ガウス・ジョルダン法
• 基本的な消去法により解を求める
• 第１ステップ
• 第２ステップ
• 最終ステップ
, ,
, 2
, 2 ,
2 2
, 22
1 1
2 12 1
11
n n
nn n
n
n n
n n
b x
a x
a
b x
a x
a
b x
a x
a x
a
= +
+
= +
+
= +
+ +
!
"
!
! 第一行をもとに
係数を消去
, , ,
,
, , 2 ,
, 2 2
, 22
, , 1 ,
, 1 1
11
n n
nn
n n
n n
b x
a
b x
a x
a
b x
a x
a
= +
+
= +
+
= +
+ +
!
"
!
!
第二行をもとに 係数を消去
*
*
* 2 2
, 22
* 1 1
11
n n
nn x b a
b x
a
b x
a
=
=
=
!
割り算のみで 解を得る
ガウス・ジョルダン法
• 右辺ｂの代わりに単位行列 Ｉ を用意し て同様の操作をすれば、最終ステップで は逆行列が求まる
• 各ステップでの計算量が同じなので、
並列化時の負荷バランスが良い
ガウス消去法
• 対角線より上の要素をゼロにしない方法
• 第１ステップ
• 第２ステップ
• 最終ステップ
, ,
, 2
, 2 ,
2 2
, 22
1 1
2 12 1
11
n n
nn n
n
n n
n n
b x
a x
a
b x
a x
a
b x
a x
a x
a
= +
+
= +
+
= +
+ +
!
"
!
!
第一行をもとに 係数を消去
, , ,
,
, 2 ,
2 2
, 22
1 1
2 12 1
11
n n
nn n n
n n
b x
a
b x
a x
a
b x
a x
a x
a
= +
+
= +
+
= +
+ +
!
"
!
!
第二行をもとに 係数を消去
この消去を
前進消去（forward elimination） とよぶ
*
*
, 2 ,
2 2
, 22
1 1
2 12 1
11
n n
nn n n
n n
b x
a
b x
a x
a
b x
a x
a x
a
=
= +
+
= +
+ +
!
"
#
#
ガウス消去法
• 前進消去後、最後の項から順に解を求めていく
この代入処理を、後退代入（backward substitution）とよぶ
*
*
, 2 ,
2 2
, 22
1 1
2 12
1 11
n n
nn
n n
n n
b x
a
b x
a x
a
b x
a x
a x
a
=
= +
+
= +
+ +
!
"
#
#
!
, /
) (
, /
1 ,
1 ,
1 1
1
*
*
+ - -
+ -
+ -
- = -
=
n n
n n
n n
nn n
n
a a
b x
a b
x
ガウス消去法
• ガウス消去法は、ガウス・ジョルダン法に比べ、
消去演算をする範囲が少ない
（基本行より下のみ）
• 演算量が低下する：
• 基本行より下のみ演算するため、並列化するとガ ウス・ジョルダン法に比べて、負荷バランスの
劣化を起こしやすい
• 並列処理に向かないと考えた専門家がいた。
• 現在はデータ分散の改良や通信の隠蔽技法、
ハードウエア能力向上から、ガウス消去法のほうが 高速である。
3 3
( 2 / 3 ) n
n ®
ピボッティング
• ガウス・ジョルダン法、ガウス消去法とも、基本行の係数がゼ ロだと、ゼロによる除算が生じ、計算が続行できない
• これを回避するため、消去する列から最も係数の大きなもの を選択して、基本行と入れ替える
（枢軸選択、ピボッティング、pivot selection）
, ,
, 2
, 2 ,
2 2
, 22
1 1
2 12
1 11
n n
nn n
n
n n
n n
b x
a x
a
b x
a x
a
b x
a x
a x
a
= +
+
= +
+
= +
+ +
!
"
!
!
第１行をもとに 係数を消去
０
ピボッティング
• ピボッティングには以下の２種の方法がある
1. 完全ピボッティング
更新対象全体から最大のものを選ぶ方法
2. 部分ピボッティング
更新対象の列または行から最大のものを選ぶ方式
• ピボッティングの手間、経験的な数値安定性から 部分ピボッティングが用いられることが多い
n n
nn n
n n
n n
n n
b x
a x
a x
a
b x
a x
a x
a
b x
a x
a x
a
= +
+ +
= +
+ +
= +
+ +
!
"
!
!
2 1
1
2 2
2 22
1 21
1 1
2 12
1 11
LU 分解法
• ガウス消去法のような消去処理を行列演算として定式化
• 連立一次方程式の行列表記：
n n
nn n
n n
n n
n n
b x
a x
a x
a
b x
a x
a x
a
b x
a x
a x
a
= +
+ +
= +
+ +
= +
+ +
!
"
!
!
2 1
1
2 2
2 22 1
21
1 1
2 12 1
11
úú úú û ù êê
êê ë é
= úú
úú û ù êê
êê ë é
= úú
úú û ù êê
êê ë é
=
n nn n
n n
n n
b b b b
x x x x
a a
a
a a
a
a a
a
A ! !
"
!
"
"
2 1 2
1
2 1
2 22
21
1 12
11
, ,
b x
A =
LU 分解法
• LU分解法では、以下の３つのステップで解を計算する
• 第１ステップ：行列AのLU分解
• 第２ステップ：前進代入
• 第３ステップ：後退代入
úú úú û ù êê
êê ë é
= úú
úú û ù êê
êê ë é
=
nn n
nn n
n u
u
u u
u U
l l
l l l l
L " !
#
#
"
!
22
1 12
11
2 1
22 21 11
,
,LU A =
b Ux
L
b x
LU b Ax
=
=
= ) (
, )
(
,
Ux c
b Lc
=
= , Lc = b
úú úú û ù êê
êê ë é
= úú úú û ù êê
êê ë é úú úú û ù êê
êê ë é
n n
nn n
n b
b b
c c c
l l
l
l l
l
!
!
"
#
!
2 1 2
1
2 1
22 21
11
：ベクトルｃ を求める
c Ux =
úú úú û ù êê
êê ë é
= úú úú û ù êê
êê ë é úú úú û ù êê
êê ë é
n n
nn n
c c c
x x x
u u
u u
u
!
!
!
"
#
2 1 2
1 22
1 12
11
：解ベクトルxを求める
LU 分解法
• 行列
A
のLU
分解 には、データアクセス の違いから以下の３種の方法が知られている1. 外積形式ガウス法（
outer-product form
）• 普通の消去法から導出
2. 内積形式ガウス法（
inner-product form
）• LU分解がなされたとして、Lの対角要素を１に 固定して導出
3. クラウト法（
Crout method
）• LU分解がなされたとして、Uの対角要素を１に 固定して導出
LU
A =
LU 分解法の種類
•
外積形式（ outer-product form ）ガウス法
• ガウス消去法と同等の操作で
LU
分解する• 第
k
列を消去したい場合、係数 を用いて を消去
n n
nn k
nk
k n
kn k
kk
n n
n n
b x
a x
a
b x
a x
a
b x
a x
a
b x
a x
a x
a
= +
+
= +
+
= +
+
= +
+ +
!
!
!
!
!
!
!
!
2 2
2 22
1 1
2 12
1 11
a
kka
k,k+1, a
k,k+2, ! , a
k,n外積形式ガウス法
•
すなわち列の消去は、
•
これを行列表記にすると、行列 L を
とすると、この消去は
n k
k i
a a
a
a
ik-
kk(
ik/
kk), = + 1 , + 2 ,...,
,
1 1
1 1
, 1
úú úú úú úú
û ù
êê êê êê êê
ë é
=
+
mk k k k
l L l
!
"
!
+1
=
kk
k
A U
L
n k
i
a a
l
ik ik kk,...,
1
), /
( +
=
-
=
外積形式ガウス法
•
一般的に
•
したがって LU 分解は
•
ここで、 は の要素の符号を反転させ たものであり、容易に得られる
• 消去作業が終われば行列Lが得られる
U A
L L
L
L
n-1 n-2!
2 1=
LU
U L
L L
L
U L
L L
L A
n n
n n
=
=
=
- - -
- -
-
- - -
) (
) (
1 1 1
2 1
2 1
1
1 1
2 2
1
!
!
-1
L
kL
kfor (k=0; k<n; k++) { dtmp = １.0 / A[k][k];
for (i=k+１; i<n; i++) { A[i][k] = A[i][k]*dtmp;
}
for (j=k+１; j<n; j++) { dakj = A[k][j];
for (i=k+１; i<n; i++) {
A[i][j] = A[i][j]–A[i][k]*dakj;
} }
L
U 注意：
Lの対角要素は １であることを仮定
（計算しない）
→Uの対角要素を 入れる
A
更新 k
k
参照
do k=1, n
dtmp = １.0d0 / A(k, k) do i=k+１, n
A(i, k) = A(i, k) * dtmp enddo
do j=k+１, n dakj = A(k, j) do i=k+１, n
A(i, j) = A(i, j)–A(i, k)*dakj enddo
enddo enddo
L
U 注意：
Lの対角要素は １であることを仮定
（計算しない）
→Uの対角要素を 入れる
A
更新 k
k
参照
外積形式ガウス法のまとめ
• 外積形式ガウス法では分解列の右側の 領域が更新される
•
right-looking アルゴリズムと呼ばれる
• 外積形式ガウス法は並列化に向く
•
処理の中心の更新領域が多い
• 負荷バランスよくデータ分散できる
•
更新処理が、分解行と分解列という少ない
データを所有するだけで、要素ごとに独立
して行える
内積形式ガウス法
• 内積形式（inner-product form）ガウス法
• LU分解がなされたと仮定した上で、行列Lの対角要素を １として導出した方法
úú úú û ù êê
êê ë é úú úú û ù êê
êê ë é
= úú úú û ù êê
êê ë é
nn n
n nn n
n n
n n
u u
u u
u
l l l a
a a
a a
a
a a
a
!
"
#
#
"
!
#
#
#
#
#
0 1
0 1
1
22
1 12
11
2 1
21
2 1
2 22
21
1 12
11
1 11
1 31
11 31
21 11
21
11 11
,...., ,
,
n
n
u a
l a
u l
a u
l
u a
=
=
=
= u
11が求まるl
21が求まる内積形式ガウス法
•
この導出作業を一般化すると、以下の二部分 に分かれる
• （I） uの導出部
• （II） （I）で得られた値を元に、L の導出部
•
まとめると
•
（ I ）
•
（ II ）
å
-=
= -
=
=
1 1 1
1 1
) ,...,
3 , 2 (
,
i j
jk ij
k ik
k k
k i
u l a
u
a u
å
-=
+ +
= -
= 1
1
) ,..., 2
, 1 (
, /
) (
k j
kk jk
ij ik
ik a l u u i k k n
l
for (k=0; k<n; k++) { for (j=0; j<k; j++) {
dajk = A[j][k];
for (i=j+１; i<n; i++) {
A[i][k]= A[i][k] –A[i][j]*dajk;
} }
A[k][k]=1.0 / A[k][k];
for (i=k+１; k<n; k++) { A[i][k]=A[i][k]*A[k][k];
} }
L
U
A
k
k
参照 更新
更新と参照
do k=1, n do j=1, k
dajk = A(j, k) do i=j+１, n
A(i, k)= A(i, k) –A(i, j) * dajk;
enddo enddo
A(k, k) =1.0d0 / A(k, k) do i=k+１, n
A(i, k)=A(i, k) * A(k, k) enddo
enddo
L
U
A
k
k
参照 更新
更新と参照
内積形式ガウス法のまとめ
•
内積形式ガウス法では、分解列の左側の領 域が主に参照される
•
left-looking
アルゴリズムと呼ばれる•
内積形式ガウス法の並列化
• 行列
A
を列方向分散（＊，Cyclic
）• 参照領域のデータがないので、通信多発
（ベクトルリダクションが毎回必要）
• 行列
A
を行方向分散（Cyclic
，＊）• 上三角行列Uの要素（データ数が少ない）を所有 すれば、独立して計算可能
クラウト法
• クラウト法（
Clout Method
）• LU分解がなされたと仮定した上で、行列Uの対角要 素を１として導出した方法（ｃｆ．内積形式ガウス法）
úú úú û ù êê
êê ë é úú úú û ù êê
êê ë é
= úú úú û ù êê
êê ë é
1 0
1 1
0
1 12
2 1
22 21
11
2 1
2 22
21
1 12
11
!
"
#
#
"
!
#
#
#
#
# n
nn n
nn n n
n
n
n u u
l l
l
l l
l
a a
a
a a
a
a a
a
n n
n n
a u
l a
u l
a u
l
a l
a l
a l
1 1
11 13
13 11
12 12
11
1 1
21 21
11 11
,...., ,
, ,
=
=
=
=
=
=
12 が求まる
u
ｌの第1列が 求まる
クラウト法
•
この計算を一般化すると、
•
L の第 k 列を求める場合
•
U の第 k 行を求める場合
å
-=
+
= -
=
11
) ,...,
1 ,
( ,
k j
jk ij
ik
ik
a l u i k k n
l
å
-=
+
= -
=
11
) ,...,
1 ,
( , /
)
(
ki
kk ij
ki kj
kj
a l u l j k k n
u
クラウト法（Ｃ言語）
A[0][0]=1.0/A[0][0];
for (j=１; j<n; j++) {
A[0][j]=A[0][j]*A[0][0]; } for (k=0; k<n; k++) {
for (j=0; j<k; j++) { dajk=A[j][k];
for (i=k; i<n; i++) {
A[i][k]=A[i][k]-A[i][j]*dajk;
} }
A[k][k]=1.0/A[k][k];
for (i=0; i<k; i++) { daki=A[k][i];
for (j=k+１; j<n; j++) {
A[k][j]=A[k][j]-daki*A[i][j];
for (j=k+１; j<n; j++) {} } A[k][j]=A[k][j]*A[k][k]; } }
L
U
A
k
k
参照
更新 参照
更新
クラウト法（ Fortran 言語）
A(1,1)=1.0d0/A(1,1) do j=2, n
A(1, j) =A(1, j) * A(1, 1) enddo do k=1, n
do j=1, k dajk=A(j, k) do i=k, n
A(i, k)=A(i, k) - A(i, j) * dajk enddo; enddo
A(k, k) =1.0d0 / A(k, k) do i=1, k
daki=A(k, i) do j=k+１, n
A(k, j)=A(k, j) – daki * A(i, j) enddo; enddo
do j=k+１, n
A(k, j)=A(k, j) * A(k, k) enddo enddo
L
U
A
k
k
参照
更新 参照
更新
クラウト法
• クラウト法では、最内ループの交換ができる
• 長さ（１～ｋ－１）のループ、長さ（ｋ－ｎ）の
ループの内、最も長いループを最内に移動可
• ベクトル計算機で実行性能が良い
• 分解列および分解行の外側に２つの参照領域
• 分散メモリ型並列計算機での実装が困難
∵どのようにデータ分割しても大量通信発生
• 共有メモリ型並列計算機では並列化可能
∵
参照領域があれば分解列と分解行は独立 に計算可能• 行列Aを小行列に分解し、その小行列単位でLU分解す る方法。LU分解と行列-行列積で実現できる。
• 具体的には （各小行列を各PEが所有）
とすると、右辺は
úú ú û ù êê
ê ë é úú ú û ù êê
ê ë é
= úú ú û ù êê
ê ë é
33 23 22
13 12
11
33 32
31
22 21
11
33 32
31
23 22
21
13 12
11
0 ~
~
~
~
~
~
~
~
~
~
~~ 0
~
~
~
~
~
~
~
~
~
U U U
U U
U L
L L
L L
L A
A A
A A
A
A A
A
33 33 23
32 13
31 33
22 32 12
31 32
11 31 31
23 22 13
21 23
22 22 12
21 22
11 21 21
13 11 13
12 11 12
11 11 11
~
~
~
~
~
~ ~
~ ,
~
~
~ ~
~ ,
~ ~ ~ ~ ~ ~ ~ ,
~ ,
~
~
~ ~
~ ,
~ ~ ~ ~ ~ ,
~ ,
~ , ~
~
~
~
U L U
L U
L A
U L U
L A
U L A
U L U
L A
U L U
L A
U L A
U L A
U L A
U L A
+ +
= +
=
=
+
= +
=
=
=
=
=
