ガウス消去法のプログラムが上手く動作しない

Question

前提・実現したいこと

ガウス消去法のプログラムを書いたのですが、上手くできません。
方程式は絶対値１以下の一様乱数によって作成されます。

発生している問題・エラーメッセージ

計算途中でnanが現れる。
問題を作り直しても上三角化途中である行が全て０になってしまう。

該当のソースコード

//////////方程式を作るプログラム//////////
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(void)
{
	int m,i,j;
	double z;
	FILE *fp;
	char fname[32];
	
	

	printf("サイズ："); scanf("%d",&m);
	
	printf("ファイル名："); scanf("%s",fname);
	fp=fopen(fname,"w");
	fprintf(fp,"%d\n",m);
	
	srand((unsigned)time(NULL));
	
	for (i=0;i<m;i++){
		for (j=0;j<m+1;j++){
			z=2*(rand() / (double)RAND_MAX) - 1;
			fprintf(fp,"%.6lf\t",z);
			
		}
		fprintf(fp,"\n");
	}
			
	fclose(fp);
	fflush(fp);
	return 0;
	
}





///////////////////ガウス消去法によって方程式を解く////////////////////////////

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//一意解を持つか否かは上三角にしてから行列式を計算して求める


int main(void)
{
	FILE *fp,*gp;
	int m,i,j;
	int piv,q;
	double **a,*b,*x;
	double max,temp,r,det;
	
	fp=fopen("af2-1.txt","r");
	fscanf(fp,"%d",&m);		//m:行列の大きさ
	
	printf("%d\n",m);
	
	//メモリ確保
	a=(double**)malloc(m*sizeof(double*));
	for (i=0;i<m;i++) {
		a[i] =(double*)malloc(m*sizeof(double));
	}
	b=(double*)malloc(m*sizeof(double));
	x=(double*)malloc(m*sizeof(double));
	
	
	//テキストファイルから行列の値を読み取り配列に格納
	//係数行列をa、既知ベクトルをb
	for(i=0;i<m;i++){
		for (j=0;j<m+1;j++){
			if (j<m){
				fscanf(fp,"%lf",&a[i][j]);	
				printf("%.2lf\t",a[i][j]);
			}else if (j==m){
				fscanf(fp,"%lf",&b[i]);
				printf("%.2lf\n",b[i]);
			}
		}
	}
	putchar('\n');
	
	//ピボット選択
	//絶対値の大きいものの行をpivに代入（fabsを用いたほうがいい？）
	for (q=0;q<m-1;q++){
		max=pow(a[q][q],2);
		piv=q;
		for (i=q+1;i<m;i++){
			if (max<pow(a[i][q],2)){
				max=pow(a[i][q],2);
				piv=i;
			}
		}	
		

		//第q行と第piv行を交換
		if (piv!=q){
			//係数行列
			for (i=q;i<m;i++){
				temp=a[q][i];
				a[q][j]=a[piv][i];
				a[piv][i]=temp;
			}
			//既知ベクトル
			temp=b[q];
			b[q]=b[piv];
			b[piv]=temp;
		}
		
		
		//上三角化
		//a[q][q]==0の時を場合分け入る？
		for (i=q+1;i<m;i++){
			r=a[i][q]/a[q][q];
			a[i][q]=0.0;
			for (j=q+1;j<m;j++){
				a[i][j]=a[i][j]-r*a[q][j];
			}
			b[i]=b[i]-r*b[q];
		}
	}
	
	//上三角行列の表示
	for (i=0;i<m;i++){
			for (j=0;j<m+1;j++){
				if (j<m){
					printf("%.2lf\t",a[i][j]);
				}else if (j==m){
					printf("%.2lf\n",b[i]);
				}
			}
	}

	//行列式が０でない時一意解を持つ
	det=a[0][0];
	for (i=1;i<m;i++){
		det*=a[i][i];
	}
	if (det==0){
		printf("一意解を持ちません。\n");
		exit(1);
	}
		
		
	//他の方のを参照したので理解していません	
	//後退代入法
	for(i = m - 1; i >= 0; i--){
         for(j = i + 1; j < m; j++){
            b[i] = b[i] - a[i][j] * b[j];
            a[i][j] = 0.0;
         }
         b[i] = b[i] / a[i][i];
         a[i][i] = 1.0;
	}	
	
	//解の表示
	for(i = 0; i < m; i++){
         printf("x[%d] = %12.8lf\n",i+1,b[i]);
    }

	fclose(fp);
	for (i=0;i<m;i++){
		free(a[i]);
	}
	free(a);
	free(b);
	return 0;
}

試したこと

計算途中をprintfして過程を追ってみましたが挙動が全くわかりません

補足情報（FW/ツールのバージョンなど）

ターミナル
ProductName: Mac OS X
ProductVersion: 10.14.6

Answer 1

行の入れ替えが間違えているみたいです。

*** 62,68 ****
              //係数行列
              for (i=q;i<m;i++){
                  temp=a[q][i];
!                 a[q][j]=a[piv][i];
                  a[piv][i]=temp;
              }
              //既知ベクトル
--- 62,68 ----
              //係数行列
              for (i=q;i<m;i++){
                  temp=a[q][i];
!                 a[q][i]=a[piv][i];
                  a[piv][i]=temp;
              }
              //既知ベクトル

Comments

[VanS]

なんと。。。ありがとうございます！！！