计算方法之用高斯列主元消去法求线性方程组

/*************************************
* 用高斯列主元消去法求线性方程组
* 
* 2*x1 + 2*x2 + 3*x3 = 3
*{4*x1 + 7*x2 + 7*x3 = 1
* -2*x1+ 4*x2 + 5*x3 = -7
*
**************************************/

#include<stdio.h>
#include<math.h>
#include<conio.h>
#include<string.h>

#define N 3

int main() {
	static double a[N][N] = { { 2, 2, 3 }, { 4, 7, 7 }, { -2, 4, 5 } };
	double b[N] = { 3, 1, -7 };
	double x[N] = { 0, 0, 0 };
	double r, s, e;
	int k, i, j, p, flag = 1;

	for (k = 0; k < N - 1; k++) {
		p = k;
		e = a[k][k];
		for (i = k; i < N; i++)
			if (fabs(a[i][k]) > e) {
				e = fabs(a[i][k]);
				p = i;
			}
		for (j = k; j < N; j++) {
			s = a[k][j];
			a[k][j] = a[p][j];
			a[p][j] = s;
		}
		s = b[k];
		b[k] = b[p];
		b[p] = s;
		if (a[k][k] == 0) {
			printf("Gauss-Method does not run!");
			flag = 0;
			break;
		} else {
			for (i = k + 1; i < N; i++) {
				r = a[i][k] / a[k][k];
				if (a[k][k] != 0) {
					for (j = k; j < N; j++)
						a[i][j] = a[i][j] - r * a[k][j];
				}
				b[i] -= r * b[k];
			}
		}
	}

	if (flag) {
		x[N - 1] = b[N - 1] / a[N - 1][N - 1];
		for (i = N - 2; i >= 0; i--) {
			s = b[i];
			for (j = i + 1; j < N; j++)
				s -= a[i][j] * x[j];
			x[i] = s / a[i][i];
		}
		for (i = 0; i < N; i++)
			printf("x[%d] = %10.7f\n", i + 1, x[i]);
	}
	return 0;
}