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にすることになる。最初であろうが、途中であろうが同じようにプログラム は書かなくてはならない。具体定期には、

$$\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;

• $ipv$行の$j$列、 $j=1,2,3,\cdots,n$を処理するために、for文を用い たループになっている。
• a[ipv][j] $\leftarrow$ inv_pivot*a[ipv][j]を、通常、C言語では a[ipv][j] *= inv_pivotと書く。あるいは、a[ipv][j] = inv_pivot*a[ipv][j]と書いても良い。前者の方が少しかっこいい。

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

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

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

$$\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}$倍したもので引けば良いのである。

この処理の実際のプログラムは次のようになる。これで、 $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つのfor文で$i$行$j$列を処理する。
• $ipv$行は対角成分を1にすることで処理が済んでいるので、この行は処理 をしてはならない。そこで、if文を用いて $i \neq ipv$の時、列の 処理をするルーチンが実行されるようになっている。$i=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になっている部分を崩さないためである。

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

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

	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(0);}
	row[ipv] = pivot_row;

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

• bigにその列の絶対値の最大が入る。
• 最大値(新しいピボット)がある行は、pivot_rowである。
• ピボットがゼロの場合、行列は特異となる。その場合、処理を中断して、 呼び出し側にゼロを返す。行列が正則な場合、1を返すことにするが、これはこの 関数の最後に書きく。
• $ipv$番目に最大値になった行を、配列row[ipv]に入れる。 これは、後で使うことにする。

2.2.2 行の交換

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

$$\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 部分ピボット選択付きガウス・ジョルダン法のソースプログ ラム

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

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

	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(0);}
	      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(1);
	}

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

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

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

これを実現するためには、以下の2つのことをすればよいだけである。

• 単位行列を作成する。
• 作成された単位行列を、係数行列 $\boldsymbol{A}$を同じ操作をする。

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}$と同じ処理 をすれば、逆行列に変換でききる。具体的には、 節のプログラムを以下のように書き換える。

• 部分ピボット選択で係数行列の行の交換を行ったならば、同じように inv_aも行の交換を行う。この部分の処理を
  

  	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;           /* -- これ追加 -- */
  

  と書き換えます。
• 対角成分を1にするためにピボット行をピボットの値で割るところでは、 同じ値で割る。これも、
  

  	a[ipv][j] *= inv_pivot;
  	inv_a[ipv][j] *= inv_pivot;          /* -- これ追加 -- */
  

  と書き換える。
• ピボット行以外のピボット列を0にするために、ピボット行の定数倍を 引くところでも同じ操作をする。これも、
  

  	a[i][j] -= temp*a[ipv][j];
  	inv_a[i][j] -= temp*inv_a[ipv][j];   /* -- これ追加 -- */
  

  と書き換える。

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(0);}	
	      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(1);
	}

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{aligned}\boldsymbol{A}^\prime = \begin{pmatrix}1 & 0 & 0... ...dast & \circledast & 0 & 0 & 0 & 1 \\ \end{pmatrix} \end{aligned}$$

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

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

としたくなる。そうすると、メモリーが半分で済む。これは、 $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];
	   }
	}

このようにすると必要なメモリーは、節のプログラ ムに比べて半分になる。さらに、計算時間も半分になる。 節では、計算結果が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}$の逆行列ということは、乗算すると下のように単位行列になる。

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

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

$$i \cdot j = \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(0);}	
	      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(1);
	}

ホームページ: Yamamoto's laboratory
著者: 山本昌志