gauss_jordan(short).cpp

数值计算的例子和参考书籍

#include <stdio.h>

#define MAX 1000

double a[MAX][MAX], b[MAX];
int n;

inline void Gauss_Jordan(void)
{
	int i, j, k, ik;
	double max, mkk, temp;
	for(k = 1; k <= n; k++)
	{
		ik = -1, max = -2000000000;
		for(i = k; i <= n; i++)//选取主元 
			if(a[i][k] > max) max = a[i][k], ik = i;
		if(ik != k)//如果ik != k，交换a | b的第k行和第ik行 
		{
			for(i = 1; i <= n; i++) {temp = a[k][i]; a[k][i] = a[ik][i]; a[ik][i] = temp;}
			temp = b[k]; b[k] = b[ik]; b[ik] = temp;
		}
		for(i = 1; i <= n; i++)//计算乘数 
			i != k ? a[i][k] = -a[i][k] / a[k][k] : mkk = 1 / a[k][k];
		for(i = 1; i <= n; i++)//消元计算 
		{
			if(i == k) continue;
			for(j = k + 1; j <= n; j++) a[i][j] = a[i][j] + a[i][k] * a[k][j];
			b[i] = b[i] + a[i][k] * b[k];
		}
		for(j = k; j <= n; j++)//计算主行 
			a[k][j] = a[k][j] * mkk;
		b[k] = b[k] * mkk;
	}
}

int main(void)
{
	freopen("in.txt", "r", stdin);
	freopen("out.txt", "w", stdout);
	int i, j;
	while(scanf("%d", &n) != EOF)
	{
		for(i = 1; i <= n; i++)
			for(j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
		for(i = 1; i <= n; i++) scanf("%lf", b + i);
		Gauss_Jordan();
		for(i = 1; i <= n; i++) printf("%lf ", b[i]);
		putchar('\n');
	}
	return 0;
}