ガウスの消去法で連立方程式を解くプログラムに途中式を書かせたい

Question

ガウスの消去法を使ってn元連立方程式を解くプログラムです。

計算過程を途中式としてoutput1.datに出力したいのですが、うまくいきません。
アドバイスをお願いしたいです。

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

void input_matrix(int n,double **a,char c,FILE *fin, FILE *fout);
void input_vector(int n,double *b,char c,FILE *fin, FILE *fout);
double **dmatrix(int nr1,int nr2,int nl1,int nl2);
void free_dmatrix(double **a,int nr1,int nr2,int nl1,int nl2);
double *dvector(int i,int j);
void free_dvector(double *a,int i);
double *gauss(int n,double **a,double *b);

int main(void){
	FILE *fin,*fout;
	double **a,*b;
	int i,n;
	
	if((fin=fopen("input1.dat","r"))==NULL){
		printf("ファイルが見つかりません:input1.dat\n");
		exit(1);
	}
	if((fout=fopen("output1.dat","w"))==NULL){
		printf("ファイルが作成できません:output1.dat\n");
		exit(1);
	}
	
	fscanf(fin,"%d",&n);
	a=dmatrix(1,n,1,n);
	input_matrix(n,a,'A',fin,fout);
	fscanf(fin,"%d",&n);
	b=dvector(1,n);
	input_vector(n,b,'b',fin,fout);
	b=gauss(n,a,b);
	
	fprintf(fout,"Ax=bの解は次の通りです\n");
	for(i=1;i<=n;i++){
		fprintf(fout,"%f\n",b[i]);
	}
	
	fclose(fin);
	fclose(fout);
	
	free_dmatrix(a,i,n,i,n);free_dvector(b,1);
	return 0;
}


void input_matrix(int n,double **a,char c, FILE *fin,FILE *fout){
	int i,j;
	
	fprintf(fout,"行列%cは次の通りです\n",c);
	for(i=1;i<=n;i++){
		for(j=1;j<=n;j++){
			fscanf(fin,"%lf",&a[i][j]);
			fprintf(fout,"%5.2f\t",a[i][j]);
		}
		fprintf(fout,"\n");
	}
}



void input_vector(int n,double *b,char c, FILE *fin,FILE *fout){
	int i;
	
	fprintf(fout,"ベクトル%cは次の通りです\n",c);
	for(i=1;i<=n;i++){
		fscanf(fin,"%lf",&b[i]);
		fprintf(fout,"%5.2f\t",b[i]);
		fprintf(fout,"\n");
	}
}



double **dmatrix(int nr1, int nr2,int nl1, int nl2){
  int i,nrow,ncol;
  double **a;

  nrow=nr2-nr1+1;
  ncol=nl2-nl1+1;
  if((a=malloc(nrow*sizeof(double*)))==NULL){
    printf("メモリが確保できません。");
    exit(1);
  }
  a=a-nr1;

  for(i=nr1;i<=nr2;i++){
    a[i]=malloc(ncol*sizeof(double));
   }
  for(i=nr1;i<=nr2;i++){
    a[i]=a[i]-nl1;
   }
  
  return a;
}


void free_dmatrix(double**a, int nr1, int nr2,int nl1, int nl2){
  int i;

  for(i=nr1;i<=nr2;i++){
    free((void*)(a[i]+nl1));
  }
  free((void*)(a+nr1));
}


double *dvector(int i,int j){
  double *a;
  if((a=malloc(((j-i+1)*sizeof(double))))==NULL){
    printf("メモリが確保できません(from dvectoer)\n");
    exit(1);
  }
  return(a-i);
}


void free_dvector(double *a,int i){
  free((void*)(a+i));
}


double *gauss(int n,double **a,double *b){
	int i,j,k,ip;
	double alpha,tmp;
	double amax,eps=pow(2.0,-50.0);
	
	for(k=1;k<=n-1;k++){
		amax=fabs(a[k][k]);ip=k;
		for(i=k+1;i<=n;i++){
			if(fabs(a[i][k])>amax){
				amax=fabs(a[i][k]);ip=i;
			}
		}
		if(amax<eps){
			printf ("入力した行列は正則ではない!\n");
		}
		
		if(ip!=k){
			for(j=k;j<=n;j++){
				tmp=a[k][j];a[k][j]=a[ip][j];a[ip][j]=tmp;
			}
			tmp=b[k];b[k]=b[ip];b[ip]=tmp;
		}
		
		for(i=k+1;i<=n;i++){
			alpha=-a[i][k]/a[k][k];
			for(j=k+1;j<=n;j++){
				a[i][j]=a[i][j]+alpha*a[k][j];
			}
		b[i]=b[i]+alpha*b[k];
		}
	}
	
	b[n]=b[n]/a[n][n];
	for(k=n-1;k>=1;k--){
		tmp=b[k];
		for(j=k+1;j<=n;j++){
			tmp=tmp-a[k][j]*b[j];
		}
		b[k]=tmp/a[k][k];
	}
	return b;
}

Comments

[episteme]

「うまくいかない」を説明してください。
計算途中の値を出力しているのはコードのどこですか?

[Rin001ss]

うまくいかないではなく「どこにfprintfすればいいかがわからない」です、すみません

[episteme]

出力先は FILE* fout じゃありませんか?

それとも「どこにfprintfすれば...」じゃなくて「どこでfprintfすれば...」ですか?

[Rin001ss]

その通りです

Answer 1

void print_matrix(int n, double** a, FILE* fout) {
  // n*n 行列 a を fout に出力する
}

void print_vector(int n, double* b, FILE* fout) {
  // 要素数 n のベクトル b を fout に出力する
}

を用意し、gauss内のテキトーなところで呼べばいい。
"テキトーなところ"とは ここで途中結果が欲しいナー と思ったところ。