解线性方程组

//Iterator.h

#ifndef Iterator_H
#define Iterator_H

#define DEFAULT_DOUBLE_PRECISION	1E-6
#define DEFAULT_MAX_ITERATION_COUNT (10240 * 128)

/************************************************************************/
/* 解线性方程组类，这里AX=b, A为n阶系数矩阵，b为n*1的方程组右值向量
 * 迭代法(比如：雅克比迭代法、高斯迭代法、超松弛迭代法)求解方程组的
 * 解时，要设置解的精度、最大迭代次数。不同的参数设置对最终解的影响
 * 很大，迭代次数过少的话，得到的解可能误差比较大；精度设置过高，很
 * 可能在最大迭代次数内无法达到所期望精度的解；精度设置过高和者迭代
 * 次数过大，计算量也会越大，但是解往往会更精确。
	
	一些特需的情况下：比如方程组是无解的，精度设置较高，此时尽管迭代
	次数很大，迭代法都不会收敛的；方程组是无解的，精度设置过低，在规
	定的次数内可能会获得近似解并且解的误差也是允许的；方程组无解，不
	管精度和最大迭代次数设置为多少，迭代法都不会收敛，无近似解(我们可
	以接受的)。

	直接利用克莱姆法则求解，会有很多局限性的，比如要求系数矩阵A可逆，
	很多时候A尽管可逆，但|A|的值很小，这样A就是一个病态矩阵，结果往往
	误差很大，可能会溢出，当然基本数据类型表示范围有限，各种计算过程中
	都存在潜在的溢出，这和给定的初始数据关系较大；不过我们可以写些大数
	加减乘除的算法来减轻溢出这个问题，使数据表示的范围更大。

	注意：程序部分代码参考了其他的博客，很多错误也没检查。
 */
/************************************************************************/
class CIterator
{
public:
CIterator ();
~CIterator ();

bool SetEnvironment ();						// 设置计算环境 (如系数矩阵等)
bool JacobianCalculate ();					// 雅克比迭代法计算方程解
bool GaussCalculate();						// 高斯迭代法计算方程解
bool RelaxationCalculate (double omiga);	// 超松弛迭代法计算方程解

void SetJacobiPrecision(double d);
void SetGuassPrecision(double d);
void SetSorPrecision(double d);

void SetMaxJacobiIterationCount(int count);
void SetMaxGuassIterationCount(int count);
void SetMaxSorIterationCount(int count);

void Init(double *coefficient, double *B, int order);

int GetJacobiIterationCount();
int GetGuassIterationCount();
int GetSorIterationCount();

bool Klemm();

template <typename T>
double det(T **mat, const int n);			//求n阶行列式的值

double ** GetMatrixA ();					// 获取系数矩阵
int GetOrder ();							// 获取系数矩阵的阶数
double * GetVectorB ();						// 获取方程组右值向量
double * GetSolution ();					// 获取方程解向量

void printResults(const char* methodName);

protected:
	bool ProcessMatrixCoefficient();

private:
double		**matrixA;			// 系数矩阵
int			order;				// 系数矩阵的阶
double		*vectorB;			// 方程组右值向量
double		*solution;			// 解向量
double		slackVariable;

double		jacobiPrecision;
double		guassPrecision;
double		sorPrecision;

int			maxJacobiIterationCount;
int			maxGuassIterationCount;
int			maxSorIterationCount;

int			jacobiIterationCount;
int			guassIterationCount;
int			sorIterationCount;
};

#endif

//Iterator.cpp
#include "Iterator.h"
#include <assert.h>
#include <iostream>
using namespace std;

CIterator :: CIterator ()
{
	matrixA = NULL;
	order = 0;
	vectorB = NULL;
	solution = NULL;

	slackVariable = 1.5;
	jacobiPrecision = DEFAULT_DOUBLE_PRECISION;
	guassPrecision = DEFAULT_DOUBLE_PRECISION;
	sorPrecision = DEFAULT_DOUBLE_PRECISION;

	maxJacobiIterationCount = DEFAULT_MAX_ITERATION_COUNT;
	maxGuassIterationCount = DEFAULT_MAX_ITERATION_COUNT;
	maxSorIterationCount = DEFAULT_MAX_ITERATION_COUNT;

	jacobiIterationCount = 0;
	guassIterationCount = 0;
	sorIterationCount = 0;
}

CIterator :: ~CIterator ()
{
	// 释放内存空间
	if (!vectorB)
	{
		delete [] vectorB;
		vectorB = NULL;
	}

	if (!solution) 
	{
		delete [] solution;
		solution = NULL;
	}

	if (matrixA!=NULL) 
	{
		for (int i = 0; i < order; i++)
		{
			delete[] matrixA[i];
			matrixA[i] = NULL;
		}

		delete[] matrixA;
		matrixA = NULL;
	}
}

void CIterator::SetJacobiPrecision(double d)
{
	jacobiPrecision = d;
}

void CIterator::SetGuassPrecision(double d)
{
	guassPrecision = d;
}

void CIterator::SetSorPrecision(double d)
{
	sorPrecision = d;
}

void CIterator::SetMaxJacobiIterationCount(int count)
{
	maxJacobiIterationCount = count;
}

void CIterator::SetMaxGuassIterationCount(int count)
{
	maxGuassIterationCount = count;
}

void CIterator::SetMaxSorIterationCount(int count)
{
	maxSorIterationCount = count;
}

bool CIterator::ProcessMatrixCoefficient()
{
	for (int i = 0; i < order; i++)
	{
		if (fabs(matrixA[i][i]) < DEFAULT_DOUBLE_PRECISION)
		{
			int j;
			for (j = 0; j < order; j++)
			{
				if (fabs(matrixA[j][i]) > DEFAULT_DOUBLE_PRECISION)
				break;
			}

			if (j >= order) return false;
			
			for (int k = 0; k < order; k++)
				matrixA[i][k] += matrixA[j][k];

			vectorB[i] += vectorB[j];
		}
	}

	return true;
}

template <typename T>
double CIterator::det(T **mat, const int n)

{
	assert(mat != NULL && n > 0);

	if (n == 1) return (double)mat[0][0];
	else if (n == 2) return (double)(mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]);
	else 
	{
		int i, j, k, flag = 1, col;
		double value = 0.0;

		T **tmpMat = (T **)malloc(sizeof(T *) * (n - 1));
		assert(tmpMat != NULL);
		memset(tmpMat, 0, sizeof(T *) * (n - 1));

		for (i = 0; i < n - 1; i++)
		{
			tmpMat[i] = (T *)malloc(sizeof(T) * (n - 1));
			assert(tmpMat[i] != NULL);
		}

		for (i = 0; i < n; i++)
		{
			for (j = 1; j < n; j++)
			{
				col = 0;
				for (k = 0; k < n; k++)
				{
					if (k != i)
					{
						tmpMat[j - 1][col++] = mat[j][k];
					}
				}
			}

			value += mat[0][i] * det(tmpMat, n - 1) * flag;
			flag = -flag;
		}

		for (i = 0; i < n - 1; i++)
		{
			if(tmpMat[i] != NULL) free(tmpMat[i]);
		}
		if(tmpMat != NULL) free(tmpMat);

		return value;
	}
}

// 设置计算环境 (如系数矩阵等)
bool CIterator::SetEnvironment ()
{
	cout << "系数矩阵阶数: ";
	cin >> order;

	matrixA = new double *[order];
	for (int i = 0; i < order; i++)
	{
		matrixA[i] = new double[order];
	}

	cout<<"系数矩阵: ";
	for (int i = 0; i < order; i++)
	{
		for (int j = 0; j<order; j++)
		cin >> matrixA[i][j];
	}
	cout << endl;

	vectorB = new double [order];
	cout << "b 向量: ";
	for (int i=0; i < order; i++)
		cin >> vectorB[i];
	cout << endl;

	solution = new double [order];

	//各种迭代法最大迭代次数和精度等我就用默认的，可以自己写代码设置，
	//接口都定义好了

	cout << "运算环境设置完成." << endl << endl;

	ProcessMatrixCoefficient();

	return true;
}

void CIterator::Init(double *coefficient, double *B, int order)
{
	int i;
	for (i = 0; i < order * order; i++)
		matrixA[i / order][i % order] = coefficient[i];

	for (i = 0; i < order; i++)
		vectorB[i] = B[i];
}

// 雅克比迭代法计算方程解
bool CIterator :: JacobianCalculate ()
{
	cout << "下面使用雅克比迭代法计算该线性方程组" << endl << endl;

	// 当前满足迭代精度控制的解分量数
	int meetPrecisionCount = 0;

	// 辅助向量，存放上一次迭代的解向量。
	double * auxiliary = new double [order];

	// 初始化解向量
	for (int i=0; i<order; i++)
	{
		auxiliary[i] = 0;
		solution[i] = 1;
	}

	// 当前迭代次数
	int itCur = 0;

	// 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环
	while (meetPrecisionCount<order && itCur < maxJacobiIterationCount)
	{
		// 当前满足迭代精度控制的解分量数清零
		meetPrecisionCount = 0;

		// 将当前解向量存入辅助向量
		for (int i=0; i<order; i++)
			auxiliary[i] = solution[i];

		// 计算新的解向量
		for (int i=0; i<order; i++)
		{
			// 暂存当前解分量
			double temp = solution[i];
			// 清零当前解分量
			solution[i] = 0;

			 // 雅克比迭代计算
			for (int j=0; j<i; j++)
				solution[i] += matrixA[i][j]*auxiliary[j];

			for (int j=i+1; j<order; j++)
				solution[i] += matrixA[i][j]*auxiliary[j];

			solution[i] = (vectorB[i]-solution[i]) / matrixA[i][i];

			// 更新当前满足迭代精度控制的解分量数
			if (fabs(temp-solution[i]) < jacobiPrecision)
				meetPrecisionCount++;
		}

		// 当前迭代次数递增
		itCur++;
	}

	delete [] auxiliary;
	// 若在规定的迭代次数内未完成计算任务，则返回错误。
	if (itCur == maxJacobiIterationCount) return false;
	
	return true;
}

// 高斯迭代法计算方程解
bool CIterator :: GaussCalculate ()
{
	// 当前满足迭代精度控制的解分量数
	int meetPrecisionCount = 0;

	// 初始化解向量
	for (int i=0; i<order; i++)
		solution[i] = 0;

	// 当前迭代次数
	int itCur = 0;

	// 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环
	while (meetPrecisionCount < order && itCur < maxGuassIterationCount)
	{
		// 当前满足迭代精度控制的解分量数清零
		meetPrecisionCount = 0;

		// 计算新的解向量
		for (int i=0; i < order; i++)
		{
			// 暂存当前解分量
			double temp = solution[i];
			// 清零当前解分量
			solution[i] = 0;

			// 高斯迭代计算
			for (int j=0; j<i; j++)
				solution[i] += matrixA[i][j]*solution[j];

			for (int j=i+1; j<order; j++)
				solution[i] += matrixA[i][j]*solution[j];

			solution[i] = (vectorB[i]-solution[i])/matrixA[i][i];

			// 更新当前满足迭代精度控制的解分量数
			if (fabs(temp-solution[i]) < guassPrecision)
				meetPrecisionCount++;
		}

		// 当前迭代次数递增
		itCur++;
	}

	// 若在规定的迭代次数内未完成计算任务，则返回错误。
	if (itCur == maxGuassIterationCount) return false;

	return true;
}

// 超松弛迭代法计算方程解
bool CIterator::RelaxationCalculate (double omiga)
{
	cout << "下面使用超松弛迭代法计算该线性方程组" << endl << endl;

	// 当前满足迭代精度控制的解分量数
	int meetPrecisionCount = 0;
	// 当前满足迭代精度控制的解分量数
	for (int i=0; i<order; i++)
		solution[i] = 0;

	// 当前迭代次数
	int itCur = 0;
	// 若当前满足迭代精度控制的解分量数等于解向量维数或者迭代次数达到最大次数则跳出循环
	while (meetPrecisionCount<order && itCur < maxSorIterationCount)
	{
		// 当前满足迭代精度控制的解分量数清零
		meetPrecisionCount = 0;
		// 计算新的解向量
		for (int i=0; i<order; i++)
		{
			// 暂存当前解分量
			double temp = solution[i];

			// 清零当前解分量
			solution[i] = 0;

			// 超松弛迭代计算
			for (int j=0; j<i; j++) 
				solution[i] += matrixA[i][j]*solution[j];

			for (int j=i+1; j<order; j++)
				solution[i] += matrixA[i][j]*solution[j];

			solution[i] = omiga*(vectorB[i]-solution[i])/matrixA[i][i] + (1-omiga)*temp;
			// 更新当前满足迭代精度控制的解分量数

			if (fabs(temp-solution[i]) < sorPrecision)
				meetPrecisionCount++;

			// 当前迭代次数递增
			itCur++;
		}
	}

	// 若在规定的迭代次数内未完成计算任务，则返回错误。
	if (itCur == maxSorIterationCount) return false;

	return true;
}

//克莱姆法则求解线性方程组的解
bool CIterator::Klemm()
{
	int i, j, k;
	double d = det(matrixA, order);
	if (abs(d) < DEFAULT_DOUBLE_PRECISION) return false;

	double **mat = NULL;
	mat = new double *[order];
	assert(mat);

	for (i = 0; i < order; i++)
	{
		mat[i] = new double[order];
		assert(mat[i]);
	}

	for (i = 0; i < order; i++)
	{
		for (j = 0; j < order; j++)
		{
			mat[j][i] = vectorB[j];
		}

		for (j = 0; j < order; j++)
		{
			for (k = 0; k < order; k++)
			{
				if (k != i)
				{
					mat[j][k] = matrixA[j][k];
				}
			}
		}

		solution[i] = det(mat, order) / d;
	}

	for (i = 0; i < order; i++) if(mat[i] != NULL) delete[] mat[i];
	if(mat != NULL) delete[] mat;

	return true;
}

int CIterator::GetJacobiIterationCount()
{
	return jacobiIterationCount;
}

int CIterator::GetGuassIterationCount()
{
	return guassIterationCount;
}

int CIterator::GetSorIterationCount()
{
	return sorIterationCount;
}

// 获取系数矩阵
double ** CIterator :: GetMatrixA ()
{
	return this->matrixA;
}

// 获取系数矩阵的阶
int CIterator :: GetOrder ()
{
	return this->order;
}

// 获取方程组右值向量
double * CIterator :: GetVectorB ()
{
	return this->vectorB;
}

// 获取方程解向量
double * CIterator :: GetSolution ()
{
	return this->solution;
}

// 打印方程解
void CIterator::printResults(const char* methodName)
{
	cout << methodName <<" 计算成功，打印方程解:  " << endl;
	for (int i=0; i < GetOrder(); i++)
		cout << "x" << i << " = " <<GetSolution()[i] << endl;
	cout<<endl;
}

//main.cpp
#include "Iterator.h"
#include <iostream>
using namespace std;

void main(void)
{
	// 定义迭代类对象
	CIterator iterator;

	// 设置计算环境 (如系数矩阵等)
	iterator.SetEnvironment();

	// 雅克比迭代法计算方程解
	if (!iterator.JacobianCalculate()) cout << "雅克比迭代法, 规定次数内未完成迭代计算" << endl;
	else iterator.printResults("雅克比迭代法");

	// 高斯迭代法计算方程解
	if (!iterator.GaussCalculate()) cout << "高斯迭代法, 规定次数内未完成迭代计算" << endl;
	else iterator.printResults("高斯迭代法");

	// 超松弛迭代法计算方程解
	double omiga = 1.5;
	if (!iterator.RelaxationCalculate(omiga)) cout << "超松弛迭代法, 规定次数内未完成迭代计算" << endl;
	else iterator.printResults("超松弛迭代法");

	//克莱姆法则求解线性方程组的解
	if (!iterator.Klemm()) cout << "系数矩阵不可逆，不能直接用克莱姆法则求解线性方程组的解" << endl;
	else iterator.printResults("克莱姆法则");
}