掃き出し法による行列の解の求め方について

[萩お]

C言語のプログラミングについての質問です。初めての投稿なのでお見苦しい点があると思います。
3×3の行列の解を求めるプログラムを作っているのですが、うまくいかず困っています。
ガウスの掃き出し法を使うのですが、枢軸を行の最大値がある列として掃きだします。
自分で作ってみたプログラムは下のようになりました。

<<ここにソースコードが掲載されていましたが、著作権侵害の申し立てがあったため削除しました>>

また、読み込んだファイルは下のようになります。
58 86 53 1422
64 6 69 1052
65 82 21 1109
いくら考えても自分では手も足も出ず困っています。
お願いします。

[かずま]

そのソースは、全角スペースが入っていたり、
{ と } の対応がとれず、コンパイルエラーになります。

質問するときは、コンパイル環境を明記してください。

Visual Studio で VC++ を使っているのなら、
ブレークポイントやステップ実行でデバッグできるでしょう。

Linux で gcc を使っているのなら、gdb や ddd でデバッグできるでしょう。
Eclipse などの統合開発環境(IDE)が使えるかもしれません。

それ以外なら、printf デバッグでしょうか？

コード:

#include <stdio.h>

#define nmax 3

#define DEBUG
#ifdef DEBUG
void print(double a[nmax][nmax+1])
{
	int i, j;
	for (i = 0; i < nmax; i++, putchar('\n'))
		for (j = 0; j <= nmax; j++)
			printf("%10.3f", a[i][j]);
	puts("---");
}
#endif

void hakidashi(double b[nmax][nmax+1], int jp[nmax])
{
	int i, j, k, l, m;
	double bx[nmax];  /*各行の最大値 */
	for (i = 0; i < nmax; i++) {
		bx[i] = 0.0;
		for (j = 0; j < nmax; j++) {
			if (b[i][j] > bx[i]) {
				jp[i] = j;
				bx[i] = b[i][j];
			} else {
				bx[i] = bx[i];
			}
		}  /*i行の最大値bx[i]を求める */
		printf("bx[%d] = %.3f\n", i, bx[i]);  // DEBUG
		for (k = 0; k < nmax+1; k++) {
			b[i][k] = b[i][k] / bx[i];
		}
		for (l = 0; l < nmax; l++) {
			if (l != i) {
				for (m = 0; m < nmax+1; m++) {
					b[l][m] = b[l][m] - b[l][jp[i]] * b[i][m];
				}
			}
		}
		print(b);  // DEBUG
	}
}

void irekae(double c[nmax][nmax+1], int jq[nmax])
{
	int n, o;
	double temp;
	for (n = 0; n < nmax; n++) {
		if (c[n][n] != 1.0) {
			for (o = 0; o < nmax+1; o++) {
				temp = c[n][o];
				c[n][o] = c[jq[n]][o];
				c[jq[n]][o] = temp;
			}
		}
	}  /*単位行列になるように入れ替える */
}

int main(void)
{
	FILE *fp;
	int i1, j1;
	int J[nmax] = { 0 };  /*各行の最大値がある列 */
	double a[nmax][nmax+1];

	fp = fopen("gyouretsu.txt", "r");
	if (!fp) { puts("fopen failed"); return 1; }  // 念のため
	for (i1 = 0; i1 < nmax; i1++) {
		fscanf(fp, "%lf%lf%lf%lf", &a[i1][0], &a[i1][1], &a[i1][2], &a[i1][3]);
	}  /*行列を読み込む */
	hakidashi(a, J);
	irekae(a, J);
	for (i1 = 0; i1 < nmax; i1++) {
		printf("%lf %lf %lf %lf\n", a[i1][0], a[i1][1], a[i1][2], a[i1][3]);
	}  /*結果を出力 */
	fclose(fp);
	return 0;
}

分からないところを具体的に質問してもらうと、
詳しい説明ができると思います。

[萩お]

すいません。コンパイル環境はVS 2017用 x64 Native Tools コマンドプロンプトを使っております。

[萩お]

分からないところは、実行結果が計算されていないことです。コンパイルしてもエラーは出ないのですが、実行するとinfやnanが出てきてしまいます。

[萩お]

質問についての補足
根本から間違っているかもしれませんが、自分でコンパイルを区切って確認していったところ、おそらくi行以外の掃き出しの部分でちゃんと計算されていないのだと思います。

[萩お]

>かずま様
infやnanが出なくなり、計算結果が出るようになりました！
ありがとうございます。
ただ、お作りしていただいたプログラムを実行してみたのですが、実行結果が
bx[0] = 58.000
0.052 0.017 1.000 1.483
49.690 1420.897 0.000 6.000
65.638 1050.879 0.000 82.000
---
bx[1] = 1420.897
0.051 0.000 1.000 1.483
0.035 1.000 0.000 0.004
28.888 0.000 0.000 82.000
---
bx[2] = 28.888
0.000 0.000 1.000 1.483
0.000 1.000 0.000 0.004
1.000 0.000 0.000 2.839
---
1.000000 0.000000 0.000000 2.838548
0.000000 1.000000 0.000000 0.004223
0.000000 0.000000 1.000000 1.482759

となり、各行の最大値が合わないです...

[かずま]

infやnanが出なくなり、計算結果が出るようになりました！

どのようにして inf や nan が出ないようにしたのかが
こちらからは全く分かりません。それを見せてもらえないと、
正しい計算結果が出るようにするための助言もできません。

ただ、お作りしていただいたプログラムを実行してみたのですが、実行結果が

私は、掃き出し法のプログラムを作ったのではありません。
printf デバッグの方法を示しただけです。

プログラムの実行順序が分かるのなら、変数がどのように変化するのかを
追っていき、思っっているものと異なる値になったところに問題があることが
分かるはずです。それがデバッグです。

掃き出し法の手順(アルゴリズム)そのものが分からないのでしょうか？

掃き出し法を愚直にコーディングすると次のようになると思います。

コード:

#include <stdio.h>   // fopen, fclose, fscanf, printf, puts, putchar
#include <stdlib.h>  // exit

#define N 3

void print(double a[N][N+1])
{
	for (int i = 0; i < N; i++, putchar('\n'))
		for (int j = 0; j <= N; j++) printf(" %9.3f", a[i][j]);
	putchar('\n');
}

void Gaussian(double a[N][N+1])
{
	double p;
	int i, j, k;
	for (i = 0; i < N; i++) {
		p = a[i][i];
		if (p == 0) {
			for (j = i + 1; j < N && a[j][i] == 0; j++) ;
			if (j == N) puts("can't solve"), exit(1);
			for (k = 0; k <= N; k++) a[i][k] += a[j][k];
			p = a[i][i];
		}
		for (k = 0; k <= N; k++) a[i][k] /= p;
		for (j = i + 1; j < N; j++)
			for (p = a[j][i], k = 0; k <= N; k++) a[j][k] -= a[i][k] * p;
		print(a);  // DEBUG
	}
	for (i = N - 1; i >= 0; i--) {
		for (j = i - 1; j >= 0; j--)
			for (p = a[j][i], k = 0; k <= N; k++) a[j][k] -= a[i][k] * p;
		print(a);  // DEBUG
	}
}

int main(void)
{
	double a[N][N+1];

	FILE *fp = fopen("gyouretsu.txt", "r");
	if (!fp) return puts("fopen failed"), 1;
	for (int i = 0; i < N; i++)
		for (int j = 0; j <= N; j++) fscanf(fp, "%lf", &a[i][j]);
	fclose(fp);
	print(a);
	Gaussian(a);
	for (int i = 0; i < N; i++) printf(" %9.3f", a[i][N]);
	putchar('\n');
}

pivot が 0 の場合の対処の仕方にはいろいろあると思いますが、
ここでは、0 でない行を見つけて、今の行に足すだけにしています。
無駄な計算もしているので改善の余地がありますが、一応答えは出るようです。

[萩お]

ご返信ありがとうございます。デバックの方法を教えていただいたのに頓珍漢な返事をしてしまい申し訳ありませんでした。デバックの方法を示していただいたプログラムではきちんと計算されるので、もしかするとスペースが開いていないなどの問題で計算が実行されなかったのかもしれません。あと掃き出し法についてなのですが、各行の絶対値の最大値がある軸を枢軸とする条件で、三角行列にしてからの解法は使いたくても使えません...　
頂いたアドバイスを参考にもう一度じっくり考えてみようと思います。

[萩お]

よく考えなおしてミスを改善することができ、自己解決できました。
ありがとうございました。

[かずま]

フォーラム(掲示板)ルールに従ってください。

d. 義務行為

　"C言語何でも質問掲示板"でのみ適用される事項
　・トピックを立て、解決した場合は「解決しました」とだけ書かず、
　　どうやって解決したか他の人に分かるように書いてからトピックを終了して下さい。

後で使用しない変数を計算しないように修正してみました。

コード:

#include <stdio.h> // fopen, fclose, fscanf, printf, puts, putchar
#include <math.h>  // fabs

#define N 3

void print(double a[N][N+1])
{
	for (int i = 0; i < N; i++, putchar('\n'))
		for (int j = 0; j <= N; j++) printf(" %9.3f", a[i][j]);
	putchar('\n');
}

int Gaussian(double a[N][N+1])
{
	double p;
	int i, j, k;
	for (i = 0; i < N; i++) {
		for (k = j = i; ++j < N; )
			if (fabs(a[j][i]) > fabs(a[k][i])) k = j;
		if (k != i)
			for (j = 0; j <= N; j++)
				p = a[i][j], a[i][j] = a[k][j], a[k][j] = p;
		else if (a[i][i] == 0) return 1;
		for (j = i; ++j < N; )
			for (p = a[j][i] / a[i][i], k = i; ++k <= N;)
				a[j][k] -= a[i][k] * p;
		//print(a);  // DEBUG
	}
	for (i = N; --i >= 0; ) {
		for (a[i][N] /= a[i][i], j = i; --j >= 0; )
			a[j][N] -= a[j][i] * a[i][N];
		//print(a);  // DEBUG
	}
	return 0;
}

int main(void)
{
	double a[N][N+1];

	FILE *fp = fopen("gyouretsu.txt", "r");
	if (!fp) return puts("fopen failed"), 1;
	for (int i = 0; i < N; i++)
		for (int j = 0; j <= N; j++) fscanf(fp, "%lf", &a[i][j]);
	fclose(fp);
	print(a);
	if (Gaussian(a)) return puts("can't solve"), 2;
	for (int i = 0; i < N; i++) printf(" %9.3f", a[i][N]);
	putchar('\n');
	return 0;
}