gauss_jordan2.cpp

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

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

#define MAX 1000

#define INF 1e80
#define eps 1e-8

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

bool Gauss_Jordan(double a[][MAX], double *b, int n)
{
	for(int k = 1; k <= n; k++)
	{
		double max = -INF, mkk;
		int ik = -1;
		for(int i = k; i <= n; i++)//选取主元 
		{
			if(a[i][k] > max)
			{
				max = a[i][k];
				ik = i;
			}
		}
		if(fabs(max) < eps) return false;
		if(ik != k)//如果ik != k，交换a | b的第k行和第ik行 
		{
			for(int i = 1; i <= n; i++)
			{
				double temp = a[k][i];
				a[k][i] = a[ik][i];
				a[ik][i] = temp;
			}
			double temp = b[k];
			b[k] = b[ik];
			b[ik] = temp;
		}
		for(int i = 1; i <= n; i++)//计算乘数 
		{
			if(i != k)
			{
				a[i][k] = -a[i][k] / a[k][k];
			}
			else
			{
				mkk = 1 / a[k][k];
			}
		}
		for(int i = 1; i <= n; i++)//消元计算 
		{
			if(i != k)
			{
				for(int 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(int j = k; j <= n; j++)//计算主行 
		{
			a[k][j] = a[k][j] * mkk;
		}
		b[k] = b[k] * mkk;
	}
	return true;
}

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(a, b, n);
		for(i = 1; i <= n; i++)
		{
			printf("%lf ", b[i]);
		}
		putchar('\n');
	}
	return 0;
}