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

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

この処理の実際のプログラムは次のようになる．これで， $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になっている部分を崩さないためである．

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

最大値は，$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;

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

• bigにその列の絶対値の最大が入る．
• 最大値(新しいピボット)がある行は，pivot_rowである．
• ピボットがゼロの場合，行列は特異となる．その場合，処理を中断して，呼び出 し側に1を返す．行列が正則な場合，0を返すことにするが，これはこの関数の 最後に書く．
• $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 部分ピボット選択付きガウス・ジョルダン法のソースプログ ラム

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

• 部分ピボット選択で係数行列の行の交換を行ったならば，同じように 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 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{aligned}\boldsymbol{A}^\prime = \begin{pmatrix}1 & 0 & 0... ...dast & \circledast & 0 & 0 & 0 & 1 \\ \end{pmatrix} \end{aligned}$$

となっているはずである．この状態は，2.3.3節の 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];
	   }
	}

このようにすると必要なメモリーは，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}$の逆行列ということは，乗算すると下のように単位行列になる．

$$\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 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
著者: 山本昌志