ガウスジョーダン法を使った逆行列の導出

ガウスジョーダン法を用いた逆行列の導出をするプログラムを教えてほしいです。

c言語でガウスジョーダン法を用いた逆行列を求めるプログラムを作成しました。
８行８列の行列
{0,0,1,0,0,0,0,0}
{1,0,1,0,0,0,2,0}
{1,1,1,0,0,0,-0.7,-0.7}
{0,1,1,0,0,0,0,-0.3}
{0,0,0,0,0,1,0,0}
{0,0,0,1,0,1,0,0}
{0,0,0,1,1,1,-1,-1}
{0,0,0,0,1,1,0,-1}
を引数にして実行したところ、
-1.000000 1.000000 0.000000 0.000000 -2.000000 2.000000 -2.000000 2.000000
-1.750000 0.750000 -0.750000 1.750000 -0.975000 0.975000 -0.975000 0.975000
1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 -1.000000 1.000000 0.000000 0.000000
-2.500000 2.500000 -2.500000 2.500000 -4.250000 3.250000 -3.250000 4.250000
0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 0.000000 0.000000 -1.000000 1.000000 -1.000000 1.000000
-2.500000 2.500000 -2.500000 2.500000 -3.250000 3.250000 -3.250000 3.250000
と出力されたのですが、これだと５列目からの値が手計算したときと異なってしまいます。
いろいろと変数を確認したりしたのですがどうにも治りませんでした。
この場合どこを見たらよいですか？

c
1/*逆行列*/
2double Gauss_Jordan(double mat[N][N],double inv_a[N][N]){
3	double inv_pivot,temp,big;
4	int ipv,i,j,pivot_row,row[N];
5
6	/*単位行列作成*/
7	for(i=0;i<N;i++){
8		for(j=0;j<N;j++){
9			if(i==j){
10				inv_a[i][j]=1.0;
11			}
12			else{
13				inv_a[i][j]=0.0;
14			}
15		}
16	}
17
18	for(ipv=0;ipv<N;ipv++){
19		/*最大値探索*/
20		big=0.0;
21		for(i=ipv;i<N;i++){
22			if(fabs(mat[i][ipv])>big){//abs=絶対値
23				big=fabs(mat[i][ipv]);
24				pivot_row=i;
25			}
26		}
27		if(big==0.0){return(0);}
28		row[ipv]=pivot_row;
29
30		/*行の入れ替え*/
31		if(ipv!=pivot_row){
32			for(i=0;i<N;i++){
33				temp=mat[ipv][i];
34				mat[ipv][i]=mat[pivot_row][i];
35				mat[pivot_row][i]=temp;
36				temp=inv_a[ipv][i];
37				inv_a[ipv][i]=inv_a[pivot_row][i];
38				inv_a[pivot_row][i]=temp;
39			}
40		}
41
42		/*ピボット行の処理(対角成分を1にする)*/
43		inv_pivot=1.0/mat[ipv][ipv];
44		for(j=0;j<N;j++){
45			mat[ipv][j]*=inv_pivot;
46			inv_a[ipv][j]*=inv_pivot;
47		}
48
49		/*ピボット行以外の処理(対角成分が1以外のところを0にする)*/
50		for(i=0;i<N;i++){
51			if(i!=ipv){
52			temp=mat[i][ipv];
53				for(j=0;j<N;j++){
54				mat[i][j]-=temp*mat[ipv][j];
55				inv_a[i][j]-=temp*inv_a[ipv][j];
56				}
57			}
58		}
59	}
60	return(1);
61}

行動規範の内容に同意します

回答1件

ベストアンサー

ことちらでそのGauss_Jordanをコピーして実行した結果は下記の通りです。
結果がちがいますね。

-1.000 +1.000 +0.000 +0.000 +2.000 -2.000 +2.000 -2.000
-1.750 +0.750 -0.750 +1.750 +2.025 -2.025 +2.025 -2.025
+1.000 +0.000 +0.000 +0.000 +0.000 +0.000 +0.000 +0.000
+0.000 +0.000 +0.000 +0.000 -1.000 +1.000 +0.000 +0.000
-2.500 +2.500 -2.500 +2.500 +5.750 -6.750 +6.750 -5.750
+0.000 +0.000 +0.000 +0.000 +1.000 +0.000 +0.000 +0.000
+0.000 +0.000 +0.000 +0.000 -1.000 +1.000 -1.000 +1.000
-2.500 +2.500 -2.500 +2.500 +6.750 -6.750 +6.750 -6.750

C
1void DispMat( double a[N][N])
2{
3	int i, j;
4	for (i = 0; i < N; i++) {
5		for (j = 0; j < N; j++) {
6			printf(" %+.3f", a[i][j]);
7		}
8		printf("\n");
9	}
10}
11
12int main(void) {
13	double a[N][N] = {
14	{ 0,0,1,0,0,0,0,0 }				,
15	{1, 0, 1, 0, 0, 0, 2, 0}		,
16	{1, 1, 1, 0, 0, 0, -0.7, -0.7}	,
17	{0, 1, 1, 0, 0, 0, 0, -0.3}		,
18	{0, 0, 0, 0, 0, 1, 0, 0}		,
19	{0, 0, 0, 1, 0, 1, 0, 0}		,
20	{0, 0, 0, 1, 1, 1, -1, -1}		,
21	{0, 0, 0, 0, 1, 1, 0, -1}		,
22	};
23	double inv[N][N];
24	Gauss_Jordan(a, inv);
25	DispMat(inv);
26}