[C++知识库]c++ 矩阵运算 Matrix(矩阵)类的粗略实现

矩阵运算在计算机图形和视觉中用着十分广泛的应用，成熟的第三方库有很多比如：Eigen、opencv、ViennaCL等。

本矩阵运算类大致实现了以下基本内容：
矩阵的初始化构造
矩阵的加法、减法、乘法；
矩阵的转置、矩阵的逆、通过增广矩阵高斯消元求解方程组的解、矩阵行列式的计算等简单线代应用；
虽较为粗糙但拿来温习线代知识还是不错的。
头文件

#pragma once
#include <cmath>
#include <assert.h>//如果方程无解或矩阵没有逆或运算不符合数学逻辑 报错
using namespace std;
#define ERROR 1e-12
#define ESP(val) (fabs(val) <= ERROR) ? 0.0 : val //过于小的数字视为0
class Mat {
private:
	double **p;//二维指针
	int rows, cols;//矩阵的尺寸
	void initMatrix();//矩阵的初始化
public:
	
	Mat(int r, int c);
	Mat(int r, int c, double val);
	Mat(const Mat & m);
	virtual ~Mat();
public:
	Mat operator=(const Mat& arr);
	Mat operator=(double *arr);
	Mat operator*(const Mat& arr)const;
	Mat operator*=(const Mat& arr);
	Mat operator+=(const Mat& arr);
	Mat operator-=(const Mat& arr);
	Mat operator+(const Mat& arr)const;
	Mat operator-(const Mat& arr)const;
public:
	void Solve(Mat& outArr)const;				//求解n元一次方程
	Mat T()const;							//转置矩阵
	Mat inv()const;								//矩阵的逆
	double det()const;							//行列式
	int row()const;
	int col()const;
	void swapLine(int r1, int r2);				//行交换
	void Show()const;

};
void Mat::initMatrix() {//分配空间
	p = new double * [rows];
	for (int i = 0; i < rows; i++) {
		p[i] = new double[cols];
	}
}
Mat::Mat(int r = 1, int c = 1) {//构造默认全为0的矩阵
	rows = r, cols = c;
	initMatrix();
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = 0.0;
		}
	}
}
Mat::Mat(int r, int c, double val) {//构造带初始值的矩阵
	rows = r, cols = c;
	initMatrix();
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = val;
		}
	}
}
Mat::Mat(const Mat & m) {
	rows = m.rows;
	cols = m.cols;
	initMatrix();
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = m.p[i][j];
		}
	}
}
Mat::~Mat() {//析构函数
	for (int i = 0; i < rows; ++i) {
		delete[] p[i];
	}
	delete[] p;
}
Mat Mat::operator = (const Mat &arr) {
	if (this == &arr) return *this;//地址相同是一个矩阵

	if (!(rows == arr.rows&&cols == arr.cols)) {	//尺寸不同重新构造
		this->~Mat();								//释放掉先前的矩阵
		rows = arr.rows;
		cols = arr.cols;
		initMatrix();
	}
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = arr.p[i][j];
		}
	}
	return *this;
}
Mat Mat::operator = (double *arr) {
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = *(arr + i*cols + j);
		}
	}
	return *this;
}
Mat Mat::operator *= (const Mat& arr) {
	if (rows != arr.cols) assert(0);
	Mat temp(rows, arr.cols);
	for (int i = 0; i < temp.rows; i++) {
		for (int j = 0; j < temp.cols; j++) {
			for (int k = 0; k < cols; k++) {
				temp.p[i][j] += p[i][k] * arr.p[k][j];
			}
			temp.p[i][j] = ESP(temp.p[i][j]);
		}
	}
	*this = temp;
	return *this;
}
Mat Mat::operator * (const Mat& arr)const {
	if (rows != arr.cols) assert(0);
	Mat temp(rows, arr.cols, 0.0);
	for (int i = 0; i < temp.rows; i++) {
		for (int j = 0; j < temp.cols; j++) {
			for (int k = 0; k < cols; k++) {
				temp.p[i][j] += p[i][k] * arr.p[k][j];
			}
			temp.p[i][j] = ESP(temp.p[i][j]);//忽略精度损失
		}
	}
	return temp;
}
Mat Mat::operator += (const Mat& arr) {
	if (!(rows == arr.rows&&cols == arr.cols)) assert(0);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] += arr.p[i][j];
		}
	}
	return *this;
}
Mat Mat::operator -= (const Mat& arr) {
	if (!(rows == arr.rows&&cols == arr.cols)) assert(0);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] -= arr.p[i][j];
		}
	}
	return *this;
}
Mat Mat::operator + (const Mat& arr)const {
	if (!(rows == arr.rows&&cols == arr.cols)) assert(0);
	Mat temp(rows, cols);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			temp.p[i][j] = p[i][j] + arr.p[i][j];
		}
	}
	return temp;
}
Mat Mat::operator - (const Mat& arr)const {
	if (!(rows == arr.rows&&cols == arr.cols)) assert(0);
	Mat temp(rows, cols);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			temp.p[i][j] = p[i][j] - arr.p[i][j];
		}
	}
	return temp;
}
void Mat::Solve(Mat& outArr)const {					//通过增广矩阵 结合 高斯消元求Ax = b的解
	if (rows != cols - 1) assert(0);				//条件不足无法求解
	Mat A(*this);									//将矩阵拷贝
	/*进行行变换转换为上三角*/
	for (int i = 0; i < A.rows; i++) {
		for (int j = i + 1; j < A.rows; j++) {
			double k;
			if (A.p[i][i] == 0) {			//如果这一个元素是0，尝试与下面的非零行交换
				for (int m = i + 1; m < A.rows; m++) {
					if (A.p[m][i] != 0) {
						A.swapLine(i, m);
						break;
					}
				}
				if (A.p[i][i] == 0) k = 0;
				else k = A.p[j][i] / A.p[i][i];
			}
			else k = A.p[j][i] / A.p[i][i];
			for (int m = i; m < A.cols; m++) {
				A.p[j][m] -= A.p[i][m] * k;
			}
		}
	}
	/*得到 “三角矩阵”*/
	/*自下而上求解*/
	outArr = Mat(A.rows, 1);
	assert(A.p[A.rows - 1][A.rows - 1]);//无解
	for (int i = A.rows - 1; i >= 0; i--) {
		double s = 0.0;
		for (int j = i + 1; j < A.cols - 1; j++) s += outArr.p[j][0] * A.p[i][j];
		if (A.p[i][i] != 0) {
			outArr.p[i][0] = (A.p[i][cols - 1] - s) / A.p[i][i];
		}
	}
}
Mat Mat::T()const {//转置矩阵
	Mat T(cols, rows);
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			T.p[j][i] = p[i][j];
		}
	}
	return T;
}
Mat Mat::inv()const {//矩阵的逆
	if (rows != cols) assert(0);
	Mat E(rows, cols, 0); //生成一个单位矩阵
	for (int i = 0; i < rows; i++) E.p[i][i] = 1.0;
	Mat A(*this);		  //将求逆的矩阵拷贝
	/*进行行变换转换为上三角*/
	for (int i = 0; i < A.rows; i++) {
		for (int j = i + 1; j < A.rows; j++) {
			double k;
			if (A.p[i][i] == 0) {//如果这一个元素是0，尝试与下面的非零行交换
				for (int m = i + 1; m < A.rows; m++) {
					if (A.p[m][i] != 0) {
						A.swapLine(i, m);
						E.swapLine(i, m);
						break;
					}
				}
				if (A.p[i][i] == 0) k = 0;
				else k = A.p[j][i] / A.p[i][i];
			}
			else k = A.p[j][i] / A.p[i][i];
			for (int m = 0; m < A.cols; m++) {

				A.p[j][m] -= A.p[i][m] * k;
				E.p[j][m] -= E.p[i][m] * k;
			}
			
		}
	}
	/*第二次转换为下三角*/
	for (int i = A.rows - 1; i >= 0; i--) {
		for (int j = i - 1; j >= 0; j--) {
			double k;
			if (A.p[i][i] == 0) k = 0;
			else k = A.p[j][i] / A.p[i][i];

			for (int m = A.rows - 1; m >= 0; m--) {
				A.p[j][m] -= A.p[i][m] * k;
				E.p[j][m] -= E.p[i][m] * k;
			}
		}
	}
	/*对角变单位阵*/
	for (int i = 0; i < A.rows; i++) {
		assert(A.p[i][i]);//无逆矩阵
		double k = A.p[i][i];
		A.p[i][i] = 1;
		for (int j = 0; j < A.cols; j++)
			E.p[i][j] /= k;
	}
	return E;
}
double Mat::det()const {
	double res = 1.0;
	Mat A(*this);
	for (int i = 0; i < A.rows; i++) {
		for (int j = i + 1; j < A.rows; j++) {
			double k;
			if (A.p[i][i] == 0) {//如果这一个元素是0，尝试与下面的非零行交换
				for (int m = i + 1; m < A.rows; m++) {
					if (A.p[m][i] != 0) {
						A.swapLine(i, m);
						res *= -1.0;
						break;
					}
				}
				if (A.p[i][i] == 0) k = 0;
				else k = A.p[j][i] / A.p[i][i];
			}
			else k = A.p[j][i] / A.p[i][i];
			for (int m = 0; m < A.cols; m++) {
				A.p[j][m] -= A.p[i][m] * k;
			}

		}
	}
	for (int i = 0; i < A.rows; i++) res *= A.p[i][i];
	return res;
}
int Mat::row()const {
	return rows;
}
int Mat::col()const {
	return cols;
}
void Mat::swapLine(int r1, int r2) {
	for (int i = 0; i < cols; i++) {
		double t = p[r1][i];
		p[r1][i] = p[r2][i];
		p[r2][i] = t;
	}
}
void Mat::Show() const {
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++)
			cout << p[i][j] << " ";
		cout << '\n';
	}
	cout << '\n';
}

test：

#include<iostream>
#include<stdlib.h>
#include<time.h>
#include"Matrix.h"
using namespace std;
int main() {
	srand((unsigned int)time(NULL));
	Mat a(10, 10), b = Mat(3, 4);
	double arr[100], arr1[12] = { 2,1,1,1,6,2,1,-1,-2,2,1,7 };
	for (int i = 0; i < 100; i++) {//随机生成一个 10*10 的矩阵
		arr[i] = rand() % 5;
	}
	a = arr;
	b = arr1;
	cout << "a矩阵：\n";
	a.Show();
	cout << "a矩阵行列式的值: " << a.det() << '\n';
	cout << "a矩阵的逆: \n";
	a.inv().Show();
	cout << "检验A*A-1 = E：\n";
	(a.inv()*a).Show();
	cout << '\n';
	cout << "b矩阵的解(b是由方程组构成的矩阵)：\n";
	Mat x;
	b.Solve(x);
	x.Show();
	cout << '\n';
	cout << "b的转置：\n";
	b.T().Show();
	cout << '\n';
	system("pause");
	return 0;
}