消元法解n元一次方程组 c++实现

昨天在找已知三个三维点求三角形外接圆的过程中找到了一篇文章求圆心方法文章连接，文章中提到了一个方法中的一步需要求多元一次方程组。在网上找到很多实现，但是效果都不好，在某些情况下都不能正确的计算，于是决定自己写一个。

一下是代码，主要就是按照手工消元计算的思路实现。

代码还有些不完善的地方，比如没有对于无解的判断，因为足够自己用了嫌麻烦没有写。

#include<iostream>
#include<math.h>
#include<fstream>
#include<stdlib.h>

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
	double pre = 1e-7;
	int n,m;
	cout<<"输入方程组介数：";
	cin>>n;	
	double a[3][4];

	double r[100];
	int i,j,k,l;
	cout<<"输入增广矩阵："<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n+1;j++){
			cin>>a[i][j];
		}
	}
	cout<<"输入方程组介数：";
	cout<<n<<endl;
	cout<<"输入增广矩阵："<<endl;
	for(i=0;i<n;i++){
		for(j=0;j<n+1;j++){
			cout<<a[i][j]<<"   ";
		}
		cout<<endl;
	}

	int posit[100]; 
	for (i = 0; i < 100; i++)
	{
		posit[i] = -1;
	}
	//i代表元,表示的是列
	for (i = 0; i < n; i++)
	{
		int zeroCount = 0;
		//标记i元第一个非零行
		for (int k = i; k < n; k++)
		{
			if (abs(a[k][i]) > pre)
			{
				posit[i] = k; 
				break;
			}
		}
		if (posit[i] == -1)
		{
			//cout<<"wujie"<<endl;
		}else
		{
			//hang
			for (j = i; j < n; j++)
			{
				if (j != posit[i])
				{
					double scale = a[j][i]/a[posit[i]][i];
					for (k = i; k < n+1; k++)
					{
						a[j][k] -= scale*a[posit[i]][k];
					}
				}
			}
			for (l = 0; l < n+1; l++)
			{
				double t = a[posit[i]][l];
				a[posit[i]][l] = a[i][l];
				a[i][l] = a[posit[i]][l];
			}
		}
		
	}

	for (i = n-1; i >= 0; i--)
	{
		for (j = i+1; j < n; j++)
		{
			a[i][n] -= r[j]*a[i][j];
		}
		//zuihouyibu
		if (abs(a[i][i]) < pre)
		{
			r[i] = 0.0;
		}else
		{
			r[i] = a[i][n]/a[i][i];
		}
	}
	cout<<"结果 "<<endl;
	for (i = 0; i < n; i++)
	{
		cout<<r[i]<<" ";
	}
	return 0;
}