n×n 行列を考える.
行列積 C = A B :
(行列 C の i 行 j 列目の要素を c_ij
と書けば)
n c_ij = Σ a_ik b_kj k=1
行列・ベクトル積 b = A x :
(ベクトル bの i 行目の要素を b_i
と書けば)
n b_i = Σ a_ik x_k k=1
連立一次方程式の解: 与えられた行列 A と ベクトル b に対して、 A x = b を満たす ベクトル x を求める。
LU分解 とは
行列 A を、下三角行列 L、上三角行列 U を 用いて A = LU
と分解すること。
ただし、 L については l_ij
= 0 (1 ≦ i < j ≦ n),
Uについては u_ij
= 0 (1 ≦ j < i ≦ n)である。
LU 分解を用いた連立一次方程式の解法: A = LU と分解した後、 L y = b を yについて解き、さらに、 U x = y を xについて解く。 L, U とも三角行列なので簡単に解ける。 ( x = U^{-1} L^{-1} b)
l_ii
=1), クラウト(Crout)型(u_ii
=1) クラウト型で考えると A=LU の A, L, Uの各成分について
j a_ij = Σ l_ik u_kj (i ≧ j) k=1 i a_ij = Σ l_ik u_kj (i < j) k=1であるから、総和の最後の項を分けると、
j-1 a_ij = Σ l_ik u_kj + l_ij (i ≧ j) k=1 i-1 a_ij = Σ l_ik u_kj + l_ii u_ij (i < j) k=1となり、これから、関係式(☆):
j-1 l_ij = a_ij - Σ l_ik u_kj (i ≧ j) k=1 i-1 u_ij = (a_ij - Σ l_ik u_kj) / l_ii (i < j) k=1が成り立つことがわかる。
a(i,j)
の次の位置にa(i+1,j)
a[i][j]
の次の位置にa[i][j+1]
l_ii
= 1)が有利だが、
本演習では C 言語を用いるので、クラウト型(u_ii
= 1)を用いている。
密行列に対するアルゴリズムのオーダはO(n^3)である。 通常,行列AをLU分解した結果のLとUは,どちらかの対角成分は 1なので,行列AをLとUで上書きする形でLU分解の結果を得ることが多い. 関係式☆に基づくアルゴリズムなら どんなアルゴリズムでも(例えば,ブロック化や再帰アルゴリズムでも) 正しい解は得られるが、単純な3重ループの場合, 最外ループについて Dongarra氏 流に分類すると次のようになる:
u_ii
= 1)なら) down-looking と呼びたいが,それは普及した呼び方とはいえない.
ガウス法と呼ばれる。
数学とは違う意味の外積を用いる。u_ii
= 1)なら、
m=1, 2, ... n について以下の関係式に基づく。
l_ij = a_ij^(m) (j=m, i=m,...,n) u_ij = a_ij^(m)/l_ii (i=m, j=m+1,...,n) a_ij^(m+1) = a_ij^(m) - l_ik u_kj (i=m+1,...,n, j=m+1,...,n, k=m)アルゴリズムは、m=1, 2, ... n-1 と進める形となり、これに k を用いると、 (添え字の使い方を 0 からに変えて)、 このようになる.
u_ii
= 1)なら) up-lookingと呼びたいが,それは普及した呼び方とはいえない.
ガウス法に基づく。内積や中間積を用いる。u_ii
= 1)なら、
m=1, 2, ... n について以下の関係式に基づく。
for 1≦s<m l_ij = a_ij^(s) (i=m, j=s) a_ij^(s+1) = a_ij^(s) - l_ik u_kj (i=m, j=s+1,...,n, k=s) l_ii = a_ii^(m) (i=m) u_ij = a_ij^(m)/l_ii (i=m, j=m+1,...,n)アルゴリズムは、m=1, 2, ... n (さらに 各 m について s = 1, 2, ... m-1) と進める形となり、これに i (さらに k)を用いると、 (添え字の使い方を 0 からに変えて)、 このようになる.
それぞれのアルゴリズムを理解するには, 2次元配列(行列)を図に書いてみて, 各計算ステップで読み出す要素と書き込む要素の対応を考えるとよい. (常人が理解するには、図を書くべきである。)
例: right-looking (down-looking?) アルゴリズムの場合 n=8のとき: k=0, 1 の場合の処理後,k = 2のとき: for(j=k+1;j<n;j++) a[k][j] /= a[k][k]; は a[][] において ........ ........ ..*!!!!! ........ ........ ........ ........ ........ (*から読み出し.各!からの読み出しとそこへの書き込み) for(i=k+1;i<n;i++) for(j=k+1;j<n;j++) a[i][j] -= a[i][k] * a[k][j]; は a[][] において ........ ........ ...***** ..*!!!!! ..*!!!!! ..*!!!!! ..*!!!!! ..*!!!!! (*から読み出し.各!からの読み出しとそこへの書き込み)
left-looking法などの用語は、
Fortran言語で Doolittle型(l_ii
= 1) を用いる場合に使われてきたものであることに注意。