ガウスのの単純消去法のプログラムです。

質問

ガウスのの単純消去法のプログラムです。
前進消去の第k段階が終わった段階でaij,biが表示されるようにしたいんですがどうすればいいでしょうか↓
よろしくお願いします。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define ERROR -1
#define EPS 1.0e-15

int gausssimp(double *, double *, double *, int);

int main(void)
{
double *a,*b,*x,val;
char s[32];
int i,j,n,result;

printf("Input size n=");
gets(s);
sscanf(s,"%d",&n);

if((a=(double *)calloc(n*n,sizeof(double)) ) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}

if((b=(double *)calloc(n,sizeof(double))) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}

if((x=(double *)calloc(n,sizeof(double))) == NULL){
fprintf(stderr,"Memory allocation error\n");
exit(-1);
}

/* A,b の成分を入力 */
printf("----- A -----\n");
for(i=0; i<n; i++){
for(j=0; j<n; j++){
printf("a(%2d,%2d)=",i+1,j+1);
gets(s);
sscanf(s,"%lf",&val);
a[n*i+j]=val;
}
}

printf("\n----- b -----\n");
for(i=0; i<n; i++){
printf("b(%2d)=",i+1);
gets(s);
sscanf(s,"%lf",&val);
b[i]=val;
}

/* Gaussの単純消去法による求解 */
result = gausssimp(x,a,b,n);

/* 解の表示 */
if(result == ERROR){
printf("ERROR occurs. pivot 0\n");
} else {
printf("\n----- solution -----\n");
for(i=0; i<n; i++){
printf("x(%2d)=%.8e\n",i+1, x[i]);
}
}

free(a);
free(b);
free(x);

return 0;
}

int gausssimp(double *x, double *a, double *b, int n)
{
int i,j,k;
double tmp,p,sum;

/* step 1: 前進消去 */
/**** 追加 ****/
/* 前進消去の各段階を終えるごとに，式がどのように変化しているかわかるように表示する */
/**************/
for(k=0; k<n-1;k++){

if(a[n*k+k] == 0.0) { /* ピボットの値が0．割り算でエラーが起きる．*/
return ERROR;
} else {
/* k+1番目以降の式から x[k] の項を消去 */
for(i=k+1; i<n; i++){
p=a[n*i+k]/a[n*k+k];
for(j=0; j<n; j++){
a[n*i+j]=a[n*i+j]-p*a[n*k+j];
}
b[i]=b[i]-p*b[k];
printf("a[%d %d]=%d b[%d]=%d",i,j,a[n*i+j]-p*a[n*k+j],i,b[i]-p*a[n*k+j]);
↑これではできませんでした。。。
}
/**********************************/
}
printf("k=%d\n",k);

}

/* step 2: 後退代入 */
for(k=n-1; k>=0; k--){
if(fabs(a[n*k+k]) < EPS) return(ERROR);

sum=0.0;
for(j=k+1; j<n; j++) sum+=a[n*k+j]*x[j];
x[k]=(b[k] - sum)/a[k*n+k];
}

return 0;
}

回答No.1

himajin100000 · ベストアンサー

結論だけ言うとこんな感じ

http://ideone.com/pzotE

======================
【割とどうでもいいこと】

gooだとインデントされずに見づらかったし、一行で括弧つけずにif (x == 0) x = 1;みたいに書かれると個人的に判りづらかったので直し、
C99かつ厳しいモードでコンパイルしてみました。
入力用の行列はソース見ながら適当に決めた。(inputのところに書いた)

http://ideone.com/eCfz8

>prog.c:97: error: format ‘%d’ expects type ‘int’, but argument 4 has type ‘double’
>prog.c:97: error: format ‘%d’ expects type ‘int’, but argument 6 has type ‘double’
>prog.c:79: error: unused variable ‘tmp’

%dってint型欲してるっつーに第4引数がdoubleだっ
%dってint型欲してるっつーに第6引数がdoubleだっ
tmpって変数、宣言しているけど、どこにも使ってねーじゃねーか

って言われたのでちまちま修正
http://ideone.com/JdE0w

gets危ないから使うな、って言われたけどとりあえず放置。

答えは正しい、か。ってことは計算過程が正しくて、表示の部分が何かおかしいっと。
そこに相当するのは…91行目から98行目の部分か。

for(i=k+1; i<n; i++){
p=a[n*i+k]/a[n*k+k];
for(j=0; j<n; j++){
a[n*i+j]=a[n*i+j]-p*a[n*k+j];
}
b[i]=b[i]-p*b[k];
printf("a[%d %d]=%lf b[%d]=%lf",i,j,a[n*i+j]-p*a[n*k+j],i,b[i]-p*a[n*k+j]);
}

ん？
何で代入した後、もう一回引く計算してんだ？
あと、ループの後のjの値ってここで使うのか？
この行が実行された時のjの値ってnだよな？

aとbの表示部を分けて、
表示ルーチンをループの中に入れて
かつ、二重に引いてる引き算をなくしてみた。
a[n*i+j]-p*a[n*k+j]→a[n*i+j]

あと、この入力方式使うと見づらいから関数の最初の処理で改行してみた。本来ココに入るべきではないような気もするけど。

http://ideone.com/pzotE

手元で計算。とりあえず期待通りの数値かな。

25 - 24 * 2 / 5 = 15.4

とまぁ、こんな思考を通って直してます。

ところで、これはC99という規格に沿ってやってます。多分昔の規格だと、変数宣言は関数の最初にしか許されてないと思うけど、
C++やC99では他の場所でも許されていて、こんなコードはエラーを吐かれます。

http://ideone.com/iUqn5

forの部分だけで使うっつってんのに7行目で使ってしまってアウトになる、と。今回の変数jでも活かせるかな、とか。
getsが良しとされない理由は面倒くさいので調べてないです。

【話してないこと。興味があれば調べてね】
●ガウスの消去法で計算できないのはa[i,i]が計算途中で0になったら、ではないです。ピボット選択とかのルーチン組んでみようぜ
●配列の配列 int x[3][5]みたいな奴使うと見やすくなるかもね、というお話。