Gauss消去法求解線性方程組

#include<iostream>
using namespace std;

#include<math.h>

void debug_print(float *c,int n);
bool ColPivot(float *,int,float *);

int main()
{
	float *x = new float[3*sizeof(float)];
	float c[3][4] ={10,-2,-1,3,
			-2,10,-1,15,
			-1,-2,5,10};

	ColPivot(c[0],3,x);

	for(int i=0;i<=2;i++ )
		cout<<"x("<<i<<")="<<x[i]<<endl;
	getchar();

	return 0;
}

void debug_print(float *c,int n)
{
	for(int i=0;i<n;i++)
	{
		for(int k=0;k<n+1;k++)
			printf("%f ",*(c+i*(n+1)+k));
		printf("\n");
	}
	printf("\n");
}

/**
 * Gauss消去法是目前求解中小規模線性方程組（即階數不要太高，例如不超過1000）最常用的方法，
 * 一般用于系數矩陣稠密（即矩陣的絕大多數元素都是非零的）而又沒有任何特殊結構的線性方程組。 
 * Gauss消去法實質上計算三角分解的方法。大緻可以分為兩步：消元和回代。
 * 1、消元和回代的算法的實作 
 * 第一步：消元（将系數矩陣經過一系列的初等變換變成上三角矩陣），方法如下： 
 * 第二步：回代（依次解出X3、 X2 、X1） 
 * 2、每次對第K行之後的方程進行消元前要進行兩項操作： 
 * ① 檢查方程組第K到第N行方程中的第K列上元素是否為0，若是，則方程組無解或得不到為唯一解。 
 * ② 為了使方程組的解具有較小的計算誤差，必須使系數矩陣主對角線上的元素的絕對值為最大。
 *    為此，在每消去第K行、第K列上主對角線元素以下的元素時，都要檢查一下主對角線上的元素
 *    是否比其下的同列元素都大，若不是最大，則把最大元素所在行調換到第K行上（最大元素作為主對角線元素）
 */

bool ColPivot( float *c, int n, float *x )
{
	for(int i=0;i<n-1;i++)
	{
		//根據2 ②判斷主對角線元素fabs是否最大
		int col = i;
		for(int j=i+1;j<n;j++)
		{
			if( fabs(*(c+j*(n+1)+i)) > fabs(*(c+col*(n+1)+i)) )
				col = j;
		}
		//如果不是，進行行交換。注：包括常數向量
		if(col!=i)
		{
			for(int k=i;k<n+1;k++)
			{
				float tmp = *(c+col*(n+1)+k);
				*(c+col*(n+1)+k) = *(c+i*(n+1)+k);
				*(c+i*(n+1)+k) = tmp;
			}
		}

		//系數為0，無解或無窮多解。結束！
		if(*(c+i*(n+1)+i) == 0)
			return false;

		//對剩餘行消元
		for(int j=i+1;j<n;j++)
		{
			float p = *(c+j*(n+1)+i) / *(c+i*(n+1)+i);
			for(int k=i;k<n+1;k++)
			{
				*(c+j*(n+1)+k) -= p * (*(c+i*(n+1)+k)); 
			}
		}
		//debug_print(c,n);
	}

	//代入求解方程
	for(int i=n-1;i>=0;i--)
	{
		for(int k=n-1;k>i;k--)
		{
			*(c+i*(n+1)+n) -= (*(c+i*(n+1)+k)) * x[k];
		}
		x[i] = *(c+i*(n+1)+n) / *(c+i*(n+1)+i);
	}

	return true;
}