数値計算第 3 回課題 (大石先生)

{G99P043-4 \atop 河邊昌彦}

出題日 : 2001/06/01
提出期限 : 2001/06/15
提出日 : 2001/06/15

1 課題

1. 教科書「数値計算」 p.31 の問題 3 を解け。以下にその問題を示す。

(1) ガウスの消去法のプログラムを作り、連立一次方程式の近似解を求めてみよ。最初はピボット選択なしで、次に部分ピボット選択をするプログラムを作れ。

(2) LU 分解法のプログラムを作り、連立一次方程式の近似解を求めてみよ。

(3) LU 分解法を利用して逆行列を求めるプログラムを作れ。

(4) 行列

A = [\begin{matrix}1 & 2 & 3 \\ 2 & 8 & 11 \\ 3 & 22 & 35 \end{matrix}]

の LU 分解を求めよ。

(5) 4 の結果を用いて

(1) b = [\begin{matrix}1 \\ 0 \\ 0 \end{matrix}]  (2) b = [\begin{matrix}1 \\ 2 \\ 3 \end{matrix}]  (3) b = [\begin{matrix}−1 \\ 0 \\ 1 \end{matrix}]

として

Ax = b

を解け。

2. 以下の性質 3, 5, 6 を証明せよ。性質 1, 2, 4 は仮定してよい。

性質 1 n × n 行列 A が ∥A∥ < 1 を満たすとする。このとき、 I を単位行列とすると、

I + A は正則行列となる。

性質 2 A を n × n 行列とする。 ∥A∥ が誘導されたノルムとすると、元のベクトルのノルム ∥ · ∥ に対して、

∥Ax∥ ≤ ∥A∥∥x∥

性質 3 A を n × n 行列、 ∥A∥ を誘導されたノルムとする。 R を n × n 行列、 I を n × n 単位行列とする。

∥RA − I∥ < 1 ならば A は正則であり、

∥A^{−1}∥ ≤ \frac{∥R∥}{1 − ∥RA − I∥}

が成り立つ。性質 4 A を n × n 行列、 ∥A∥ を誘導されたノルムとする。このとき、 λ を A の任意の固有値とすると、

|λ| ≤ ∥A∥

が成り立つ。

性質 5 (ゲルシュゴリンの定理)A を n×n(複素) 行列とする。このとき、行列 A の固有値は

K = \{λ ∈ C | |λ − a_{ii}| ≤ ∑_{\begin{matrix}k=1 \\ k≠i \end{matrix}}^{n} |a_{ik}|\}

の中にある。

性質 6 IEEE754 の下での浮動小数点演算の下で議論する。最近点での丸めのモードで考える。

u = 2^{−53}

F で normalized 及び denormalized 及び 0 の double の集合を考える。

x, y ∈ F
∘ ∈ {+, −, ∗, /}
|δ| ≤ u
x ∘ y = (x ∘ y)(1 + δ)

となる。これを利用すると、

x, y ∈ F^{n}

とすると、 s = x · y (実数としての計算) に対して、

near(); t = 0; for ( i = 0 ; i ≤ n ; i++ )

t += x_{i} ∗ y_{i};

として計算した値 t との間に、

|s − t| ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|
γ_{n} = \frac{nu}{1 − nu}

という関係がある。

2 ガウスの消去法 (1-1)

行列やベクトルに関しては、柏木先生の方の課題でライブラリ (A.1) を作ってあるので、それを使うことにする。

ガウスの消去法により連立一次方程式を解く関数はピボット選択無しが A.2.1、ピボット選択有りが A.2.2 となる。それぞれ、

     template< size_t Dim, class T >
     void gaussian_elimination( Vector< Dim, T > &result,
                                Matrix< Dim, Dim, T > a, Vector< Dim, T > b );
     template< size_t Dim, class T >
     bool gaussian_elimination_pivot( Vector< Dim, T > &result,
                                      Matrix< Dim, Dim, T > a, Vector< Dim, T > b );

において、 a を係数行列、 b を定数ベクトルとする連立一次方程式の解を、 result に返す。

また、ピボット選択をする場合、ピボットの係数の絶対値が 2^{−50} よりも小さい場合は、失敗として false を返す。実際には、この閾値は、係数行列の各要素の精度等に応じて決まるので、定数ではない。しかし、ここでは簡単のため、 2^{0} 程度の大きさの数を扱うということを前提にして、定数とした。

3 LU 分解法 (1-2)

3.1 LU 分解

LU 分解 (A.3.1) は、ガウスの消去法の a に関する部分だけを抜き出せばよい。ただしそのままでは、下三角行列と上三角行列が両方とも a に重ね書きされているので、それを分ける必要がある。また、ピボット選択も行うので、行の置換の情報も返さなければならない。

LU 分解を行う関数は次のようになっている。

     template< size_t Dim, class T >
     bool lu_decomposition( Matrix< Dim, Dim, T > &lower, Matrix< Dim, Dim, T > &upper,
                            Vector< Dim, size_t > &pivot, Matrix< Dim, Dim, T > a );

行列 a の LU 分解を求め、 lower に下三角行列を、 upper に上三角行列を返す。また、行の置換情報を pivot に返すが、その意味は、第 i 列消去の時に i 行目と pivot[i] を入れ換えた、ということである。従って、この行置換の操作を表す行列を P とすると、

P A = LU

である。そして、ガウスの消去法と同様に、ピボット選択に失敗すると false を返す。

3.2 連立一次方程式の解法

係数行列の LU 分解が求まれば、それを使って連立一次方程式を解くことができる。連立一次方程式を

Ax = b

とすると、 LU 分解により、

LUx = P b

である。従って、 P b を求め、 Ly = P b の解 y を求め、最後に Ux = y の解 x を求めればよい。このようにして連立一次方程式を解く関数が A.3.2 である。

     template< size_t Dim, class T >
     void lu_solve( Vector< Dim, T > &x, const Matrix< Dim, Dim, T > &lower,
                    const Matrix< Dim, Dim, T > &upper, const Vector< Dim, size_t > &pivot,
                    Vector< Dim, T > b );

LU 分解で得られた lower, upper, pivot と、定数ベクトルの b を与えると、その連立一次方程式の解を x に返す。

4 逆行列 (1-3)

LU 分解法は逆行列を求めるのにも使うことができる。

行列 A の逆行列を X = [x_{1} x_{2} ⋯ x_{n}] とし、 I = [e_{1} e_{2} ⋯ e_{n}] とすると、

[e_{1} e_{2} ⋯ e_{n}] = I = AX = A [x_{1} x_{2} ⋯ x_{n}] = [Ax_{1} Ax_{2} ⋯ Ax_{n}]

であるので、 1 ≤ i ≤ n に対して

Ax_{i} = e_{i}

である。従って、各 x_{i} を LU 分解法で求めれば、 A の逆行列を得ることができる。これを行う関数が A.3.3 であり、

     template< size_t Dim, class T >
     bool inverse( Matrix< Dim, Dim, T > &x, const Matrix< Dim, Dim, T > &a );

となっている。

5 実行例 (1-4, 1-5)

5.1 結果

ガウスの消去法 (ピボット選択無し/有り)、 LU 分解法、 LU 分解法による逆行列の関数を実際に使ってみているのが A.5 のプログラムである。これを実行すると、以下のような出力を得る。 (cygwin の gcc 2.95 でオプション “-lstdc++” を付けてコンパイルした)

     a =
           0.00000     -15.5232      14.7239      6.85938      3.85791
         -0.359863      5.89575     -6.69263     -12.0139     -15.7734
           10.7073      12.0906      9.57666     -3.34082     0.783936
           14.2158     -10.0654      12.6670      8.43091      14.1641
           10.8904     -10.3730      14.1353     -8.00977     -1.12402
     b =
          -12.8843     -12.8826     -4.56714     -10.1074      5.65894

gaussian elimination (w/o pivot selection)
a\b =
   -NaN         -NaN         -NaN         -NaN         -NaN

gaussian elimination (with pivot selection)
a\b =
-2.33233     0.103442     0.160586     -4.30643      4.12048

LU decomposition
l =
1.00000      0.00000      0.00000      0.00000      0.00000
0.753194      1.00000      0.00000      0.00000      0.00000
0.00000    -0.789109      1.00000      0.00000      0.00000
0.766075    -0.135330     0.300717      1.00000      0.00000
-0.0253143     0.286753    -0.432631     0.602355      1.00000
u =
14.2158     -10.0654      12.6670      8.43091      14.1641
0.00000      19.6718    0.0359533     -9.69093     -9.88436
0.00000      0.00000      14.7522    -0.787826     -3.94192
0.00000      0.00000      0.00000     -15.5430     -12.1270
0.00000      0.00000      0.00000      0.00000     -6.98115
pivot =
      3            2            3            4            4
p =
0.00000      0.00000      0.00000      1.00000      0.00000
0.00000      0.00000      1.00000      0.00000      0.00000
1.00000      0.00000      0.00000      0.00000      0.00000
0.00000      0.00000      0.00000      0.00000      1.00000
0.00000      1.00000      0.00000      0.00000      0.00000
p^-1 * l * u =
0.00000     -15.5232      14.7239      6.85938      3.85791
-0.359863      5.89575     -6.69263     -12.0139     -15.7734
10.7073      12.0906      9.57666     -3.34082     0.783936
14.2158     -10.0654      12.6670      8.43091      14.1641
10.8904     -10.3730      14.1353     -8.00977     -1.12402

x = a\b =
-2.33233     0.103442     0.160586     -4.30643      4.12048

a^-1 =
-0.00887746    0.0932906  -0.00282690     0.103496   -0.0374130
-0.000941757   -0.0168583    0.0520111   -0.0231048   -0.0215337
0.0489903   -0.0323072    0.0500915   -0.0508225    0.0160246
0.0879430     0.111761    0.0195316    0.0889775    -0.131657
-0.0879181    -0.143243   -0.0166250   -0.0572035    0.0862831

a * a^-1 =
     1.00000 -3.90313e-18  6.47269e-17 -4.37340e-17 -2.87314e-17
-3.64942e-16      1.00000 -3.47216e-17 -7.52436e-17  2.20527e-16
-1.14959e-16  2.54178e-17      1.00000  8.37207e-17 -6.84199e-17
-5.85469e-18 -1.40296e-16 -1.62901e-17      1.00000  1.42464e-16
-3.11789e-16 -1.15657e-16 -7.34716e-17 -8.56249e-17      1.00000

a * a^-1 =
[     0.99999999999998856      1.0000000000000113 ]
 [ -1.3433698597964394e-14  1.3322676295501878e-14 ]
  [ -4.5241588253475129e-15  4.6768144912334719e-15 ]
   [ -1.0463852007092100e-14  1.0547118733938987e-14 ]
    [ -1.1601830607332886e-14  1.1768364061026659e-14 ]
[ -1.2656542480726785e-14  1.2434497875801753e-14 ]
 [     0.99999999999998535      1.0000000000000147 ]
  [ -4.8849813083506888e-15  4.7739590058881731e-15 ]
   [ -1.1768364061026659e-14  1.1546319456101628e-14 ]
    [ -1.2656542480726785e-14  1.3100631690576847e-14 ]
[ -1.6764367671839864e-14  1.6694978732800791e-14 ]
 [ -2.0539125955565396e-14  1.9900747716405931e-14 ]
  [     0.99999999999999289      1.0000000000000073 ]
   [ -1.6035783811929605e-14  1.5827616994812388e-14 ]
    [ -1.7624790515924360e-14  1.7291723608536813e-14 ]
[ -2.5757174171303632e-14  2.5757174171303632e-14 ]
 [ -3.1086244689504383e-14  3.1086244689504383e-14 ]
  [ -1.0269562977782698e-14  1.0241807402167069e-14 ]
   [     0.99999999999997569      1.0000000000000244 ]
    [ -2.6645352591003757e-14  2.6645352591003757e-14 ]
[ -1.9026447084513620e-14  1.8651746813702630e-14 ]
 [ -2.2565282975506307e-14  2.2787327580431338e-14 ]
  [ -7.6917638924811627e-15  7.5842110369706006e-15 ]
   [ -1.7902346272080649e-14  1.7708057242771247e-14 ]
    [     0.99999999999998035      1.0000000000000195 ]

solutions
     1.00000      0.00000      0.00000
    0.666667      1.00000      0.00000
    0.333333     0.800000      1.00000
     3.00000      22.0000      35.0000
     0.00000     -6.66667     -12.3333
     0.00000      0.00000      1.20000
     1.58333     -1.54167     0.833333
     1.00000 -6.93889e-17  5.55112e-17
    -1.66667      1.33333    -0.666667

a の (1, 1) 成分が 0 なので、ピボット選択無しではガウスの消去法はできない。また、 P^{−1}LU = A となっているようなので、正しく LU 分解されているといえるであろう。

逆行列を求め、元の行列と掛けてみると、多少の誤差が出た。そこで、区間演算クラス (A.4) を使って、その誤差を評価してみた。その結果によると、単位行列 I はその区間の中に入っていることが分かる。従って、アルゴリズムに誤りが無いと見なしてよいであろう。

出力の最後の部分は、問題 1-4, 1-5 の解になっている。以下でそれを確かめる。

5.2 検証 (1-4)

問題 1-4 を手で解いてみる。 U はその作られ方から、 A を何度か基本変形 (行を定数倍、行の交換、行の定数倍を他の行に加算) することによって得られる。そこで、それらの基本変形を表す基本行列を S_{1}, ⋯, S_{m} とすると、

S_{m} ⋯ S_{1}A = U
P A = LU
∴ L^{−1}P = S_{m} ⋯ S_{1}
(L^{−1}P)^{−1} = P^{−1}L = (S_{m} ⋯ S_{1})^{−1}

L は下三角行列で、 P は行交換行列なので、 P^{−1}L から L と P を求めることができる。従って、基本変形を “∼” で書くと、

[A I] = [\begin{matrix}1 & 2 & 3 & 1 & 0 & 0 \\ 2 & 8 & 11 & 0 & 1 & 0 \\ 3 & 22 & 35 & 0 & 0 & 1 \end{matrix}] ∼ [\begin{matrix}3 & 22 & 35 & 0 & 0 & 1 \\ 2 & 8 & 11 & 0 & 1 & 0 \\ 1 & 2 & 3 & 1 & 0 & 0 \end{matrix}]

∼ [\begin{matrix}3 & 22 & 35 & 0 & 0 & 1 \\ 0 & −20/3 & −37/3 & 0 & 1 & −2/3 \\ 0 & −16/3 & −26/3 & 1 & 0 & −1/3 \end{matrix}]

∼ [\begin{matrix}3 & 22 & 35 & 0 & 0 & 1 \\ 0 & −20/3 & −37/3 & 0 & 1 & −2/3 \\ 0 & 0 & 6/5 & 1 & −4/5 & 1/5 \end{matrix}]

となって、

U = [\begin{matrix}3 & 22 & 35 \\ 0 & −20/3 & −37/3 \\ 0 & 0 & 6/5 \end{matrix}]  L^{−1}P = [\begin{matrix}0 & 0 & 1 \\ 0 & 1 & −2/3 \\ 1 & −4/5 & 1/5 \end{matrix}]

である。次に、

[L^{−1}P I] = [\begin{matrix}0 & 0 & 1 & 1 & 0 & 0 \\ 0 & 1 & −2/3 & 0 & 1 & 0 \\ 1 & −4/5 & 1/5 & 0 & 0 & 1 \end{matrix}] ∼ [\begin{matrix}1 & −4/5 & 1/5 & 0 & 0 & 1 \\ 0 & 1 & −2/3 & 0 & 1 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \end{matrix}]

∼ [\begin{matrix}1 & −4/5 & 0 & −1/5 & 0 & 1 \\ 0 & 1 & 0 & 2/3 & 1 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \end{matrix}]

∼ [\begin{matrix}1 & 0 & 0 & 1/3 & 4/5 & 1 \\ 0 & 1 & 0 & 2/3 & 1 & 0 \\ 0 & 0 & 1 & 1 & 0 & 0 \end{matrix}]

なので、

P^{−1}L = [\begin{matrix}1/3 & 4/5 & 1 \\ 2/3 & 1 & 0 \\ 1 & 0 & 0 \end{matrix}]

よって、

L = P [\begin{matrix}1/3 & 4/5 & 1 \\ 2/3 & 1 & 0 \\ 1 & 0 & 0 \end{matrix}] = [\begin{matrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{matrix}] [\begin{matrix}1/3 & 4/5 & 1 \\ 2/3 & 1 & 0 \\ 1 & 0 & 0 \end{matrix}] = [\begin{matrix}1 & 0 & 0 \\ 2/3 & 1 & 0 \\ 1/3 & 4/5 & 1 \end{matrix}]

とすればよい。つまり、

L = [\begin{matrix}1 & 0 & 0 \\ 2/3 & 1 & 0 \\ 1/3 & 4/5 & 1 \end{matrix}] U = [\begin{matrix}3 & 22 & 35 \\ 0 & −20/3 & −37/3 \\ 0 & 0 & 6/5 \end{matrix}] P = [\begin{matrix}0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{matrix}]

である。

5.3 検証 (1-5)

(1)

b = [\begin{matrix}1 \\ 0 \\ 0 \end{matrix}]

なので、

P b = [\begin{matrix}0 \\ 0 \\ 1 \end{matrix}]

である。 Ly = P b の解を求めると、

[L P b] = [\begin{matrix}1 & 0 & 0 & 0 \\ 2/3 & 1 & 0 & 0 \\ 1/3 & 4/5 & 1 & 1 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{matrix}] = [I y]

次に、 Ux = y の解を求めると、

[U y] = [\begin{matrix}3 & 22 & 35 & 0 \\ 0 & −20/3 & −37/3 & 0 \\ 0 & 0 & 6/5 & 1 \end{matrix}] ∼ [\begin{matrix}3 & 22 & 0 & −175/6 \\ 0 & −20/3 & 0 & 185/18 \\ 0 & 0 & 1 & 5/6 \end{matrix}]

∼ [\begin{matrix}3 & 22 & 0 & −175/6 \\ 0 & 1 & 0 & −37/24 \\ 0 & 0 & 1 & 5/6 \end{matrix}]

∼ [\begin{matrix}3 & 0 & 0 & 19/4 \\ 0 & 1 & 0 & −37/24 \\ 0 & 0 & 1 & 5/6 \end{matrix}]

∼ [\begin{matrix}1 & 0 & 0 & 19/12 \\ 0 & 1 & 0 & −37/24 \\ 0 & 0 & 1 & 5/6 \end{matrix}] = [I x]

よって、

x = [\begin{matrix}19/12 \\ −37/24 \\ 5/6 \end{matrix}]

である。

(2) 同様に解くと、

[L P b] = [\begin{matrix}1 & 0 & 0 & 3 \\ 2/3 & 1 & 0 & 2 \\ 1/3 & 4/5 & 1 & 1 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 3 \\ 0 & 1 & 0 & 0 \\ 0 & 4/5 & 1 & 0 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 3 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}] ∼ [I y]

[U y] = [\begin{matrix}3 & 22 & 35 & 3 \\ 0 & −20/3 & −37/3 & 0 \\ 0 & 0 & 6/5 & 0 \end{matrix}] ∼ [\begin{matrix}3 & 0 & 0 & 3 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}] ∼ [I x]

∴ x = [\begin{matrix}1 \\ 0 \\ 0 \end{matrix}]

(3) 同様に、

[L P b] = [\begin{matrix}1 & 0 & 0 & 1 \\ 2/3 & 1 & 0 & 0 \\ 1/3 & 4/5 & 1 & −1 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 1 \\ 0 & 1 & 0 & −2/3 \\ 0 & 4/5 & 1 & −4/3 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 1 \\ 0 & 1 & 0 & −2/3 \\ 0 & 0 & 1 & −4/5 \end{matrix}] ∼ [I y]

[U y] = [\begin{matrix}3 & 22 & 35 & 1 \\ 0 & −20/3 & −37/3 & −2/3 \\ 0 & 0 & 6/5 & −4/5 \end{matrix}] ∼ [\begin{matrix}3 & 22 & 0 & −73/3 \\ 0 & −20/3 & 0 & 80/9 \\ 0 & 0 & 1 & −2/3 \end{matrix}] ∼ [\begin{matrix}1 & 0 & 0 & 5/3 \\ 0 & 1 & 0 & −4/3 \\ 0 & 0 & 1 & −2/3 \end{matrix}] ∼ [I x]

∴ x = [\begin{matrix}5/3 \\ −4/3 \\ −2/3 \end{matrix}]

5.4 まとめ

数値計算により得られた解と、手で計算して得れた解がほぼ一致することが分かる。従って、厳密な証明はしていないが、数値計算によってほぼ正確な解が得られると言ってもよいであろう。

6 性質 3 (2-3)

補題 6.1 一般の正方行列 A が正則でないならば、 A は零因子である。証明 (補題 6.1) A を n 次非正則正方行列とする。即ち、

rank A < n である。また、線形写像の次元定理より

dim Ker A + rank A = n であるので、

dim Ker A = n − rank A > 0 となる。従って、 Ker A には零でないベクトルが含まれている。そこで、

B = [b_{1} ⋯ b_{n}]
o ≠ b_{i} ∈ Ker A

とすると、 B ≠ O で、

AB = A[b_{1} ⋯ b_{n}]
= [Ab_{1} ⋯ Ab_{n}]
= [o ⋯ o]
= O

であるので、 A は零因子である。 (証明終) 補題 6.2 A, B を n 次正方行列とする。 AB が正則ならば B は正則である。証明 (補題 6.2) 背理法により証明する。

AB が正則であって、かつ B は正則でないとする。すると、補題 6.1 より、 B は零因子であるので、 BC = O となる零行列でない C が存在する。また、 AB は正則なので、その逆行列 (AB)^{−1} が存在する。従って、

I = (AB)^{−1} · AB 両辺に右から C を掛けると、

C = (AB)^{−1}AB · C
= (AB)^{−1}A · O
= O

となって矛盾である。従って、 B は正則である。 (証明終)

補題 6.3 ∥ · ∥ を誘導されたノルムとし、 A, B を n 次正方行列とする。すると、

∥AB∥ ≤ ∥A∥∥B∥

である。証明 (補題 6.3) 性質 2 より、 n 次元ベクトル x に対して、

∥ABx∥ = ∥A(Bx)∥ ≤ ∥A∥∥Bx∥ ≤ ∥A∥∥B∥∥x∥ である。そして、誘導されたノルムの定義より、

∥AB∥ = sup_{∥x∥=1} ∥ABx∥ ≤ sup_{∥x∥=1} ∥A∥∥B∥∥x∥ = ∥A∥∥B∥

となる。 (証明終) 系 6.1 ∥ · ∥ を誘導されたノルムとする。 A, B を n 次正方行列とし、 B は正則であるとする。この時、

∥AB∥ ≥ \frac{∥A∥}{∥B^{−1}∥}

である。証明 (系 6.1) 補題 6.3 で、 A = AB, B = B^{−1} と置き直せばよい。

∥AB∥∥B^{−1}∥ ≥ ∥AB · B^{−1}∥
∴ ∥AB∥ ≥ \frac{∥A∥}{∥B^{−1}∥}

(証明終) 証明 (性質 3) まず、 A が正則であることを示す。 ∥RA − I∥ < 1 と性質 1 より、

RA − I + I = RA は正則である。従って、補題 6.2 より A は正則である。

ここで、

R = R′ + A^{−1} とおくと、

RA − I = (R′ + A^{−1})A − I = R′A であり、条件から

∥R′A∥ < 1 となる。系 6.1 を使うと、

1 > ∥R′A∥ ≥ \frac{∥R′∥}{∥A^{−1}∥}

従って、

0 < 1 − ∥R′A∥ ≤ 1 − \frac{∥R′∥}{∥A^{−1}∥}
\frac{1}{1 − ∥R′A∥} ≥ \frac{1}{1 − ∥R′∥/∥A^{−1}∥} = \frac{∥A^{−1}∥}{∥A^{−1}∥ − ∥R′∥}

ノルムの公準である三角不等式を使うと、

\frac{∥R∥}{1 − ∥RA − I∥} = \frac{∥R′ + A^{−1}∥}{1 − ∥R′A∥} ≥ ∥R′ + A^{−1}∥ · \frac{∥A^{−1}∥}{∥A^{−1}∥ − ∥R′∥}

≥ (∥A^{−1}∥ − ∥R′∥) · \frac{∥A^{−1}∥}{∥A^{−1}∥ − ∥R′∥}

= ∥A^{−1}∥

よって、

∥A^{−1}∥ ≤ \frac{∥R∥}{1 − ∥RA − I∥}

(証明終)

7 性質 5 (2-5)

定義 7.1 n 次複素正方行列 A = (a_{ij}) が n > 1, a_{11} ≠ 0 とする。 A に対して、第一列目に関する列の消去を行った行列を A′ とする。即ち、

r_{j} = \frac{a_{1j}}{a_{11}}

とすると、

A′ = [\begin{matrix}a_{11} & 0 & ⋯ & 0 & ⋯ & 0 \\ a_{21} & a_{22} − r_{2}a_{21} & ⋯ & a_{2j} − r_{j}a_{21} & ⋯ & a_{2n} − r_{n}a_{21} \\ ⋮ & ⋮ & ... & ⋮ & & ⋮ \\ a_{i1} & a_{i2} − r_{2}a_{i1} & ⋯ & a_{ij} − r_{j}a_{i1} & ⋯ & a_{in} − r_{n}a_{i1} \\ ⋮ & ⋮ & & ⋮ & ... & ⋮ \\ a_{n1} & a_{n2} − r_{2}a_{n1} & ⋯ & a_{nj} − r_{j}a_{n1} & ⋯ & a_{nn} − r_{n}a_{n1} \end{matrix}]

である。ここで、 A′ を次のようにブロック化する。

A′ = [\begin{matrix}a_{11} & O \\ C & B \end{matrix}]

O は 1 × (n − 1) の零行列であり、 C は (n − 1) × 1 行列、 B は (n − 1) 次正方行列である。一般の A から B を得る操作を ρ とし、

B = ρ(A)

と書くことにする。補題 7.1 定義 7.1 において、

det A = |a_{11}| det ρ(A) である。証明 (補題 7.1) 行列式の双線形性から

det A = det A′ である。また、 O が零行列であるから

det A′ = det [a_{11}] · det ρ(A)
= |a_{11}| det ρ(A)

よって、

det A = |a_{11}| det ρ(A) (証明終) 定義 7.2 (行対角優位条件)n 次複素正方行列 A = (a_{ij}) に対して、

P (A) ⟺ 1 ≤ ∀i ≤ n : ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} |a_{ij}| < |a_{ii}|

と述語 P を定義する。この P を行対角優位条件という。系 7.1 P (A) ならば、 |a_{ii}| > 0 である。補題 7.2 n > 1 とすると、 n 次複素正方行列 A において、

P (A) ⇒ P (ρ(A)) である。証明 (補題 7.2) P (A) を仮定する。定義 7.1 より、

ρ(A) = (a_{ij} − r_{j}a_{i1})  (2 ≤ i, j ≤ n)

である。従って、

P (ρ(A)) ⟺ 2 ≤ ∀i ≤ n : ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij} − r_{j}a_{i1}| < |a_{ii} − r_{i}a_{i1}|

となる。 2 ≤ ∀i ≤ n に対して、

\begin{aligned}|a_{ii} − r_{i}a_{i1}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij} − r_{j}a_{i1}| & ≥ |a_{ii}| − |r_{i}a_{i1}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} (|a_{ij}| + |r_{j}a_{i1}|) \\ & = |a_{ii}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij}| − ∑_{j=2}^{n} |r_{j}a_{i1}| \end{aligned}  (7.1)

である。ここで、

∑_{j=2}^{n} |r_{j}a_{i1}| = ∑_{j=2}^{n} |\frac{a_{1j}}{a_{11}} a_{i1}|

= \frac{|a_{i1}|}{|a_{11}|} · ∑_{j=2}^{n} |a_{1j}|

< \frac{|a_{i1}|}{|a_{11}|} · |a_{11}|

= |a_{i1}|

なので、 (7.1) は、

\begin{aligned}|a_{ii}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij}| − ∑_{j=2}^{n} |r_{j}a_{i1}| & > |a_{ii}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij}| − |a_{i1}| \\ & = |a_{ii}| − ∑_{\begin{matrix}j=1 \\ j≠i \end{matrix}}^{n} |a_{ij}| \\ & > 0 \end{aligned}  (7.2)

となる。 (7.1) と (7.2) から、

|a_{ii} − r_{i}a_{i1}| − ∑_{\begin{matrix}j=2 \\ j≠i \end{matrix}}^{n} |a_{ij} − r_{j}a_{i1}| > 0

であるので、 P (ρ(A)) が示された。 (証明終) 補題 7.3 n 次複素正方行列 A に対して、

P (A) ⇒ det A ≠ 0 である。証明 (補題 7.3) P (A) を仮定する。 A に対して補題 7.2 を繰り返し適用すると、

1 ≤ ∀i ≤ n − 1 : P (A) ⇒ P (ρ^{i}(A)) となる。従って、 ρ^{i}(A) の (1, 1) 成分を a_{i} とすると、系 7.1 より

0 ≤ ∀i ≤ n − 1 : |a_{i}| > 0 (7.3) また、補題 7.1 から

det A = |a_{0}| det ρ(A) = |a_{0}||a_{1}| det ρ^{2}(A) = ⋯ = (∏_{i=0}^{n−2} |a_{i}|) det ρ^{n−1}(A) = ∏_{i=0}^{n−1} |a_{i}| (7.4)

よって、 (7.3) と (7.4) より

det A ≠ 0 である。 (証明終)

証明 (性質 5) 背理法により証明する。 λ を A の固有値とすると、ゲルシュゴリンの定理の否定は、

¬∃i 1 ≤ i ≤ n ∧ |λ − a_{ii}| ≤ ∑_{\begin{matrix}k=1 \\ k≠i \end{matrix}}^{n} |a_{ik}| ⟺ ∀i 1 ≤ i ≤ n ⇒ |λ − a_{ii}| > ∑_{\begin{matrix}k=1 \\ k≠i \end{matrix}}^{n} |a_{ik}|

となる。即ち、 A − λI が行対角優位条件を満たすということであるので、

P (A − λI) である。従って、補題 7.3 より

det(A − λI) ≠ 0 であるが、これは、固有方程式

det(A − λI) = 0 と矛盾する。よって、ゲルシュゴリンの定理が成り立つ。 (証明終)

8 性質 6 (2-6)

定義 8.1 s, t を次のように帰納的に表す。

\{ \begin{matrix}s_{0} & = 0 \\ s_{i} & = s_{i−1} + x_{i}y_{i} & (1 ≤ i ≤ n) \\ s & = s_{n} \end{matrix}
\{ \begin{matrix}t_{0} & = 0 \\ t_{i} & = t_{i−1} + (x_{i} ∗ y_{i}) & (1 ≤ i ≤ n) \\ t & = t_{n} \end{matrix}

補題 8.1 1 ≤ i ≤ n について、

|t_{i}| ≤ (1 + u)^{i} ∑_{j=1}^{i} |x_{j}y_{j}|

である。証明 (補題 8.1) 数学的帰納法で証明する。

i = 1 に関しては、

t_{1} = t_{0} + (x_{1} ∗ y_{1}) = 0 + (x_{1} ∗ y_{1}) = x_{1} ∗ y_{1}

であるので、 |δ| ≤ u なる δ が存在して、

t_{1} = x_{1} ∗ y_{1} = x_{1}y_{1}(1 + δ)

である。従って、

|t_{1}| = |x_{1}y_{1}(1 + δ)| ≤ |x_{1}y_{1}|(1 + |δ|) = |x_{1}y_{1}|(1 + u)

となって、 i = 1 で成立する。次に、 i = k − 1 (k ≥ 2) で成立を仮定する。すると、 |δ|, |δ′| ≤ u なる δ, δ′ が存在して、

|t_{k}| = |t_{k−1} + (x_{k} ∗ y_{k})|
= |(t_{k−1} + x_{k}y_{k}(1 + δ))(1 + δ′)|
≤ |1 + δ′||t_{k−1}| + |1 + δ′||1 + δ||x_{k}y_{k}|
≤ (1 + u)|t_{k−1}| + (1 + u)^{2}|x_{k}y_{k}|

ここで、 i = k − 1 での結論を使うと、

|t_{k}| ≤ (1 + u) ((1 + u)^{k−1} ∑_{j=1}^{k−1} |x_{j}y_{j}|) + (1 + u)^{2}|x_{k}y_{k}|

≤ (1 + u)^{k} ∑_{j=1}^{k−1} |x_{j}y_{j}| + (1 + u)^{k}|x_{k}y_{k}|

= (1 + u)^{k} ∑_{j=1}^{k} |x_{j}y_{j}|

よって、 i = k の時も成り立つ。以上により、 1 ≤ i ≤ n に対して、

|t_{i}| ≤ (1 + u)^{i} ∑_{j=1}^{i} |x_{j}y_{j}|

である。 (証明終)

補題 8.2 1 ≤ i ≤ n に関して、

|s_{i} − t_{i}| ≤ ((1 + u)^{i} − 1) ∑_{j=1}^{i} |x_{j}y_{j}|

である。

証明 (補題 8.2) 同様に数学的帰納法で証明する。 i = 1 の時は、 |δ| ≤ u なる δ が存在して、

|s_{1} − t_{1}| = |x_{1}y_{1} − (x_{1} ∗ y_{1})|
= |x_{1}y_{1} − x_{1}y_{1}(1 + δ)|
= |δ||x_{1}y_{1}|
≤ u|x_{1}y_{1}|

よって、 i = 1 で成立する。 i = k − 1 (k ≥ 2) での成立を仮定する。 |δ|, |δ′| ≤ u なる δ, δ′ が存在して、

|s_{k} − t_{k}| = |(s_{k−1} + x_{k}y_{k}) − (t_{k−1} + (x_{k} ∗ y_{k}))|
= |(s_{k−1} + x_{k}y_{k}) − (t_{k−1} + x_{k}y_{k}(1 + δ))(1 + δ′)|
= |s_{k−1} − t_{k−1} − δ′t_{k−1} − δ(1 + δ′)x_{k}y_{k}|
≤ |s_{k−1} − t_{k−1}| + |δ′||t_{k−1}| + |δ||1 + δ′||x_{k}y_{k}|
≤ |s_{k−1} − t_{k−1}| + u|t_{k−1}| + u(1 + u)|x_{k}y_{k}|

である。ここで、 i = k − 1 の時の結果と、補題 8.1 を使うと、

|s_{k} − t_{k}| ≤ (((1 + u)^{k−1} − 1) ∑_{j=1}^{k−1} |x_{j}y_{j}|) + (u(1 + u)^{k−1} ∑_{j=1}^{k−1} |x_{j}y_{j}|) + u(1 + u)|x_{k}y_{k}|

= ((1 + u)^{k−1} − 1 + u(1 + u)^{k−1}) ∑_{j=1}^{k−1} |x_{j}y_{j}| + u(1 + u)|x_{k}y_{k}|

= ((1 + u)^{k} − 1) ∑_{j=1}^{k−1} |x_{j}y_{j}| + ((1 + u)^{2} − 1) |x_{k}y_{k}|

≤ ((1 + u)^{k} − 1) ∑_{j=1}^{k−1} |x_{j}y_{j}| + ((1 + u)^{k} − 1) |x_{k}y_{k}|

= ((1 + u)^{k} − 1) ∑_{j=1}^{k} |x_{j}y_{j}|

従って、 i = k の時も成立する。

以上より、 1 ≤ i ≤ n で成立することが示された。 (証明終) 補題 8.3 0 ≤ n < u^{−1} について、

γ_{n} = \frac{nu}{1 − nu} ≥ (1 + u)^{n} − 1

である。特に、等号が成り立つのは n = 0 の時のみである。証明 (補題 8.3) 次のような写像 f と数列 {α_{i}}_{i≥0} を考える。

f(x) = nu(1 + x)

\{ \begin{matrix}α_{0} \\ α_{i} = f(α_{i−1}) & (i ≥ 1) \end{matrix}

α_{i} は一般の i(≥ 1) について、

α_{i} = ∑_{j=1}^{i}(nu)^{j} + (nu)^{i}α_{0}

と書ける。

何故ならば、 α_{1} = f(α_{0}) = nu + nuα_{0} であり、 α_{i−1} = ∑_{j=1}^{i−1}(nu)^{j} + (nu)^{i−1}α_{0} ならば、

α_{i} = nu + nu (∑_{j=1}^{i−1}(nu)^{j} + (nu)^{i−1}α_{0}) = nu + ∑_{j=2}^{i}(nu)^{j} + (nu)^{i}α_{0} = ∑_{j=1}^{i}(nu)^{j} + (nu)^{i}α_{0}

となるからである。ところで、

f(γ_{n}) = nu(1 + γ_{n}) = nu (1 + \frac{nu}{1 − nu}) = nu · \frac{1}{1 − nu} = γ_{n}

であるから、 γ_{n} は f の不動点である。従って、 α_{0} = γ_{n} とすると、

α_{i} = γ_{n}  (i ≥ 0)

である。よって、

γ_{n} = α_{n} = ∑_{j=1}^{n}(nu)^{j} + (nu)^{n}γ_{n}

≥ ∑_{j=1}^{n}(nu)^{j}

≥ ∑_{j=1}^{n} n(n − 1)(n − 2) ⋯ (n − j + 1)u^{j}

= ∑_{j=1}^{n} \frac{n!}{(n − j)!}u^{j}

≥ ∑_{j=1}^{n} \frac{n!}{j!(n − j)!}u^{j}

= ∑_{j=1}^{n}_{n}C_{j}u^{j}

= ∑_{j=0}^{n}_{n}C_{j}1^{n−j}u^{j} − 1

= (1 + u)^{n} − 1

となる。また、これらの等号が成り立つのは n = 0 の時のみである。 (証明終)

証明 (性質 6) 補題 8.2 より、

|s − t| ≤ ((1 + u)^{n} − 1) ∑_{i=1}^{n} |x_{i}y_{i}|

また、補題 8.3 より、

((1 + u)^{n} − 1) ∑_{i=1}^{n} |x_{i}y_{i}| ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|

従って、

|s − t| ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|

である。 (証明終)

なお、補題 8.3 により、 0 < n < u^{−1} に関しては、

|s − t| ≤ ((1 + u)^{n} − 1) ∑_{i=1}^{n} |x_{i}y_{i}|

の方が良い誤差評価を与えることが分かる。

9 考察

問題 1 の数値計算の結果は、それだけではどれほどの精度があるのか知ることはできない。従って、ガウスの消去法などの誤差はどれくらいであるか、といったことや、或いは、真の解が別の方法 (手計算など) で求まっている場合に、数値計算の結果が妥当なものであるか、といったことを知るためには誤差の評価が必要である。

個々の場合に関しては、区間演算によって誤差を評価することが可能であるが、計算を繰り返していくと区間の幅が広がってしまうため、必ずしも有効とは限らない。なお、区間演算のクラス (A.4) は、今回間に合わせで作ったものである。従って、無駄な計算や丸めモード変更が多く、特に丸めモードの変更は一演算ごとに行っているので、効率が悪い。

A ソースコード

A.1 行列、ベクトル

template< size_t Dim, class T = double > class Vector;
template< size_t Dim, class T = double > class RowVector;
template< size_t Dim, class T = double > class ColumnVector;
template< size_t Row, size_t Col, class T = double > class Matrix;
///////////////////////////////////////////////////////////////////////////////
// utility functions
template< class T > inline std::complex< T >    conjugate( const std::complex< T > &c )
{ return std::conjg( c ); }
inline double   conjugate( double d ){ return d; }
template< class T > T   random( void );
template<> double       random< double >( void )
{ return ( rand() & 1 ? 1 : -1 ) * ( rand() & 0xffff ) / ( double )0x1000; }
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
class Vector
{
        typedef Vector< Dim, T >        Self;
public:
        // types and constants
        typedef T       Element, Scalar;
        enum { dimension = Dim };

        // generators
        static Self     rand( void );
        static Vector< 0, T >   empty( void ){ return Vector< 0, T >(); }
        // constructors
        Vector( void );
        Vector( const Self &v );
        // access grant operators
        Element         &operator[]( size_t n ){ return elements[n]; }
        const Element   &operator[]( size_t n ) const { return elements[n]; }
        // unary operators
        Self            operator+( void ) const;
        Self            operator-( void ) const;
        // assignment operators
        Self            &operator+=( const Self &v1 );
        Self            &operator-=( const Self &v1 );
        Self            &operator*=( Scalar s );
        Self            &operator/=( Scalar s );
private:
        Element elements[ ( Dim > 0 ? Dim : 1 ) ];
};
// binary operators
template< size_t Dim, class T >
Vector< Dim, T >        operator+( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 );
template< size_t Dim, class T >
Vector< Dim, T >        operator-( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 );
template< size_t Dim, class T >
Vector< Dim, T >        operator*( const Vector< Dim, T > &v1, Vector< Dim, T >::Scalar s );
template< size_t Dim, class T >
Vector< Dim, T >        operator*( Vector< Dim, T >::Scalar s, const Vector< Dim, T > &v1 );
template< size_t Dim, class T >
Vector< Dim, T >        operator/( const Vector< Dim, T > &v1, Vector< Dim, T >::Scalar s );
// equality operators
template< size_t Dim, class T >
bool                    operator==( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 );
template< size_t Dim, class T >
bool                    operator!=( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 );
// vector functions
template< size_t Dim, class T >
Vector< Dim, T >::Scalar        inner_product( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 );
template< size_t Dim, class T >
Vector< Dim, T >::Scalar        sqrabs( const Vector< Dim, T > &v1 );
template< size_t Dim, class T >
Vector< Dim, T >::Scalar        abs( const Vector< Dim, T > &v1 );
template< size_t Dim, class T >
Vector< Dim, T >                normalize( const Vector< Dim, T > &v1 );
template< size_t Dim, class T >
Vector< Dim, T >                conjugate( const Vector< Dim, T > &v1 );
// inserters
template< size_t Dim, class T >

Vector< Dim + 1, T >    operator<<( const Vector< Dim, T > &v1, Vector< Dim, T >::Element e );
template< size_t Dim, class T >
std::ostream            &operator<<( std::ostream &ostr, const Vector< Dim, T > &v );
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
class RowVector
{
typedef RowVector< Dim, T >     Self;
typedef Vector< Dim, T >        Vec;
typedef ColumnVector< Dim, T >  ColVec;
typedef Matrix< 1, Dim, T >     Mat;
public:
// types and constants
typedef Vec::Element    Element;
typedef Vec::Scalar     Scalar;
enum { dimension = Dim };
// generators
static Self     rand( void ){ return Vec::rand(); }
static RowVector< 0, T >        empty( void ){ return Vec::empty(); }
// constructors
RowVector( void ){}
RowVector( const Self &v ) : vec( v.vec ){}
// implicit converters
RowVector( const Vec &v ) : vec( v ){}
RowVector( const Mat &m ) : vec( m[0] ){}
operator Vec( void ) const { return vec; }
operator Mat( void ) const { return ( Matrix< 0, Dim, T >() << *this ); }
// explicit converters
friend Self     row_vector( const Vec &v ){ return v; }
friend Self     row_vector( const Mat &m ){ return m; }
//      friend Vec      vector( const Self &v ){ return v; }
friend Mat      matrix( const Self &v ){ return v; }
// access grant operators
Element         &operator[]( size_t n ){ return vec[n]; }
const Element   &operator[]( size_t n ) const { return vec[n]; }
// unary operators
Self            operator+( void ) const { return ( +vec ); }
Self            operator-( void ) const { return ( -vec ); }
// assignment operators
Self            &operator+=( const Self &v1 ){ vec += v1.vec; return *this; }
Self            &operator-=( const Self &v1 ){ vec -= v1.vec; return *this; }
Self            &operator*=( Scalar s ){ vec *= s; return *this; }
Self            &operator/=( Scalar s ){ vec /= s; return *this; }
// binary operators
friend Self     operator+( const Self &v1, const Self &v2 ){ return ( v1.vec + v2.vec ); }
friend Self     operator-( const Self &v1, const Self &v2 ){ return ( v1.vec - v2.vec ); }
friend Self     operator*( const Self &v1, Scalar s ){ return ( v1.vec * s ); }
friend Self     operator*( Scalar s, const Self &v1 ){ return ( s * v1.vec ); }

friend Self     operator/( const Self &v1, Scalar s ){ return ( v1.vec / s ); }
// equality operators
friend bool     operator==( const Self &v1, const Self &v2 ){ return ( v1.vec == v2.vec ); }
friend bool     operator!=( const Self &v1, const Self &v2 ){ return ( v1.vec != v2.vec ); }
// vector functions
friend Scalar   inner_product( const Self &v1, const Self &v2 )
                        { return inner_product( v1.vec, v2.vec ); }
friend Scalar   sqrabs( const Self &v1 ){ return sqrabs( v1.vec ); }
friend Scalar   abs( const Self &v1 ){ return abs( v1.vec ); }
friend Self     normalize( const Self &v1 ){ return normalize( v1.vec ); }
friend ColVec   transpose( const Self &v ){ return Vec( v ); }
friend ColVec   conjugate( const Self &v ){ return conjugate( Vec( v ) ); }
// inserters
friend RowVector< Dim + 1, T >  operator<<( const Self &v1, Scalar s )
                                        { return ( v1.vec << s ); }
friend std::ostream     &operator<<( std::ostream &ostr, const Self &v )
                                { return ( ostr << matrix( v ) ); }
private:
Vec     vec;
};
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
class ColumnVector
{
typedef ColumnVector< Dim, T >  Self;
typedef Vector< Dim, T >        Vec;
typedef RowVector< Dim, T >     RowVec;
typedef Matrix< Dim, 1, T >     Mat;
public:
// types and constants
typedef Vec::Element    Element;
typedef Vec::Scalar     Scalar;
enum { dimension = Dim };
// generators
static Self     rand( void ){ return Vec::rand(); }
static ColumnVector< 0, T >     empty( void ){ return Vec::empty(); }
// constructors
ColumnVector( void ){}
ColumnVector( const Self &v ) : vec( v.vec ){}
// implicit converters
ColumnVector( const Vec &v ) : vec( v ){}
ColumnVector( const Mat &m ) : vec( transpose( m )[0] ){}
operator Vec( void ) const { return vec; }
operator Mat( void ) const { return ( Matrix< Dim, 0, T >() << *this ); }
// explicit converters
friend Self     column_vector( const Vec &v ){ return v; }
friend Self     column_vector( const Mat &m ){ return m; }
//      friend Vec      vector( const Self &v ){ return v; }

friend Mat      matrix( const Self &v ){ return v; }
// access grant operators
Element         &operator[]( size_t n ){ return vec[n]; }
const Element   &operator[]( size_t n ) const { return vec[n]; }
// unary operators
Self            operator+( void ) const { return ( +vec ); }
Self            operator-( void ) const { return ( -vec ); }
// assignment operators
Self            &operator+=( const Self &v1 ){ vec += v1.vec; return *this; }
Self            &operator-=( const Self &v1 ){ vec -= v1.vec; return *this; }
Self            &operator*=( Scalar s ){ vec *= s; return *this; }
Self            &operator/=( Scalar s ){ vec /= s; return *this; }
// binary operators
friend Self     operator+( const Self &v1, const Self &v2 ){ return ( v1.vec + v2.vec ); }
friend Self     operator-( const Self &v1, const Self &v2 ){ return ( v1.vec - v2.vec ); }
friend Self     operator*( const Self &v1, Scalar s ){ return ( v1.vec * s ); }
friend Self     operator*( Scalar s, const Self &v1 ){ return ( s * v1.vec ); }
friend Self     operator/( const Self &v1, Scalar s ){ return ( v1.vec / s ); }
// equality operators
friend bool     operator==( const Self &v1, const Self &v2 ){ return ( v1.vec == v2.vec ); }
friend bool     operator!=( const Self &v1, const Self &v2 ){ return ( v1.vec != v2.vec ); }
// vector functions
friend Scalar   inner_product( const Self &v1, const Self &v2 )
                        { return inner_product( v1.vec, v2.vec ); }
friend Scalar   sqrabs( const Self &v1 ){ return sqrabs( v1.vec ); }
friend Scalar   abs( const Self &v1 ){ return abs( v1.vec ); }
friend Self     normalize( const Self &v1 ){ return normalize( v1.vec ); }
friend RowVec   transpose( const Self &v ){ return Vec( v ); }
friend RowVec   conjugate( const Self &v ){ return conjugate( Vec( v ) ); }
// inserters
friend ColumnVector< Dim + 1, T >       operator<<( const Self &v1, Scalar s )
                                                { return ( v1.vec << s ); }
friend std::ostream     &operator<<( std::ostream &ostr, const Self &v )
                                { return ( ostr << matrix( v ) ); }
private:
Vec     vec;
};
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
class Matrix
{
typedef Matrix< Row, Col, T >   Self;
typedef RowVector< Col, T >     RowVec;
public:
// types and constants
typedef T       Element, Scalar;
enum { row = Row, column = Col };

        // generators
        static Self     identity( void );
        static Self     ones( void );
        static Self     rand( void );
        static Matrix< 0, Col > empty( void ){ return Matrix< 0, Col >(); }
        // constructors
        Matrix( void ){}
        Matrix( const Self &m );
        // unary operators
        Self            operator+( void ) const;
        Self            operator-( void ) const;
        // access grant operators
        RowVec          &operator[]( size_t n ){ return rows[n]; }
        const RowVec    &operator[]( size_t n ) const { return rows[n]; }
        // assignment operators
        Self            &operator+=( const Self &m1 );
        Self            &operator-=( const Self &m1 );
        Self            &operator*=( Scalar s );
        Self            &operator/=( Scalar s );
private:
        RowVec  rows[ ( Row > 0 ? Row : 1 ) ];
};
// binary operators
template< size_t Row, size_t Col, class T >
Matrix< Row, Col, T >       operator+( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 );
template< size_t Row, size_t Col, class T >
Matrix< Row, Col, T >       operator-( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 );
template< size_t Row, size_t ColRow, size_t Col, class T >
Matrix< Row, Col, T >        operator*( const Matrix< Row, ColRow, T > &m1,
                                        const Matrix< ColRow, Col, T > &m2 );
template< size_t Row, size_t Col, class T >
RowVector< Col, T >         operator*( const RowVector< Row, T > &v1, const Matrix< Row, Col, T > &m2 );
template< size_t Row, size_t Col, class T >
ColumnVector< Row, T >      operator*( const Matrix< Row, Col, T > &m1, const ColumnVector< Col, T > &v2 );
template< size_t Row, size_t Col, class T >
Matrix< Row, Col, T >       operator*( const Matrix< Row, Col, T > &m1, Matrix< Row, Col, T >::Scalar s );
template< size_t Row, size_t Col, class T >
Matrix< Row, Col, T >       operator*( Matrix< Row, Col, T >::Scalar s, const Matrix< Row, Col, T > &m1 );
template< size_t Row, size_t Col, class T >
Matrix< Row, Col, T >       operator/( const Matrix< Row, Col, T > &m1, Matrix< Row, Col, T >::Scalar s );
// equality operators
template< size_t Row, size_t Col, class T >
bool                        operator==( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 );
template< size_t Row, size_t Col, class T >
bool                        operator!=( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 );
// matrix functions
template< size_t Row, size_t Col, class T >
Matrix< Col, Row, T >       transpose( const Matrix< Row, Col, T > &m1 );
template< size_t Row, size_t Col, class T >
Matrix< Col, Row, T >       conjugate( const Matrix< Row, Col, T > &m1 );

template< size_t RowCol, class T >
Matrix< RowCol, RowCol, T >          pow( const Matrix< RowCol, RowCol, T > &m1, size_t e );
// inserters
template< size_t Row, size_t Col, class T >
Matrix< Row + 1, Col, T >   operator<<( const Matrix< Row, Col, T > &m1, const RowVector< Col, T > &v );
template< size_t Row, size_t Col, class T >
Matrix< Row, Col + 1, T >   operator<<( const Matrix< Row, Col, T > &m1, const ColumnVector< Row, T > &v );
template< size_t Row, size_t Col, class T >
std::ostream                &operator<<( std::ostream &ostr, const Matrix< Row, Col, T > &m1 );
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim, T > Vector< Dim, T >::rand( void )
{
        Self    v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = random< T >();
        return v;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim, T >::Vector( void )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] = Element();
}
template< size_t Dim, class T >
inline Vector< Dim, T >::Vector( const Self &v )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] = v.elements[i];
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim, T > Vector< Dim, T >::operator+( void ) const
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = +elements[i];
        return v;
}
template< size_t Dim, class T >
inline Vector< Dim, T > Vector< Dim, T >::operator-( void ) const
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = -elements[i];

        return v;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim, T > &Vector< Dim, T >::operator+=( const Self &v1 )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] += v1[i];
        return *this;
}
template< size_t Dim, class T >
inline Vector< Dim, T > &Vector< Dim, T >::operator-=( const Self &v1 )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] -= v1[i];
        return *this;
}
template< size_t Dim, class T >
inline Vector< Dim, T > &Vector< Dim, T >::operator*=( Scalar s )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] *= s;
        return *this;
}
template< size_t Dim, class T >
inline Vector< Dim, T > &Vector< Dim, T >::operator/=( Scalar s )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                elements[i] /= s;
        return *this;
}
///////////////////////////////////////////////////////////////////////////////
// Vector + Vector
template< size_t Dim, class T >
inline Vector< Dim, T > operator+( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = v1[i] + v2[i];
        return v;
}
// Vector - Vector
template< size_t Dim, class T >
inline Vector< Dim, T > operator-( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = v1[i] - v2[i];
        return v;
}

// Vector * Scalar
template< size_t Dim, class T >
inline Vector< Dim, T > operator*( const Vector< Dim, T > &v1, Vector< Dim, T >::Scalar s )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = v1[i] * s;
        return v;
}
// Scalar * Vector
template< size_t Dim, class T >
inline Vector< Dim, T > operator*( Vector< Dim, T >::Scalar s, const Vector< Dim, T > &v1 )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = s * v1[i];
        return v;
}
// Vector / Scalar
template< size_t Dim, class T >
inline Vector< Dim, T > operator/( const Vector< Dim, T > &v1, Vector< Dim, T >::Scalar s )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = v1[i] / s;
        return v;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline bool     operator==( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                if ( !( v1[i] == v2[i] ) )
                        return false;
        return true;
}
template< size_t Dim, class T >
inline bool     operator!=( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 )
{
        for ( size_t i = 0 ; i < Dim ; ++i )
                if ( v1[i] != v2[i] )
                        return true;
        return false;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim, T >::Scalar inner_product( const Vector< Dim, T > &v1, const Vector< Dim, T > &v2 )
{
        Vector< Dim, T >::Scalar        s = Vector< Dim, T >::Scalar();
        for ( size_t i = 0 ; i < Dim ; ++i )

                s += v1[i] * v2[i];
        return s;
}
template< size_t Dim, class T >
inline Vector< Dim, T >::Scalar sqrabs( const Vector< Dim, T > &v1 )
{
        return inner_product( v1, v1 );
}
template< size_t Dim, class T >
inline Vector< Dim, T >::Scalar abs( const Vector< Dim, T > &v1 )
{
        return sqrt( sqrabs( v1 ) );
}
template< size_t Dim, class T >
inline Vector< Dim, T > normalize( const Vector< Dim, T > &v1 )
{
        return ( v1 / abs( v1 ) );
}
template< size_t Dim, class T >
inline Vector< Dim, T > conjugate( const Vector< Dim, T > &v1 )
{
        Vector< Dim, T >        v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = conjugate( v1[i] );
        return v;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline Vector< Dim + 1, T >     operator<<( const Vector< Dim, T > &v1, Vector< Dim, T >::Element e )
{
        Vector< Dim + 1, T >    v;
        for ( size_t i = 0 ; i < Dim ; ++i )
                v[i] = v1[i];
        v[Dim] = e;
        return v;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Dim, class T >
inline std::ostream     &operator<<( std::ostream &ostr, const Vector< Dim, T > &v )
{
        return ( ostr << RowVector< Dim, T >( v ) );
}
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    Matrix< Row, Col, T >::identity( void )

{
        const int       n = ( Row < Col ? Row : Col );
        Self    m;
        for ( size_t i = 0 ; i < n ; ++i )
                m[i][i] = T( 1 );
        return m;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    Matrix< Row, Col, T >::ones( void )
{
        Self    m;
        for ( size_t i = 0 ; i < Row ; ++i )
                for ( size_t j = 0 ; j < Col ; ++j )
                        m[i][j] = T( 1 );
        return m;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    Matrix< Row, Col, T >::rand( void )
{
        Self    m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = RowVec::rand();
        return m;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >::Matrix( const Matrix< Row, Col, T > &m )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                rows[i] = m.rows[i];
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    Matrix< Row, Col, T >::operator+( void ) const
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = +rows[i];
        return m;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    Matrix< Row, Col, T >::operator-( void ) const
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = -rows[i];
        return m;
}
///////////////////////////////////////////////////////////////////////////////

template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    &Matrix< Row, Col, T >::operator+=( const Self &m1 )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                rows[i] += m1[i];
        return *this;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    &Matrix< Row, Col, T >::operator-=( const Self &m1 )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                rows[i] -= m1[i];
        return *this;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    &Matrix< Row, Col, T >::operator*=( Scalar s )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                rows[i] *= s;
        return *this;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    &Matrix< Row, Col, T >::operator/=( Scalar s )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                rows[i] /= s;
        return *this;
}
///////////////////////////////////////////////////////////////////////////////
// Matrix + Matrix
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    operator+( const Matrix< Row, Col, T > &m1,
                                           const Matrix< Row, Col, T > &m2 )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = m1[i] + m2[i];
        return m;
}
// Matrix - Matrix
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    operator-( const Matrix< Row, Col, T > &m1,
                                           const Matrix< Row, Col, T > &m2 )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = m1[i] - m2[i];
        return m;
}
// Matrix * Matrix
template< size_t Row, size_t ColRow, size_t Col, class T >

inline Matrix< Row, Col, T >    operator*( const Matrix< Row, ColRow, T > &m1,
                                           const Matrix< ColRow, Col, T > &m2 )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                for ( size_t j = 0 ; j < Col ; ++j )
                        for ( size_t k = 0 ; k < ColRow ; ++k )
                                m[i][j] += m1[i][k] * m2[k][j];
        return m;
}
// Vector * Matrix
template< size_t Row, size_t Col, class T >
inline RowVector< Col, T >      operator*( const RowVector< Row, T > &v1,
                                           const Matrix< Row, Col, T > &m2 )
{
        return ( matrix( v1 ) * m2 );
}
// Matrix * Vector
template< size_t Row, size_t Col, class T >
inline ColumnVector< Row, T >   operator*( const Matrix< Row, Col, T > &m1,
                                           const ColumnVector< Col, T > &v2 )
{
        return ( m1 * matrix( v2 ) );
}
// Matrix * Scalar
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    operator*( const Matrix< Row, Col, T > &m1,
                                           Matrix< Row, Col, T >::Scalar s )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = m1[i] * s;
        return m;
}
// Scalar * Matrix
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    operator*( Matrix< Row, Col, T >::Scalar s,
                                           const Matrix< Row, Col, T > &m1 )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = s * m1[i];
        return m;
}
// Matrix / Scalar
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col, T >    operator/( const Matrix< Row, Col, T > &m1,
                                           Matrix< Row, Col, T >::Scalar s )
{
        Matrix< Row, Col, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = m1[i] / s;
        return m;

}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline bool     operator==( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                if ( !( m1[i] == m2[i] ) )
                        return false;
        return true;
}
template< size_t Row, size_t Col, class T >
inline bool     operator!=( const Matrix< Row, Col, T > &m1, const Matrix< Row, Col, T > &m2 )
{
        for ( size_t i = 0 ; i < Row ; ++i )
                if ( m1[i] != m2[i] )
                        return true;
        return false;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline Matrix< Col, Row, T >    transpose( const Matrix< Row, Col, T > &m1 )
{
        Matrix< Col, Row, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                for ( size_t j = 0 ; j < Col ; ++j )
                        m[j][i] = m1[i][j];
        return m;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Col, Row, T >    conjugate( const Matrix< Row, Col, T > &m1 )
{
        Matrix< Col, Row, T >   m;
        for ( size_t i = 0 ; i < Row ; ++i )
                for ( size_t j = 0 ; j < Col ; ++j )
                        m[j][i] = conjugate( m1[i][j] );
        return m;
}
template< size_t RowCol, class T >
inline Matrix< RowCol, RowCol, T >      pow( const Matrix< RowCol, RowCol, T > &m1, size_t e )
{
        Matrix< RowCol, RowCol, T >     m, t;
        t = m1;
        m = ( e & 1 ? m1 : Matrix< RowCol, RowCol, T >::identity() );
        for ( e >>= 1 ; e ; e >>= 1 )
        {
                t = t * t;
                if ( e & 1 )
                        m = m * t;
        }

        return m;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline Matrix< Row + 1, Col, T >        operator<<( const Matrix< Row, Col, T > &m1,
                                                    const RowVector< Col, T > &v )
{
        Matrix< Row + 1, Col, T >       m;
        for ( size_t i = 0 ; i < Row ; ++i )
                m[i] = m1[i];
        m[Row] = v;
        return m;
}
template< size_t Row, size_t Col, class T >
inline Matrix< Row, Col + 1, T >        operator<<( const Matrix< Row, Col, T > &m1,
                                                    const ColumnVector< Row, T > &v )
{
        Matrix< Row, Col + 1, T >       m;
        for ( size_t i = 0 ; i < Row ; ++i )
        {
                for ( size_t j = 0 ; j < Col ; ++j )
                        m[i][j] = m1[i][j];
                m[i][Col] = v[i];
        }
        return m;
}
///////////////////////////////////////////////////////////////////////////////
template< size_t Row, size_t Col, class T >
inline std::ostream     &operator<<( std::ostream &ostr, const Matrix< Row, Col, T > &m1 )
{
        for ( size_t i = 0 ; i < Row ; ++i )
        {
                for ( size_t j = 0 ; j < Col ; ++j )
                        ostr << " " << std::setw( 6 + 6 ) << std::setprecision( 6 )
                             << std::setiosflags( std::ios::showpoint ) << m1[i][j];
                ostr << std::endl;
        }
        return ostr;
}
///////////////////////////////////////////////////////////////////////////////

A.2 ガウスの消去法

A.2.1 ピボット選択無し

template< size_t Dim, class T >

inline void     gaussian_elimination( Vector< Dim, T > &result,
                      Matrix< Dim, Dim, T > a, Vector< Dim, T > b )
{
for ( int k = 0 ; k < Dim - 1 ; ++k )
{
for ( int i = k + 1 ; i < Dim ; ++i )
{
        const T r = a[i][k] / a[k][k];
        for ( int j = k + 1 ; j < Dim ; ++j )
                a[i][j] -= r * a[k][j];
        b[i] -= r * b[k];
}
}
for ( int i = Dim - 1 ; i >= 0 ; --i )
{
T       t = T();
for ( int j = i + 1 ; j < Dim ; ++j )
        t += a[i][j] * b[j];
b[i] = ( b[i] - t ) / a[i][i];
}
result = b;
}

A.2.2 ピボット選択有り

template< size_t Dim, class T >
inline bool     gaussian_elimination_pivot( Vector< Dim, T > &result,
                            Matrix< Dim, Dim, T > a, Vector< Dim, T > b )
{
for ( int k = 0 ; k < Dim ; ++k )
{
int     p = k;
const T e = ldexp( 1.0, -50 );
// pivot
for ( int i = k + 1 ; i < Dim ; ++i )
        if ( fabs( a[p][k] ) < fabs( a[i][k] ) )
                p = i;
if ( fabs( a[p][k] ) < e )
        return false;
if ( p != k )
{
        std::swap( a[k], a[p] );
        std::swap( b[k], b[p] );
}
for ( int i = k + 1 ; i < Dim ; ++i )
{
        const T r = a[i][k] / a[k][k];
        for ( int j = k + 1 ; j < Dim ; ++j )
                a[i][j] -= r * a[k][j];
        b[i] -= r * b[k];
}
}

for ( int i = Dim - 1 ; i >= 0 ; --i )
{
T       t = T();
for ( int j = i + 1 ; j < Dim ; ++j )
        t += a[i][j] * b[j];
b[i] = ( b[i] - t ) / a[i][i];
}
result = b;
return true;
}

A.3 LU 分解法

A.3.1 LU 分解

template< size_t Dim, class T >
inline bool     lu_decomposition( Matrix< Dim, Dim, T > &lower, Matrix< Dim, Dim, T > &upper,
                  Vector< Dim, size_t > &pivot, Matrix< Dim, Dim, T > a )
{
for ( int k = 0 ; k < Dim ; ++k )
{
const T e = ldexp( 1.0, -50 );
int     p = k;
for ( int i = k + 1 ; i < Dim ; ++i )
        if ( fabs( a[p][k] ) < fabs( a[i][k] ) )
                p = i;
if ( fabs( a[p][k] ) < e )
        return false;
pivot[k] = p;
if ( p != k )
        std::swap( a[k], a[p] );
for ( int i = k + 1 ; i < Dim ; ++i )
{
        const T r = a[i][k] / a[k][k];
        a[i][k] = r;
        for ( int j = k + 1 ; j < Dim ; ++j )
                a[i][j] -= r * a[k][j];
}
}
for ( int i = 0 ; i < Dim ; ++i )
{
for ( int j = 0 ; j < i ; ++j )
{
        lower[i][j] = a[i][j];
        upper[i][j] = 0.0;
}
lower[i][i] = 1.0;
upper[i][i] = a[i][i];

for ( int j = i + 1 ; j < Dim ; ++j )
{
        lower[i][j] = 0.0;
        upper[i][j] = a[i][j];
}
}
return true;
}

A.3.2 LU 分解法

template< size_t Dim, class T >
inline void     lu_solve( Vector< Dim, T > &x, const Matrix< Dim, Dim, T > &lower,
          const Matrix< Dim, Dim, T > &upper, const Vector< Dim, size_t > &pivot,
          Vector< Dim, T > b )
{
for ( int i = 0 ; i < Dim ; ++i )
std::swap( b[i], b[pivot[i]] );
for ( int k = 0 ; k < Dim - 1 ; ++k )
for ( int i = k + 1 ; i < Dim ; ++i )
        b[i] -= lower[i][k] * b[k];
for ( int i = Dim - 1 ; i >= 0 ; --i )
{
T       t = T();
for ( int j = i + 1 ; j < Dim ; ++j )
        t += upper[i][j] * b[j];
b[i] = ( b[i] - t ) / upper[i][i];
}
x = b;
}

A.3.3 逆行列

template< size_t Dim, class T >
inline bool     inverse( Matrix< Dim, Dim, T > &x, const Matrix< Dim, Dim, T > &a )
{
Matrix< Dim, Dim, T >   lower, upper;
Vector< Dim, size_t >   pivot;
if ( !lu_decomposition( lower, upper, pivot, a ) )
return false;
for ( int i = 0 ; i < Dim ; ++i )
{
Vector< Dim, T >        b, xi;
b[i] = 1.0;
lu_solve( xi, lower, upper, pivot, b );
for ( int j = 0 ; j < Dim ; ++j )
        x[j][i] = xi[j];

}
return true;
}

A.4 区間演算

template< class T >
class Interval
{
typedef Interval< T >   Self;
static const unsigned short     MASK = 0x0c00, UP = 0x0800, DOWN = 0x0400;
public:
T       lower, upper;
Interval( void ) : lower( T() ), upper( T() ){}
Interval( T v ) : lower( v ), upper( v ){}
Interval( T l, T u ) : lower( l ), upper( u ){}
bool    contains( T v ) const { return ( lower <= v && v <= upper ); }

static void     up( void )
{
        unsigned short  cw;
        asm( "fstcw %0" : "=g" (cw) );
        cw &= ~MASK;
        cw |= UP;
        asm( "fldcw %0" : : "g" (cw) );
}
static void     down( void )
{
        unsigned short  cw;
        asm( "fstcw %0" : "=g" (cw) );
        cw &= ~MASK;
        cw |= DOWN;
        asm( "fldcw %0" : : "g" (cw) );
}
};
////////////////////////////////////////////////////////////////////////////////
template< class T >
inline Interval< T >    operator+( const Interval< T > &x, const Interval< T > &y )
{
Interval< T >   z;
z.down();
z.lower = x.lower + y.lower;
z.up();
z.upper = x.upper + y.upper;

return z;
}
template< class T >
inline Interval< T >    operator-( const Interval< T > &x, const Interval< T > &y )
{
Interval< T >   z;
z.down();
z.lower = x.lower - y.upper;
z.up();
z.upper = x.upper - y.lower;
return z;
}
template< class T >
inline Interval< T >    operator*( const Interval< T > &x, const Interval< T > &y )
{
T       lower[4], upper[4];
Interval< T >   z;
z.down();
lower[0] = x.lower * y.lower;
lower[1] = x.lower * y.upper;
lower[2] = x.upper * y.lower;
lower[3] = x.upper * y.upper;
z.lower = *std::min_element( lower, lower + 4 );
z.up();
upper[0] = x.lower * y.lower;
upper[1] = x.lower * y.upper;
upper[2] = x.upper * y.lower;
upper[3] = x.upper * y.upper;
z.upper = *std::max_element( upper, upper + 4 );
return z;
}
template< class T >
inline Interval< T >    operator/( const Interval< T > &x, const Interval< T > &y )
{
T       lower[4], upper[4];
Interval< T >   z;
if ( y.contains( 0.0 ) )
        z = Interval< T >( -1.0 / 0.0, 1.0 / 0.0 );
else
{
        z.down();
        lower[0] = x.lower / y.lower;
        lower[1] = x.lower / y.upper;
        lower[2] = x.upper / y.lower;
        lower[3] = x.upper / y.upper;
        z.lower = *std::min_element( lower, lower + 4 );

                z.up();
                upper[0] = x.lower / y.lower;
                upper[1] = x.lower / y.upper;
                upper[2] = x.upper / y.lower;
                upper[3] = x.upper / y.upper;
                z.upper = *std::max_element( upper, upper + 4 );
        }
        return z;
}
////////////////////////////////////////////////////////////////////////////////
template< class T >
inline Interval< T >    &operator+=( Interval< T > &x, const Interval< T > &y )
{
        x = x + y;
        return x;
}
template< class T >
inline Interval< T >    &operator-=( Interval< T > &x, const Interval< T > &y )
{
        x = x - y;
        return x;
}
template< class T >
inline Interval< T >    &operator*=( Interval< T > &x, const Interval< T > &y )
{
        x = x * y;
        return x;
}
template< class T >
inline Interval< T >    &operator/=( Interval< T > &x, const Interval< T > &y )
{
        x = x / y;
        return x;
}
////////////////////////////////////////////////////////////////////////////////
template< class T >
inline bool     operator<( const Interval< T > &x, const Interval< T > &y )
{
        return ( x.upper < y.lower );
}
template< class T >
inline bool     operator<=( const Interval< T > &x, const Interval< T > &y )
{
        return ( x.upper <= y.lower );
}
template< class T >
inline bool     operator>( const Interval< T > &x, const Interval< T > &y )
{

        return ( x.lower > y.upper );
}
template< class T >
inline bool     operator>=( const Interval< T > &x, const Interval< T > &y )
{
        return ( x.lower >= y.upper );
}
////////////////////////////////////////////////////////////////////////////////
template< class T >
inline T        minabs( const Interval< T > &x )
{
        if ( x.contains( 0.0 ) )
                return 0.0;
        else
                return std::min( fabs( x.lower ), fabs( x.upper ) );
}
template< class T >
inline T        maxabs( const Interval< T > &x )
{
        return std::max( fabs( x.lower ), fabs( x.upper ) );
}
template< class T >
inline Interval< T >    fabs( const Interval< T > &x )
{
        return Interval< T >( minabs( x ), maxabs( x ) );
}
////////////////////////////////////////////////////////////////////////////////
template< class T >
inline std::ostream     &operator<<( std::ostream &ostr, const Interval< T > &x )
{
        ostr    << std::setw( 0 ) << "[ "
                << std::setw( 6 + 17 ) << std::setprecision( 17 )
                << std::setiosflags( std::ios::showpoint ) << x.lower
                << " " << std::setw( 6 + 17 ) << std::setprecision( 17 )
                << std::setiosflags( std::ios::showpoint ) << x.upper << " ]";
        return ostr;
}
////////////////////////////////////////////////////////////////////////////////

A.5 メインプログラム

void    test( void )
{
        const size_t    Dim = 5;
        typedef Matrix< Dim, Dim >      Mat;
        typedef Vector< Dim >           Vec;
        Mat     a, l, u, p, inv;
        Vec     x, b;

Vector< Dim, size_t >   pivot;
Matrix< Dim, Dim, Interval< double > >  ia, ix;

a = Mat::rand();
a[0][0] = 0.0;
b = Vec::rand();
cout    << "a =" << endl << a << endl
        << "b =" << endl << b << endl;

cout << endl;
cout << "gaussian elimination (w/o pivot selection)" << endl << endl;
gaussian_elimination( x, a, b );
cout << "a\\b =" << endl << x << endl;
cout << endl;
cout << "gaussian elimination (with pivot selection)" << endl << endl;
if ( gaussian_elimination_pivot( x, a, b ) )
        cout << "a\\b =" << endl << x << endl;
else
        cout << "failure." << endl << endl;

cout << endl;
cout << "LU decomposition" << endl << endl;
lu_decomposition( l, u, pivot, a );
p = Mat::identity();
for ( int i = 0 ; i < Dim ; ++i )
        std::swap( p[i], p[pivot[i]] );
cout
        << "l =" << endl << l << endl
        << "u =" << endl << u << endl
        << "pivot =" << endl << pivot << endl
        << "p =" << endl << p << endl
        << "p^-1 * l * u =" << endl << transpose( p ) * l * u << endl
        << endl;
lu_solve( x, l, u, pivot, b );
cout
        << "x = a\\b =" << endl << x << endl
        << endl;
inverse( inv, a );
cout    << "a^-1 =" << endl << inv << endl
        << "a * a^-1 =" << endl << a * inv << endl
        << endl;
for ( int i = 0 ; i < Dim ; ++i )
        for ( int j = 0 ; j < Dim ; ++j )
                ia[i][j] = a[i][j];
inverse( ix, ia );
cout    << "a * a^-1 =" << endl << ia * ix << endl
        << endl;

}
////////////////////////////////////////////////////////////////////////////////
void    solve( void )
{
Matrix< 3, 3 >  a, l, u;
Vector< 3 >     b1, b2, b3, x1, x2, x3;
Vector< 3, size_t >     p;
cout << endl;
cout << "solutions" << endl << endl;
a[0][0] = 1.0;  a[0][1] = 2.0;  a[0][2] = 3.0;
a[1][0] = 2.0;  a[1][1] = 8.0;  a[1][2] = 11.0;
a[2][0] = 3.0;  a[2][1] = 22.0; a[2][2] = 35.0;
b1[0] = 1.0;    b1[1] = 0.0;    b1[2] = 0.0;
b2[0] = 1.0;    b2[1] = 2.0;    b2[2] = 3.0;
b3[0] = -1.0;   b3[1] = 0.0;    b3[2] = 1.0;
lu_decomposition( l, u, p, a );
cout << l << endl << u << endl;
lu_solve( x1, l, u, p, b1 );
cout << x1 << endl;
lu_solve( x2, l, u, p, b2 );
cout << x2 << endl;
lu_solve( x3, l, u, p, b3 );
cout << x3 << endl;
}
////////////////////////////////////////////////////////////////////////////////
int     main( void )
{
srand( time( NULL ) );
test();
solve();
return 0;
}
////////////////////////////////////////////////////////////////////////////////