le4 parallel programming

LU分解 (LU decomposition)

数学的定義

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)

LU分解の分類

Doolittle型(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

が成り立つことがわかる。

配列の添字と連続アクセス

Fortran言語とC言語で異なる。

Fortran言語では、a[i,j]の次の位置にa[i+1,j]
C言語では a[i][j] の次の位置にa[i][j+1]
Doolittle型, クラウト型のどちらが有利 (メモリアクセスが連続的になるか)が異なる

LU分解の3種類のアルゴリズム

密行列に対するアルゴリズムのオーダはO(n^3)である。通常，行列AをLU分解した結果のLとUは，どちらかの対角成分は 1なので，行列AをLとUで上書きする形でLU分解の結果を得ることが多い．関係式☆に基づくアルゴリズムならどんなアルゴリズムでも正しい解は得られるが、 3重ループの最外ループについて Dongarra氏流に分類すると次のようになる:

right-looking法: 配列の添え字によってはdown-looking。ガウス法と呼ばれる。数学とは違う意味の外積を用いる。
クラウト型(u_ii = 1)なら:
```
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)
```
アルゴリズムはこのようになる．
left-looking法: 配列の添え字によっては up-looking。ガウス法に基づく。内積や中間積を用いる。
クラウト型(u_ii = 1)なら:
```
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)
u_ij = a_ij^(m)/l_ii                   (i=m, j=m+1,...,n)
```
アルゴリズムはこのようになる．
Crout法: upもleftも見る。Doolittle型でも同様となる。
rigt-looking法では、 L と U について求まった最新の要素を使って A のまだ分解が済んでいない部分を更新してしまう。
更新をすぐにはせずに、関係式☆を直接使って L と U の次の要素を求めるようにしたものである。
アルゴリズムはこのようになる．

それぞれのアルゴリズムを理解するには， 2次元配列(行列)を図に書いてみて，各計算ステップで読み出す要素と書き込む要素の対応を考えるとよい．

実用上の考慮

実用上は:

部分ピボッティング(partial pivoting)
スケーリング(scaling)
残差補正

を考慮する必要があるが、本実験演習では無視することにする。無視するためには入力となる行列Aがある意味で「よい」性質を持つ必要がある。

並列プログラミング, 先頭ページへ

Masahiro Yasugi: yasugi@kuis.kyoto-u.ac.jp