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という関数で絶対値 を比較することにより求めます。

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

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

この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(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}$と同じ処理 をすれば、逆行列に変換できます。具体的には、 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(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}$$

となっているはずです。この状態は、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(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
著者: 山本昌志