C++写矩阵求逆

systemvue是一款专用的雷达，通信仿真软件，其底层的封装的模块库是用c++写的。

听博士哥哥说该软件目前处于推广阶段，由于目前模块库匮乏，所以比较蛋疼，就在这蛋疼的时候，我导师接到了这块烫手的山芋，然后更蛋疼的是，把这个写库的任务分给了我。。。。于是乎我又重拾C++大法，被迫暂停大数据比赛（我的穿衣搭配啊），全力转到写库上面来。最近接到博士哥哥的一个任务，让我写一个矩阵求逆的模块，给systemvue用。不管怎么样，写起来再说。

#include <iostream>

.......

然后怎么办。。。

脑子里面搜索着矩阵求逆的算法，想到以前本科老师讲过有个克莱默算法，求矩阵的逆的，先求矩阵的行列式，然后按照公式套一下，这是第一种方法。然后想到了算法导论上面貌似有矩阵运算的讨论，打开一开果然有。第28章都是关于矩阵运算的，有求解线性方程的问题，有矩阵求逆的问题，也有对称正定矩阵和最小二乘逼近 。

上面讲的大致算法是通过将矩阵求逆运算转换成线性方程求解问题，比如求A矩阵的逆，AX=E，E为单位阵，X(i)为X的第i列，e(i)为E的第i列，有A*X(i) = e(i)。X为A的逆矩阵，写到这些相信大家明白了，A的逆矩阵X即为n个方程的解组成的，换句话说这n个解是X矩阵的列向量。至于怎么求解线性方程组，算法导论提供了一种LUP分解的方法。

LUP分解：找三个n*n矩阵L，U和P，满足PA=LU，P为置换矩阵，讲A = P(逆)LU，在求解Ax=b时，两边同乘以P，即得LUx=Pb，L为下三角矩阵，U为上三角矩阵，x可以通过简单的循环过程求出。该算法的关键点在于找到LUP这三个矩阵。具体的矩阵求解过程在这里就不详细说明了，有兴趣的同学可以在算法导论上仔细研读。

下面贴下我的C++代码。由于好久没有写过C++代码了，所以写的并不是很好，数组定义应该用vector容器，不应该定长，而且就算定长也不应该用魔数来写，起码定义const int 常量。这里各位不要学不负责任的我，坏坏的说一句，要是输入模块的矩阵维数超过100的话，我的程序就要出错了。。。。

类声明：

#include <iostream>
#include <cstdio>

class Matrix{
public:
	Matrix();
	~Matrix();
	int get_size() const;
	void inverse();
	void show() const;
	void initialize();
	void lupDecomposition(int (*p)[100], double (*left)[100], double (*up)[100]);

private:
	double data[100][100];
	int rowsize, colsize;
	double* lupSolve(double (*left)[100], double (*up)[100], int *position, double *b);
};

类函数定义：

#include "Matrix.h"

Matrix::Matrix()
{
	rowsize = 3;
	colsize = 3;
	for (int i = 0; i < rowsize; i++) 
		for (int j = 0; j < colsize; j++) {
			data[i][j] = (i + j ) * j + 1;
		}	
}

Matrix::~Matrix()
{
}

void Matrix::initialize()
{
	using namespace std;
	int myrowsize = 0,	mycolsize = 0;
	cout << "input the size of Matrix (row, col) : " << endl;
	cin >> myrowsize >> mycolsize;
	rowsize = myrowsize;
	colsize = mycolsize;
	cout << "input the value of Matrix !" << endl;
	for (int i = 0; i < myrowsize; i++) 
		for (int j = 0; j < mycolsize; j++) {
			cin >> data[i][j];
		}
}

void Matrix::show() const
{
	printf("the rowsize is %d, the colsize is %d\n", rowsize, colsize);
	for (int i = 0; i < rowsize; i++) {
		for (int j = 0; j < colsize; j++)
			printf("%f ", data[i][j]);
		printf("\n");
	}
}

int Matrix::get_size() const
{
	return colsize;
}

void Matrix::lupDecomposition(int (*p)[100], double (*left)[100], double (*up)[100])
{
	int n = get_size();
	int position[100] = {0};

	//初始化position数组
	for (int i = 0; i < n; i++) {
		position[i] = i;
	}

	//定义矩阵a， 并且赋值
	double a[100][100] = {0};
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			a[i][j] = data[i][j];

	for (int k = 0; k < n; k++) {
		double max_row_data = -100000;
		int max_row = k;
		for (int i = k; i < n; i++) {
			if (abs(a[i][k]) > max_row_data) {
				max_row_data = abs(a[i][k]);
				max_row = i;
			}
		}

		if (max_row_data == 0) { 
			printf("this Matrix is singular!");
			return;
		}

		int tmp = 0;
		tmp = position[k];
		position[k] = position[max_row];
		position[max_row] = tmp;

		for (int i = 0; i < n; i++) {
			double tmp = 0;
			tmp = a[k][i];
			a[k][i] = a[max_row][i];
			a[max_row][i] = tmp;
		}

		for (int i = k+1; i < n; i++) {
			a[i][k] /= a[k][k];
			for (int j = k+1; j < n; j++)
				a[i][j] -= a[i][k]*a[k][j];
		}
	}

	//给矩阵P赋值
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++) {
			if (j == position[i])
				p[i][j] = 1;
			else
				p[i][j] = 0;
		}

	//给矩阵left赋值
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++) {
			if (i == j)
				left[i][j] = 1;
			else if (i > j)
				left[i][j] = a[i][j];
			else
				left[i][j] = 0;
		}

	//给矩阵up赋值
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++) {
			if (i <= j)
				up[i][j] = a[i][j];
			else
				up[i][j] = 0;
		}

}


double* Matrix::lupSolve(double (*left)[100], double (*up)[100], int *position, double *b)
{
	int n = get_size();
	double x[100] = {0};
	double y[100] = {0};

	for (int i = 0; i < n; i++) {
		double partSum = 0;
		for (int j = 0; j <= i-1; j++)
			partSum += left[i][j]*y[j];
		y[i] = b[position[i]] - partSum;
	}

	for (int i = n-1; i >= 0; i--) {
		double partSum = 0;
		for (int j = i+1; j <= n-1; j++)
			partSum += up[i][j] * x[j];
		x[i] = (y[i] - partSum) / up[i][i];
	}
	return x;
}

void Matrix::inverse()
{
	double left[100][100] = {{1, 0, 0}, {0.2, 1, 0}, {0.6, 0.5, 1}};
	double up[100][100] = {{5, 6, 3}, {0, 0.8, -0.6}, {0, 0, 2.5}};
	int position[100] = {2, 0, 1};
	int p[100][100] = {0};

	/*double *x;
	x = lupSolve(left, up, position, b);
	for (int i = 0; i < 3; i++)
		printf("%f ", x[i]);*/

	lupDecomposition(p, left, up);

	//得到position，和得到单位矩阵额e[n][n]
	int n = get_size();
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			if (p[i][j] == 1)
				position[i] = j;
	double e[100][100] = {0};
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			if (i == j)
				e[i][j] = 1;
			else
				e[i][j] = 0;

	//构造逆矩阵
	double result[100][100] = {0};
	for (int i = 0; i < n; i++) {
		double *x = NULL;
		x = lupSolve(left, up, position, e[i]);
		for (int j = 0; j < n; j++)
			result[j][i] = x[j];
	}

	//打印result数组
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++)
			printf("%f ", result[i][j]);
		printf("\n");
	}

	return;
}