2 実際のプログラム(手取り足取り)

ここでは,実際のプログラムを作成するときの考え方を示す.最初に, 何にも考えていないガウス・ジョルダン法から出発し,少しずつ機能を追加し て,最終的にパッケージとして完成した関数を作成する事にする.以下の順序でプロ グラムをブラッシュアップする.
  1. $ \boldsymbol{A}$を単位行列にする.そして, $ \boldsymbol{x}$を求める.
  2. 行を交換するだけのピボット選択(部分選択)の機能を追加する.
  3. 逆行列を計算するルーチンを追加する.
  4. 不要な配列を排除して,メモリー効率を上げる改造をする.


2.1 素朴なガウス・ジョルダン法(行列の対角化のみ)

まずは,行列の対角化(単位行列に変換)のみのプログラムを作成する.ピボット選択や逆 行列は考えない.係数行列 $ \boldsymbol{A}$を単位行列に,非同次項 $ \boldsymbol{b}$は解ベ クトル $ \boldsymbol{x}$に変換することのみを取り扱う.

この節のプログラムは,ピボット選択がないため,実用上問題を含んでいる.対角成分に ゼロが現れた場合,計算ができなくなる.さらに,行列が特異な場合でも,同様なことが 生じる.このようなとき,ゼロで割ることになるので,実行時エラーが発生する.あるい は,大きな計算誤差を伴った解になる.しかし,最初の学習では,これは気にしない-- エラーの処理のルーチンを書かない--ことにする.ガウス・ジョルダン法のプログラム の学習は簡単な方が良い.しかし,諸君が学習ではなく実際に使うプログラムを組むとき, ピボット選択は必要不可欠である--ということを忘れてはならない.

2.1.1 最外殻のループ

$ N \times N$の行列 $ \boldsymbol{A}$を対角化することを考える.この場合, $ a_{11},\,a_{22},\,a_{33},\,\cdots,\,a_{NN}$と同じ手順で対角化を進めることになる. ということは,$ N$回のループが必要である.プログラムでは,以下のループ構造が一番 外側で回ることになる.ループの回数が分かっている場合,for文を使うのが普通である.
	for(ipv=1 ; ipv <= n ; ipv++){

	   対角化の処理

	}
ここで,ipvは対角化する要素 $ a_{ipv\,ipv}$の添え字を表す. nは行列のサイズである.

2.1.2 対角成分を1に(ピボット行の処理)

次の処理は,対角成分を $ a_{ipv\,ipv}=1$にすることである. $ a_{ipv-1\,ipv-1}$間では 対角化できており,次の成分を対角化するということです.もし,$ ipv=1$ならば最初の 対角成分を1にすることになる.最初であろうが,途中であろうが同じようにプログラム は書かなくてはならない.具体的には,

$\displaystyle \boldsymbol{A} = \begin{pmatrix}1 & 0 & \ast & \ast & \ast & \ast...
...\ast & \ast & \ast & \ast \\ 0 & 0 & \ast & \ast & \ast & \ast \end{pmatrix}\\ $    

と変形したいのである.係数行列 $ \boldsymbol{A}$を変形させれば,同じ操作により非同次項 $ \boldsymbol{b}$も処理しなくてならない.

数学では,対角成分を1にするために,その行を対角成分で割る.しかし,コンピューター のプログラムでは予め逆数を計算して,それを乗じた方が良い.コンピューターは除算よ りも乗算の方が得意なのでできるだけ,除算は使わないようにプログラムを記述すべきで ある.非同次項 $ \boldsymbol{b}$の演算は1回であるが,係数行列 $ \boldsymbol{A}$は列毎なので$ N$回の演 算が必要になる.対角成分を1にする処理は,次のようにする.

	inv_pivot = 1.0/a[ipv][ipv];

	for(j=1 ; j <= n ; j++){
	   a[ipv][j] *= inv_pivot;
	}

	b[ipv] *= inv_pivot;

これでピボットのある行の処理は終わり.

2.1.3 ピボットのある列を0に(ピボット行以外の行の処理)

次は,ピボットがある行以外の処理である.それは,ピボットがある列を全てゼ ロにすることに他ならない.要するに次のように,係数行列 $ \boldsymbol{A}$を変形する.

$\displaystyle \boldsymbol{A} = \begin{pmatrix}1 & 0 & a_{1\,ipv} & \ast & \ast ...
...} \\ 0 & 0 & 0 & \ast^{\prime} & \ast^{\prime} & \ast^{\prime} \\ \end{pmatrix}$    

このように係数行列 $ \boldsymbol{A}$を変形させる.もちろん,方程式の解 $ \boldsymbol{x}$ が変わらないように,非同次項 $ \boldsymbol{b}$も同じ操作をする.

このように変形するのは簡単である.例えば,$ i$行を処理する場合を考える. $ i$行を,ピボットのある$ ipv$行を $ a_{i\,ipv}$倍したもので引けば良いのである.

\begin{displaymath}\begin{array}{@{\,}lcccccccccccl@{\,}} \text{$i$行} & \Righta...
...qquad & \multicolumn{2}{r}{b_i-a_{i\,ipv}b_{ipv}}\\ \end{array}\end{displaymath}    

この処理の実際のプログラムは次のようになる.これで, $ i=1,2,3,\cdots,n$ $ j=1,2,3,\cdots,n$の全ての $ \boldsymbol{A}$の成分を処理をする.

	for(i=1 ; i<=n ; i++){
	   if(i != ipv){
	      temp = a[i][ipv];
	      for(j=1 ; j<=n ; j++){
	         a[i][j] -= temp*a[ipv][j];
	      }
	      b[i] -= temp*b[ipv];
	   }
	}

これで対角化の処理はおしまい.

2.1.4 素朴なガウス・ジョルダン法のソースプログラム

以上をまとめると,ピボット選択は無く逆行列も求めないのガウス・ジョルダ ン法が完成する.ここで,ひとつ気が付いてほしいのは解ベクトル $ \boldsymbol{x}$ のための配列は不要ということである.非同次項の配列が解になっている.従っ て,この最も素朴なガウス・ジョルダン法のプログラムは次のようになる.

	/* ========== ガウスジョルダン法の関数 =================*/
	void gauss_jordan(int n, double a[][100], double b[])
	{

	   int ipv, i, j;
	   double inv_pivot, temp;

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 対角成分=1(ピボット行の処理) ---- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }
	}


2.2 ピボット選択機能追加(行交換)

先ほどの素朴なガウス・ジョルダン法は,爆弾を抱えた関数になっている.もし,対角成 分にゼロが現れたら,ゼロで割ることになり処理が破綻するのである.そこで,それを解 決するピボット選択が登場するのである.もっとも,この問題の解決ばかりでなく,解の 精度も向上する.

ピボット選択,ここでは行の交換のみの部分選択を考える.その処理は,

の2つの部分から構成される.これらは,先のプログラムの対角成分を1にする処理の前に 入れなくてはならない.

2.2.1 最大値探索

行交換のみを行う部分選択の場合,ピボットはピボット行以下の最大値とする.これは, 今まで処理した,対角成分が1になっている部分を崩さないためである.

\begin{displaymath}\begin{array}{rrcccccc@{\hspace{5pt}}l} \ldelim\{{2}{59pt}[行...
... \ast & \\ & &0 &0 & \ast & \ast & \ast & \ast & \\ \end{array}\end{displaymath}    

最大値は,$ ipv$行よりも下の行で,$ ipv$列の最大値は,以降に示すプログラ ムにより探すことができる.最大値は,fabsという関数で絶対値 を比較することにより求める.

ここで,もし最大値がゼロの場合,行列は特異(行列式がゼロ)ということにな り,解は一意的に決まらない.その場合,関数の値として1を返し,その ことをコールした側に伝えるのが良い.

	big=0.0;
	for(i=ipv ; i<=n ; i++){
	   if(fabs(a[i][ipv]) > big){
	      big = fabs(a[i][ipv]);
	      pivot_row = i;
	   }
	}
	if(big == 0.0) return 1;
	row[ipv] = pivot_row;

このプログラムは,以下のことを行っている.

2.2.2 行の交換

ピボットとすべき値がある行(pivot_row)がわかったので,$ ipv$行 と入れ替える.

$\displaystyle \boldsymbol{A} = \begin{pmatrix}1 & 0 & \ast & \ast & \ast & \ast...
...ircledcirc & \circledcirc \\ 0 & 0 & \ast & \ast & \ast & \ast \\ \end{pmatrix}$    

入れ替えは,係数行列 $ \boldsymbol{A}$と非同時項 $ \boldsymbol{b}$の両方について行う. 変数を入れ替える場合,一時的な変数を記憶する場所が必要で,ここでは, tempという変数を使っている.

ただし,もともと最大の値が$ ipv$行にある場合は,行の入れ替えは行わない.

	if(ipv != pivot_row){

	   for(i=1 ; i<=n ; i++){
	      temp = a[ipv][i];
	      a[ipv][i] = a[pivot_row][i];
	      a[pivot_row][i] = temp;
	   }

	   temp = b[ipv];
	   b[ipv] = b[pivot_row];
	   b[pivot_row] = temp;
	
	}


2.2.3 部分ピボット選択付きガウス・ジョルダン法のソースプログ ラム

この2.2節をまとめ, 2.1節の素朴なガウス・ジョルダン法とあわせ ると,部分ピボット選択付きのガウス・ジョルダン法が完成する.

	/* ========= ガウスジョルダン法の関数====================== */

	int gauss_jordan(int n, double a[][MAXN+10], double b[])
	{

	   int ipv, i, j;
	   double inv_pivot, temp;
	   double big;
	   int pivot_row, row[MAXN+10];

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 最大値探索 ---------------------------- */
	      big=0.0;
	      for(i=ipv ; i<=n ; i++){
	         if(fabs(a[i][ipv]) > big){
	            big = fabs(a[i][ipv]);
	            pivot_row = i;
	         }
	      }
	      if(big == 0.0) return 1;
	      row[ipv] = pivot_row;

	      /* ---- 行の入れ替え -------------------------- */
	      if(ipv != pivot_row){
	         for(i=1 ; i<=n ; i++){
	            temp = a[ipv][i];
	            a[ipv][i] = a[pivot_row][i];
	            a[pivot_row][i] = temp;
	         }
	         temp = b[ipv];
	         b[ipv] = b[pivot_row];
	         b[pivot_row] = temp;
	      }

	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }

	   return 0;
	}

2.3 逆行列計算ルーチンの追加

逆行列を計算するルーチンは難しそうではあるが,実は単純である.単位行 列を,係数行列 $ \boldsymbol{A}$と同じ処理をすれば,それは求められる.線形代数の授業を思い 出せ.係数行列 $ \boldsymbol{A}$が単位行列に変形されたならば,元の単位行列は逆行列 に変換されるのである.これは,簡単に実現できる.

\begin{equation*}\begin{aligned}\boldsymbol{A} = \begin{pmatrix}\ast & \ast & \a...
...ar & \star & \star & \star & \star \\ \end{pmatrix} \end{aligned}\end{equation*}

これを実現するためには,以下の2つのことをすればよい.

2.3.1 単位行列の作成

まず,単位行列を作成する.これは,以下のようにすればよいであろう.

	for(i=1 ; i<=n ; i++){
	   for(j=1 ; j<=n ; j++){
	      if(i == j){
	         inv_a[i][j]=1.0;
	      }else{
	         inv_a[i][j]=0.0;
	      }
	   }
	}

2.3.2 逆行列の計算

単位行列inv_aができたので,これを係数行列 $ \boldsymbol{A}$と同じ処理 をすれば,逆行列に変換でききる.具体的には, 2.2.3節のプログラムを以下のように書き換える.


2.3.3 逆行列計算付きガウス・ジョルダン法のソースプログラム

いままで述べたことを全て網羅したガウス・ジョルダン法の計算の関数は,次 のようになる.
	/* ========= ガウスジョルダン法の関数====================== */

	int gauss_jordan(int n, double a[][MAXN+10], double b[],
	                 double inv_a[][MAXN+10]){

	   int ipv, i, j;
	   double inv_pivot, temp;
	   double big;
	   int pivot_row, row[MAXN+10];

	   /* ---- 単位行列作成 ---------------------------- */
	   for(i=1 ; i<=n ; i++){
	      for(j=1 ; j<=n ; j++){
	         if(i==j){
	            inv_a[i][j]=1.0;
	         }else{
	            inv_a[i][j]=0.0;
	         }
	      }
	   }

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 最大値探索 ---------------------------- */
	      big=0.0;
	      for(i=ipv ; i<=n ; i++){
	         if(fabs(a[i][ipv]) > big){
	            big = fabs(a[i][ipv]);
	            pivot_row = i;
	         }
	      }
	      if(big == 0.0) return 1;
	      row[ipv] = pivot_row;

	      /* ---- 行の入れ替え -------------------------- */
	      if(ipv != pivot_row){
	         for(i=1 ; i<=n ; i++){
	            temp = a[ipv][i];
	            a[ipv][i] = a[pivot_row][i];
	            a[pivot_row][i] = temp;          	
	            temp = inv_a[ipv][i];
	            inv_a[ipv][i] = inv_a[pivot_row][i];
	            inv_a[pivot_row][i] = temp;
	         }
	         temp = b[ipv];
	         b[ipv] = b[pivot_row];
	         b[pivot_row] = temp;
	      }

	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	         inv_a[ipv][j] *= inv_pivot;          
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	               inv_a[i][j] -= temp*inv_a[ipv][j];                	
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }

	   return 0;
	}

2.4 メモリー,計算効率の改善

昔,といってもそんなに過去のことではない.プログラマーは出来るだけメモリーを大事 に使った.使えるメモリーが限られていたので,その資源を有効に活用しなくてはならな かった.いまでこそ,パソコンで1G Byteのメモリーを使うのは何でもないが,たった10 年ほど前では状況は異なっていた.当時,メインフレームと言われた大型のコンピューター でさえ,1つのプログラムが使える領域は10M Byte程度であった.

メモリーと合わせて,計算効率も重要であった.大規模な計算になると,計算が終了する まで何日も費やす場合がある.そのような場合,プログラムの改良により,速度が10%アッ プするとかなりのメリットがあるのである.

そこで,ここではメモリーの効率的な利用を考える.ただし,ここは少し難しい.

2.4.1 メモリーと計算の効率化

これまでの計算過程を考える.$ i$$ ipv$列までの処理が完了したとき,係数行列 $ \boldsymbol{A}$を示す配列A[i][j]と逆行列 $ \boldsymbol{A}^{-1^\prime}$ が最終的に格納さ れる配列inv_a[i][j]の状態を見る.それぞれは,

\begin{equation*}\begin{aligned}\boldsymbol{A}^\prime = \begin{pmatrix}1 & 0 & 0...
...dast & \circledast & 0 & 0 & 0 & 1 \\ \end{pmatrix} \end{aligned}\end{equation*}

となっているはずである.この状態は,2.3.3節の i=4,ipv=3が完了したときである.この行列を良く見ると,係数行列 $ \boldsymbol{A}^\prime$では$ i$$ ipv$列までは,役に立つ情報をもっていないことが分かる 3.同様に, $ \boldsymbol{A}^{-1^\prime}$行列は, $ \boldsymbol{A}^\prime$では$ i$$ ipv$列以降 は役に立つ情報が無いことが分かる.これらの情報として役に立たない成分$ {0,1}$は, メモリーの無駄遣いなので,

\begin{equation*}\begin{aligned}\boldsymbol{A}^\prime = \begin{pmatrix}\circleda...
...ledast & \ast & \ast & \ast & \ast \\ \end{pmatrix} \end{aligned}\end{equation*}

としたくなる.そうすると,メモリーが半分で済む.これは,$ n=10000$の行列とすると, 800M Byteの節約になる.

これを実現するのは,簡単である.次のようにプログラムを書けばよい.

	/* ---- 対角成分=1(ピボット行の処理) ---------- */
	inv_pivot = 1.0/a[ipv][ipv];
	a[ipv][ipv]=1.0;                /* --- この行を追加 --- */
	for(j=1 ; j <= n ; j++){
	   a[ipv][j] *= inv_pivot;
	}
	b[ipv] *= inv_pivot;

	/* ---- ピボット列=0(ピボット行以外の処理) ---- */
	for(i=1 ; i<=n ; i++){
	   if(i != ipv){
	      temp = a[i][ipv];
	      a[i][ipv]=0.0;	        /* --- この行を追加 --- */
	      for(j=1 ; j<=n ; j++){
	         a[i][j] -= temp*a[ipv][j];
	      }
	      b[i] -= temp*b[ipv];
	   }
	}

このようにすると必要なメモリーは,2.3.3節のプログラ ムに比べて半分になる.さらに,計算時間も半分になる. 2.3.3節では,計算結果が0$ 1$の場合も計算し ていたが,このプログラムではそれを省いている.

2.4.2 逆行列の列の入れ替え

これで,メモリーと計算の効率化が図れたが,このままでは $ \boldsymbol{A}^{-1}$を示すinv_a[i][j]の配列は, $ \boldsymbol{A}$の逆行列になっ ていない.ピボット選択により,行を入れ替えた行列 $ \boldsymbol{A}^\prime$の逆 行列になっている.そこで,元の行列 $ \boldsymbol{A}$の逆行列にするために, $ \boldsymbol{A}^{-1^\prime}$の列を入れ替える必要がある.なぜ,列を入れ替え る必要があるか?.それは,以下のように考える.

当然, $ \boldsymbol{A}$の逆行列ということは,乗算すると下のように単位行列になる.

$\displaystyle \left( \begin{array}{@{\,}ccccccc@{\,}} \ast & \ast & \ast & \ast...
...ap{\smash{\Huge$0$}}\quad} & & & & \ddots & \\ & & & & & &1 \end{array} \right)$    

すなわち,元の行列の行ベクトルと逆行列の列ベクトルの内積が,

(元の行列の$ i$行の行ベクトル)$\displaystyle \cdot$(逆行列の$ j$列の列ベクトル)$\displaystyle = \delta_{ij}$    

という関係を満たしていることを意味する. $ \delta_{ij}$は,クロネッカーの記号で, $ i=j$のときその値は$ 1$$ i\neq j$のときその値は0という意味である.

このようなわけで,元の行列の行を入れ替えた場合,その逆行列は元の行列の逆行列の列 を入れ替えたものになる.従って,ピボット選択により係数行列の行を入れ替えると,逆 行列の列を入れ替える必要が生じる.実際にプログラムでは,以下のようにする.

	/* ---- 列の入れ替え(逆行列) -------------------------- */
	for(j=n ; j>=1 ; j--){
	   if(j != row[j]){
	      for(i=1 ; i<=n ; i++){
	         temp = a[i][j];
	         a[i][j]=a[i][row[j]];
	         a[i][row[j]]=temp;
	      }
	   }
	}


2.4.3 効率化したガウス・ジョルダン法のソースプログラム

これで,ほとんどガウス・ジョルダン法の計算の関数は,完成である.この関数 は,かなり実用に使えるであろう.残っている問題は, くらいと思う.

これ以上,改良するのは大変なので,ほとんど問題なく使える関数のプログラ ムを以下に示す.

	/* ========= ガウスジョルダン法の関数====================== */

	int gauss_jordan(int n, double a[][MAXN+10], double b[]){

	   int ipv, i, j;
	   double inv_pivot, temp;
	   double big;
	   int pivot_row, row[MAXN+10];

	   for(ipv=1 ; ipv <= n ; ipv++){

	      /* ---- 最大値探索 ---------------------------- */
	      big=0.0;
	      for(i=ipv ; i<=n ; i++){
	         if(fabs(a[i][ipv]) > big){
	            big = fabs(a[i][ipv]);
	            pivot_row = i;
	         }
	      }
	      if(big == 0.0) return 1;
	      row[ipv] = pivot_row;

	      /* ---- 行の入れ替え -------------------------- */
	      if(ipv != pivot_row){
	         for(i=1 ; i<=n ; i++){
	            temp = a[ipv][i];
	            a[ipv][i] = a[pivot_row][i];
	            a[pivot_row][i] = temp;
	         }
	         temp = b[ipv];
	         b[ipv] = b[pivot_row];
	         b[pivot_row] = temp;
	      }

	      /* ---- 対角成分=1(ピボット行の処理) ---------- */
	      inv_pivot = 1.0/a[ipv][ipv];
	      a[ipv][ipv]=1.0;
	      for(j=1 ; j <= n ; j++){
	         a[ipv][j] *= inv_pivot;
	      }
	      b[ipv] *= inv_pivot;

	      /* ---- ピボット列=0(ピボット行以外の処理) ---- */
	      for(i=1 ; i<=n ; i++){
	         if(i != ipv){
	            temp = a[i][ipv];
		    a[i][ipv]=0.0;
	            for(j=1 ; j<=n ; j++){
	               a[i][j] -= temp*a[ipv][j];
	            }
	            b[i] -= temp*b[ipv];
	         }
	      }

	   }

	   /* ---- 列の入れ替え(逆行列) -------------------------- */
	   for(j=n ; j>=1 ; j--){
	      if(j != row[j]){
	         for(i=1 ; i<=n ; i++){
	            temp = a[i][j];
	            a[i][j]=a[i][row[j]];
	            a[i][row[j]]=temp;
	         }
	      }
	   }

	   return 0;
	}

ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年11月1日


no counter