LAPACK・BLASを用いた連立1次方程式の数値計算
理化学研究所計算科学研究センター 特別研究員
寺尾 剛史

2

今回の目標
連立1次方程式の数値解法の実装（密行列）LU分解
残差反復改良を用いた混合精度数値計算法
BLAS・LAPACKを用いたコードの作成に慣れる
密行列に対する高性能なコードの作成

3

連立1次方程式の数値計算
 正則な𝑛 × 𝑛行列𝐴、 𝑛次ベクトル𝑏に対して𝐴𝑥 = 𝑏を満たす、ベクトル𝑥を求めたい。
 数値計算法は、計算速度と精度が求められる。
 計算速度はわかりやすい（実行時間を計ればよい）。
 精度とは、近似解𝑥に対して
ො‖𝑥ො − 𝑥‖がどれだけ小さいかである。（どう計算する？）
 精度保証付き数値計算
 残差ベクトルのノルム‖𝑏 − 𝐴𝑥‖が代用されることが多い。
ො4

連立1次方程式について
 連立1次方程式は、様々な分野に応用される。
 密行列：
LU分解、コレスキー分解、QR分解
 疎行列：
定常反復法：Jacobi法、Gauss-Seidel法、SOR法
非定常反復法：CG法、ICCG法、Bi-CG法、GMRES法
 問題の計算コストの大部分を占める場合がある。
 問題や係数行列で有効な手法が異なる。

5LU分解 𝑛 × 𝑛の正則な行列𝐴に対するLU分解は、上三角行列𝐿と下三角行列𝑈を用いて𝐴=𝐿×𝑈2 3の形に分解する。（計算回数： 𝑛 flops）
3

 以上のように分解できたとき、連立1次方程式𝐴𝑥 = 𝑏は𝑥 = 𝑈 −1 (𝐿−1 𝑏)で計算できる。 （計算回数：𝑂(𝑛2 ) flops）
 LU分解を用いた連立1次方程式の解法では、LU分解に多くのコストが掛かる。

6LU分解（軸交換無し）
 𝑛 = 3の場合

7

LU分解の例（軸交換無し）
10 0567𝐴 = 10 20 23 = 𝑙21 1 0𝑙31 𝑙32 115 50 67
5 671 0 00 𝑢22 𝑢23= 2 1 03 𝑙32 10 0 𝑢33
1= 231= 230140140010015 6 70 8 90 0 𝑢33
5 6 70 8 90 0 10𝑢1100𝑢12𝑢220𝑢13𝑢23𝑢33

𝑢22 = 20 − 2 ∗ 6 = 𝟖
𝑢23 = 23 − 2 ∗ 7 = 𝟗50 − 3 ∗ 6𝑙32 ==𝟒𝑢22
𝑢33 = 67 − 3 ∗ 7 − 4 ∗ 9 = 𝟏𝟎

8

プロクラムの高速化
 高速な数値計算を行うにはどうする？
逐次実行の高速化• 遅いメモリへのアクセスの削減• キャッシュ・レジスタブロッキング

並列計算• 計算の並列化• OpenMP、MPI 個人で行うには限界がある（手間もかかる）
 HPCのプロが作った数値計算ライブラリを用いることも重要！

9

BLAS・LAPACKとは BLASとは、Basic Linear Algebra Subprogramsの略
 行列やベクトルの基本的な計算を行う関数群
 行列積やベクトルのノルム計算など

 LAPACKとは、Linear Algebra PACKageの略
 連立1次方程式や固有値問題など

 IEEE754の単精度、倍精度からなる実数、複素数を要素とするベクトルや行列をサポート
 密行列に対するコードの高速化や簡略化ができる。

10void sgesv_(N, NRHS, A, LDA, IPIV, B, LDB,INFO);

 単精度でLU分解を用いて𝐵 ← 𝐴−1 𝐵を解く。
 N(入力) – INTEGER
 行列Aの行列サイズ

 NRHS（入力）ーINTEGER 右辺ベクトルBの数 A（入力・出力）ーSINGLE PRECISION （入力）係数行列
 （出力）LU分解の計算結果

 LDA（入力）ーINTEGER 行列Aの最初の次元数。max 𝑛, 1 である必要がある。

 IPIV（出力）ーINTEGER LU分解での軸交換の情報 B（入力・出力）ーSINGLE PRECISION （入力）右辺ベクトル
 （出力） 𝐴−1 𝐵の数値解
 LDB（入力）ーINTEGER
 行列Bの最初の次元数。max 𝑛, 1 である
必要がある。
 INFO（出力）ーINTEGER ＝０：正常 ＜０：INFO=-iのときi行目の引数がおかしい ＞０：INFO=iならばU(i,i)=0で計算できない11

BLAS・LAPACKの注意点
 引数はすべてポインタで渡す。
 Fortranで記述されているため、C言語では二次元配列の順番に注意する
 一次元配列にすれば問題ない

C： A[2][3]
1 2 3
4 5 6

Fortran：
1 3 5
2 4 6

 今回のサンプルプログラムは、１次元配列を用いているので、気にする必要はない。
 コンパイラオプションを追加する必要がある。
 OpenBLAS: -lopenblas, LAPACK;: -llapack , MKL: -mkl など

12

サンプルプログラム
 サンプルを以下からコピーしてください。
/home/guest32/LinSys 次を実行してコピー：>> cp –r ~guest32/LinSys ./LinSys>> cd LinSys>> ls
ex01.cex02.crun.sh 次でコンパイル：>> icc –mkl –qopenmp –O3 ex01.c
 次で実行：>> pjsub run.sh

13

課題１
1. ex01.cを確認してください
2. 様々な行列サイズで残差と計算性能の比較を行ってください
3. 様々なスレッド数で実行してください 実行結果（例）：(Double) N: 1000, Residual: 2.444267e-12, Time: 1.453111e-01(Single) N: 1000, Residual: 7.654402e-04, Time: 1.061270e-01

混合精度計算を用いた連立1次方程式の解法
1415

混合精度数値計算
複数の精度を用いた数値計算方法
例えば、（Single）＋（Double）で（Double）と同程度の精度で数値計算

主な目的は：
計算性能の向上
省エネ化

今回はLU分解を用いた混合精度計算法を紹介する

16

連立1次方程式に対する反復改良法
 連立1次方程式𝐴𝑥 = 𝑏の数値解を𝑥とする
ො また、𝑥の誤差をΔ𝑥ො≔ 𝑥 − 𝑥とする
ො
 このとき、𝐴 𝑥ො + Δ𝑥 = 𝑏Δ𝑥 = 𝐴−1 𝑏 − 𝐴𝑥ොが成り立つ。
 ここで、𝑏 − 𝐴𝑥には高精度な行列ベクトル積が必要ො LU分解の計算結果は1回目で計算出来ているため、計算量は𝑂(𝑛2 ) flopsで計算できる。（LU分解は2/3𝑛3 flops）

17

混合精度計算を用いた連立1次方程式
 LU分解を用いた混合精度計算のイメージ
計算時間

LU分解法（double）

LU分解法（single）

反復改良（double）

 残差反復が効率的に行われる場合、計算時間を削減できる。

18

混合精度の利点と欠点
 利点：
主要なコストであるLU分解のコストを低減できる
GPUなどの低精度演算が高速な計算機ではさらに有効

 欠点：
扱える条件数𝜅 𝐴 ≔ 𝜎𝑚𝑎𝑥 (𝐴)/𝜎𝑚𝑖𝑛 (𝐴)が変化
倍精度の場合は𝜅 𝐴 ≈ 1016 、単精度＋倍精度の場合は
𝜅 𝐴 ≈ 107 が限界（単精度のLU分解に依存）
行列やベクトルを保存するためのメモリが必要void dgemv_(TRANS, M, N, alpha, A, LDA, x,INCX, beta, y,INCY);

19 倍精度で𝑦 ← 𝛼𝐴𝑥 + 𝛽𝑦を計算。
 TRANS(入力)：CHAR ‘N’:そのまま ‘T’: 転置 M（入力）ーINTEGER
 op(A)の行数
 N（入力）ーINTEGER
 op(A)の列数

 alpha（入力）ーDOUBLE PRECISION スカラ値を指定

 A（入力）ーDOUBLE PRECISION 行列Aの配列

 LDA（入力）ーINTEGER
 行列Aの最初の次元数。max 𝑛, 1 である必要がある。

 x（入力）ーDOUBLE PRECISION ベクトルxの配列
 INCX（入力）ーINTEGER
 ベクトルxのインクリメント幅

 beta（入力） スカラ値を指定

 y（入力・出力）
 ベクトルyの配列

 INCY（入力）ーINTEGER
 ベクトルyのインクリメント幅20void dgetrs_(TRANS, N, NRHS, A,LDA,IPIV, B,LDB, INFO); LU分解の計算結果を用いて𝐵 ← (𝐿𝑈)−1 𝐵を解く。
 N(入力) – INTEGER
 行列Aの行列サイズ

 NRHS（入力）ーINTEGER 右辺ベクトルBの数 A（入力）ーDOUBLE PRECISION
 （出力）LU分解の計算結果

 LDA（入力）ーINTEGER
 行列Aの最初の次元数。max 𝑛, 1 である必要がある。

 IPIV（入力）ーINTEGER LU分解での軸交換の情報 B（入力・出力）ーDOUBLE PRECISION （入力）右辺ベクトル
 （出力） 𝐴−1 𝐵の数値解
 LDB（入力）ーINTEGER
 行列Aの最初の次元数。max 𝑛, 1 である必要がある。
 INFO（出力）ーINTEGER ＝０：正常 ＜０：INFO=-iのときi行目の引数がおかしい ＞０：INFO=iならばU(i,i)=0で計算できない21

void dcopy_(N, X, INCX, Y, INCY);
 倍精度配列XをYにコピー
 N(入力) – INTEGER
 ベクトルXのサイズ

 X（入力）ーDOUBLE PRECISION
 コピーする配列

 INCX（入力）ーINTEGER Xのインクリメント幅 Y（出力）ーDOUBLE PRECISION
 コピーされる配列

 INCY（入力）ーINTEGER Yのインクリメント幅22void daxpy_(N, alpha, X, INCX, Y, INCY); 𝑦 ← 𝛼𝑥 + 𝑦を計算
 N(入力) – INTEGER
 ベクトルXのサイズ

 alpha（入力） スカラ値を指定
 X（入力）ーDOUBLE PRECISION
 コピーする配列

 INCX（入力）ーINTEGER Xのインクリメント幅 Y（出力）ーDOUBLE PRECISION
 コピーされる配列

 INCY（入力）ーINTEGER Yのインクリメント幅

23アルゴリズム
 混合精度計算を用いた連立1次方程式を解くアルゴリズム

1. 𝑥ො ≈ 𝐴−1 𝑏 ：LU分解を用いて単精度で連立1次方程式を解く(sgesv)
2. 𝑟 ← 𝑏：配列のコピー（dcopy）
3. 𝑟 ← 𝑟 − 𝐴𝑥：行列ベクトル積を倍精度計算(dgemv)
ො
4. 𝑟 ← 𝐴−1 𝑟 ：LU分解の計算結果を用いて連立1次方程式を解く(dgetrsまたはsgetrs)
5. 𝑥ො ← 𝑥ො + 𝑟：数値解を更新する(daxpyまたはfor文)

6. 2～5を十分に繰り返す
 Aやbなどの情報を保存する必要があることに注意！

24

課題2
ex02.cを確認してください反復回数ごとの残差と計算性能の比較を行ってください倍精度や単精度で計算した場合と混合精度で計算した場合の比較を行ってください25

反復改良のコード例

単精度でLU分解法で解く
dlag2s_(&n, &n, A, &n, As, &n, &info);

反復改良
for (i = 0; i < 30; i++) {dlag2s_(&n, &one, b, &n, xs, &n, &info);

dcopy_(&n, b, &one, r, &one);// solve the linear systemdgemv_(&chrN, &n, &n, &dmone, A, &n, x, &one, &done, r, &one);sgesv_(&n, &one, As, &n, IPIV, xs, &n, &info);

rnrm = fabs(r[idamax_(&n, r, &one) - 1]);

slag2d_(&n, &one, xs, &n, x, &n, &info);

xnrm = fabs(x[idamax_(&n, x, &one) - 1]);
if (rnrm < sqrt(n) * EPS * anrm * xnrm) {printf("No. of iter. : %d¥n", i + 1);
break;}

dlag2s_(&n, &one, r, &n, xs, &n, &info);
sgetrs_(&chrN, &n, &one, As, &n, IPIV, xs, &n, &info);slag2d_(&n, &one, xs, &n, r, &n, &info);
daxpy_(&n, &done, r, &one, x, &one);}

数値実験の例

26

Single：3.11E-02

残差の推移（n=10000）
1.00E+00

1.00E-01

Double：9.47E-11

3.11E-02

1.00E-02

Mix：4.46E-11

1.00E-03
2.66E-05

1.00E-04

4回の反復で倍精度と同等の精度で計算1.00E-05
1.00E-06
5.67E-08

1.00E-07
1.00E-08

3.43E-10

1.00E-09
1.00E-10

4.46E-11

4.48E-1145

1.00E-110

123

反復回数

行列サイズ、条件数で変動