計算機実験 (講義 L3)
計算機実験 (講義 L3)
藤堂眞治
2015-05-13/14
1 連立一次方程式の直接解法
2 反復法
3 C 言語における行列の取り扱い
4 LAPACK ライブラリ
1 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
連立一次方程式の現れる例
微分方程式の初期値問題の陰解法
非線形連立方程式に対するニュートン法
x′ = x −
( ∂f(x) )−1
∂x
f(x)
偏微分方程式の境界値問題の差分法による求解
ベクトルに逆行列を掛ける代わりに連立一次方程式を解く場
合が多い
x = A−1 b ⇒ Ax = b
2 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ポアソン方程式の境界値問題
二次元ポアソン方程式
∂ 2 u(x, y ) ∂ 2 u(x, y )
+
= f (x, y )
∂x 2
∂y 2
0 ≤ x ≤ 1, 0 ≤ y ≤ 1
ディリクレ型境界条件: u(x, y ) = g (x, y ) on ∂Ω
有限差分法により離散化
▶
▶
▶
▶
x 方向、y 方向をそれぞれ n 等分: (xi , yj ) = (i/n, j/n)
(n + 1)2 個の格子点の上で u(xi , yj ) = uij が定義される
そのうち 4n 個の値は境界条件で定まる
ポアソン方程式を中心差分で近似 (h = 1/n)
ui,j+1 − 2uij + ui,j−1
ui+1,j − 2uij + ui−1,j
+
= fij
h2
h2
残り (n − 1)2 個の未知数に対する連立一次方程式
3 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ポアソン方程式の境界値問題
ノイマン型境界条件の場合
▶
▶
境界上で u(x, y ) の微分が定義される。
例) ∂u(0, y )/∂x = h(0, y )
境界条件を差分近似で表す
u1j − u0j
= h0j
h
j = 1 · · · (n − 1)
(n + 1)2 − 4 個の未知数に対して、ポアソン方程式の差分近
似とあわせて、合計 (n − 1)2 + 4(n − 1) = (n + 1)2 − 4 個の
連立一次方程式
4 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
逆行列の「間違った」求め方
線形代数の教科書に載っている公式
A−1 =
Ã
|A|
|A|: A の行列式、Ã: A の余因子行列
n × n 行列の行列式を定義通り計算すると、計算量∼O(n!)
したがって、上の方法で逆行列を計算すると、計算量∼O(n!)
n = 100 の場合: n! ≈ 9.3 × 10157
5 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
逆行列の「正しい」求め方
連立一次方程式 Ax = ej を全ての ej について解く
Gauss の消去法による連立一次方程式の解法: 計算量∼O(n3 )
Gauss の消去法の途中で出てくる下三角行列 (L) と上三角行
列 (U) 行列を再利用 (LU 分解) すれば、逆行列全体を求める
ための計算量も O(n3 )
行列式も O(n3 ) で計算可
n = 100 の場合: n3 = 106 ≪ 9.3 × 10157
6 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法
解くべき連立方程式
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
(1)
an1 x1
(1)
an2 x2
(1)
an3 x3
(1)
ann xn
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a21 x1 + a22 x2 + a23 x3 + · · · + a2n xn = b2
a31 x1 + a32 x2 + a33 x3 + · · · + a3n xn = b3
···
+
+
+ ··· +
(1)
= bn
ある行を定数倍しても、方程式の解は変わらない
ある行の定数倍を他の行から引いても、方程式の解は変わら
ない
7 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法
(1)
(1)
1 行目を mi1 = ai1 /a11 倍して、i 行目 (i ≥ 2) から引く
(1)
(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
(2)
(2)
(2)
(2)
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a22 x2 + a23 x3 + · · · + a2n xn = b2
a32 x2 + a33 x3 + · · · + a3n xn = b3
···
(2)
(2)
(2)
(2)
an2 x2 + an3 x3 + · · · + ann xn = bn
ここで
(2)
(1)
(1)
i ≥ 2, j ≥ 2
(1)
i ≥2
aij = aij − mi1 a1j
(2)
bi
(1)
= bi
− mi1 b1
8 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法
(2)
(2)
2 行目を mi2 = ai2 /a22 倍して、i 行目 (i ≥ 3) から引く
(1)
(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
(3)
(3)
(3)
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a22 x2 + a23 x3 + · · · + a2n xn = b2
a33 x3 + · · · + a3n xn = b3
···
(3)
(3)
(3)
an3 x3 + · · · + ann xn = bn
ここで
(3)
(2)
(2)
i ≥ 3, j ≥ 3
(2)
i ≥3
aij = aij − mi2 a2j
(3)
bi
(2)
= bi
− mi2 b2
9 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法
最終的には、左辺が右上三角形をした連立方程式となる
(1)
(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
(3)
(3)
(3)
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a22 x2 + a23 x3 + · · · + a2n xn = b2
a33 x3 + · · · + a3n xn = b3
···
(n−1)
an−1,n−1 xn−1
+
(n−1)
an−1,n xn
(n)
ann xn
(n−1)
= bn−1
(n)
= bn
これを下から順に解いていけばよい (後退代入)
10 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
練習問題
次の連立方程式をガウスの消去法で (手で) 解け
x1
18
1 4 7
2 5 8 x2 = 24
3 6 10
x3
31
11 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法のコード
for ( k = 0; k < n ; ++ k ) {
for ( i = k + 1; i < n ; ++ i ) {
for ( j = k + 1; j < n ; ++ j ) {
a [ i ][ j ] -= a [ k ][ j ] * a [ i ][ k ] / a [ k ][ k ];
}
b [ i ] -= b [ k ] * a [ i ][ k ] / a [ k ][ k ];
}
}
for ( k = n -1; k >= 0; --k ) {
for ( j = k + 1; j < n ; ++ j ) {
b [ k ] -= a [ k ][ j ] * b [ j ];
}
b [ k ] /= a [ k ][ k ];
}
C 言語では配列の添字が 0 から始まることに注意
12 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ピボット選択
(k)
ガウスの消去法の途中で akk が零になると、計算を先に進め
ることができなくなる
行を入れ替えても、方程式の解は変わらない ⇒ k 行以降で、
(k)
aik が非零の行と入れ替える (ピボット選択)
(k)
実際のコードでは、情報落ちを防ぐため、akk が零でない場
(k)
合でも、aik の絶対値が最大の行と入れ替える
ピボット選択が必要となる例
8
1 2 3
x1
3 6 4 x2 = 19
x3
23
4 6 7
13 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
ガウスの消去法の行列表示
(k)
(k)
akk を用いた aik (i > k) の消去は、方程式の両辺から
1
0 1
.
0 0 . .
. .
. .
1
. .
Mk =
.. ..
−mk+1,k
. .
. .
.. ..
−mk+2,k
. .
..
.. ..
.
0 0 ...
−mnk
1
..
0
..
.
.
1
0 ... 0 1
を掛けるのと等価: Mk A(k) = A(k+1) 、Mk b(k) = b(k+1)
14 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
LU 分解
Mk の逆行列
Lk = Mk−1
1
0 1
.
0 0 . .
. .
. .
1
. .
=
.. ..
mk+1,k
. .
. .
.. ..
mk+2,k
. .
..
.. ..
.
0 0 ...
mnk
1
..
0
..
.
.
1
0 ... 0 1
から L = L1 L2 · · · Ln を定義すると、L は下三角行列、また
R = A(n) (上三角行列) とすると、A = LU
15 / 31
計算機実験 (講義 L3)
連立一次方程式の直接解法
LU 分解
LU 分解による連立一次方程式の解法
▶
▶
▶
方程式は Ax = LUx = b と書ける
まず、Ly = b を解いて、y を求める (前進代入)
次に、Ux = y を解いて、x を求める (後退代入)
計算量はガウスの消去法と変わらない
一度 LU 分解をしておけば、異なる b に対する解も簡単に求
められる
行列式は U の対角成分の積で与えられる
16 / 31
計算機実験 (講義 L3)
反復法
直接法と反復法
直接法: 連立方程式を有限回数 (∼ n3 ) の手間で直接解く
反復法: Ax = b を、等価な x = ϕ(x) = Mx + c の形に変形
し、適当な初期値 x0 から出発して、x(k+1) = ϕ(x(k) ) を繰り
返して解を求める
▶
▶
▶
欠点: 有限回数では終わらない (あらかじめ定めた収束条件が
満たされるまで反復)
利点: 行列ベクトル積 Mx が計算できさえすればよい。特に
M が疎行列の場合には、Mx は非常に高速に計算できる可能
性がある。メモリの点でも有利
利点: 直接法に比べて、プログラムも比較的単純になる場合
が多い
17 / 31
計算機実験 (講義 L3)
反復法
反復法
行列 A を対角行列 D 、左下三角行列 E 、右上三角行列 F の和
に分解
Ax = (D + E + F )x = b
ヤコビ法: 対角成分以外を右辺に移す
x(k+1) = D −1 (b − (E + F )x(k) ) = −D −1 (E + F )x(k) + D −1 b
ガウスザイデル法: ヤコビ法で右辺の x の値として、各段階
ですでに得られている最新のものを使う
x(k+1) = D −1 (b − E x(k+1) − F x(k) )
x(k+1) = −(D + E )−1 F x(k) + (D + E )−1 b
18 / 31
計算機実験 (講義 L3)
反復法
反復法
SOR (Successive Over-Relaxation) 法: ガウスザイデル法にお
ける修正量に 1 より大きな値 (ω) を掛け、補正を加速
ξ (k+1) = D −1 (b − E x(k+1) − F x(k) )
x(k+1) = x(k) + ω(ξ (k+1) − x(k) )
ξ (k+1) を消去すると
x(k+1) =(I + ωD −1 E )−1 {(1 − ω)I − ωD −1 F }x(k)
+ ω(D + ωE )−1 b
反復法は常に収束するとは限らない
行列 A が対角優位、あるいは正定値対称行列の場合には収束
が保証される
19 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
一次元配列
(静的) 一次元配列 (ハンドブック 3.3.1 節)
double v [10];
v [0] = 1.0;
v [1] = 2.0;
...
要素数はコンパイル時にすでに決まっている定数でなければ
ならない
(動的) 一次元配列 (ハンドブック 3.11 節)
double * v ; /* ポ イ ン タ */
v = ( double *) malloc (( size_t )(10 * sizeof ( double ));
...
free ( v ); /* 確 保 し た 領 域 を 開 放 */
実行時に要素数を指定可能
20 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
ポインタと一次元配列
一次元配列を表す変数は、(実は) 最初の要素を指すポインタ
(ハンドブック 3.5.3 節)
▶
▶
▶
▶
▶
v と &v[0] は等価
(v+2) と &v[2] は等価
*v と v[0] は等価
*(v+2) と v[2] は等価
(v+2)[3] は?
C 言語では配列の添字は 0 から始まることに注意
double v[10]; と宣言した場合、v[0] ∼ v[9] の 10 個の要
素を持つ配列が作られる。v[10] は存在しない。値を代入し
たり参照しようとするとエラーとなる
21 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
二次元配列
C 言語では、二次元配列は一次元配列の先頭をさす (ポイン
タ) の配列として表される (と理解しておけば良い)
m[i] は、要素 m[i][0] を指すポインタ
▶
▶
▶
▶
▶
▶
▶
m と &m[0] は等価 (&m[0][0] ではない)
m[0] と &m[0][0] は等価
m[2] と &m[2][0] は等価
(m+2) と &m[2] は等価
(*(m+2))[3] と *(*(m+2)+3) と m[2][3] は等価
*(m+2)[3] と *((m+2)[3]) と *(m[5]) と m[5][0] は等価
[] は*よりも強い
22 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
動的二次元配列の確保
各行を表す配列とそれぞれの先頭アドレスを保持する配列の
二種類が必要
double ** a ;
m = 10;
n = 10;
a = ( double **) malloc (( size_t )( m * sizeof ( double *));
for ( int i = 0; i < m ; ++ i )
a [ i ] = ( double *) malloc (( size_t )( n * sizeof ( double ));
各行を保持する配列が、メモリ上で連続に確保される保証は
ない
LAPACK を使うときに問題となる
23 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
動的二次元配列の確保
二次元配列の要素を格納する長い配列を用意する
double ** a ;
m = 10;
n = 10;
a = ( double **) malloc (( size_t )( m * sizeof ( double *));
a [0] = ( double *) malloc (( size_t )( m * n * sizeof ( double ));
for ( int i = 1; i < m ; ++ i )
a [ i ] = a [i -1] + n ;
開放は逆の順序で行う
free ( a [0]);
free ( a );
24 / 31
計算機実験 (講義 L3)
C 言語における行列の取り扱い
動的二次元配列の確保
二次元配列の確保 (alloc dmatrix)、開放 (free dmatrix)、
出力 (print dmatrix)、読み込み (read dmatrix) を行うた
めのユーティリティ関数を準備
実習用ワークステーションの
/home/public/ce2015/ex3/matrix util.h に置いてある
使用例
# include " matrix_util . h "
...
double ** mat ;
mat = alloc_dmatrix (m , n );
...
free_dmatrix ( mat );
25 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
LAPACK (Linear Algebra PACKage)
線形計算のための高品質な数値計算ライブラリ
▶
▶
▶
▶
http://www.netlib.org/lapack
線形方程式、固有値問題、特異値問題、線形最小二乗問題など
(FFT 高速フーリエ変換は入っていない)
LAPACK 自体は Fortran で書かれている
ほぼ全ての PC、ワークステーション、スーパーコンピュー
タで利用可 (インストール済)
Netlib でソースが公開されているリファレンス実装は遅いが、
それぞれのベンダー (Intel、Fujitsu、etc) による最適化された
LAPACK が用意されている場合が多い (MKL、SSL2、etc)
LAPACK を使うことにより、高速で信頼性が高く、ポータブ
ルなコードを書くことが可能になる
26 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
LAPACK による連立一次方程式の求解
LU 分解を行うサブルーチン dgetrf
http://www.netlib.org/lapack/explore-html/d3/d6a/
dgetrf_8f.html
Fortran による関数宣言
subroutine dgetrf ( integer M , integer N ,
double precision , dimension ( lda , *) A ,
integer LDA , integer , dimension (*) IPIV ,
integer INFO )
A: 左辺の行列、N,M: 次元、IPIV: 選択されたピボット行のリ
スト、lda: M と同じで良い
27 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
LAPACK による連立一次方程式の求解
C 言語から呼び出すための関数宣言を作成 (ハンドブック
3.6.4 節)
void dgetrf_ ( int *M , int *N , double *A ,
int * LDA , int * IPIV , int * INFO );
関数名は全て小文字。関数名の最後に
(下線) を付ける
LU 分解の例
M = 10;
N = 10;
LDA = 10;
dgetrf_ (& M , &N , & A [0][0] , & LDA , & IPIV [0] , & INFO );
完全なソースコードは、実習用ワークステーションの
/home/public/ce2015/ex3/lu decomp.c に置いてある
28 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
C から Fortran のライブラリを呼び出す際の注意事項
スカラーも配列も全てポインタ渡しとする
C 言語は添字が 0 から始まる。Fortran は 1 から始まる
C と Fortran で、二次元配列のメモリ上での並びが違う
C は row-major: a[0][0], a[0][1], a[0][2], · · ·
Fortran は column-major: a(1,1), a(2,1), a(3,1), · · ·
C で作成した行列を Fortran に渡すと転置されてしまう
コンパイル時には-llapack オプションを指定し、LAPACK
ライブラリをリンクする (ハンドブック 3.1.6 節)
29 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
計算機を使う上で大事なこと
数学公式と数値計算アルゴリズムは別物
逆行列を明示的に作らない
再利用できるものをうまく利用する
計算量 (コスト) のスケーリングに気をつける
LAPACK などの汎用ライブラリを積極的に利用する
異言語間結合には少し注意が必要
30 / 31
計算機実験 (講義 L3)
LAPACK ライブラリ
実習・講義予定
実習 EX2
▶
▶
▶
バージョン管理システム
摩擦のあるバネの問題、解の精度
レポート No.1
実習 EX3
▶
▶
▶
ガウスの消去法、ピボット選択
LAPACK の利用
速度の測定
講義 L4
▶
▶
▶
固有値問題
直接法と反復法
Krylov 部分空間法
31 / 31
© Copyright 2024 ExpyDoc
