c++ 矩阵运算 Matrix(矩阵)类的粗略实现

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

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

#pragma once
#include<iostream>
#include <cmath>
#include <assert.h>//如果方程无解或矩阵没有逆或运算不符合数学逻辑 报错
using namespace std;

#define ERROR 1e-10
#define ESP(val) (fabs(val) <= ERROR) ? 0.0 : val	//过于小的数字视为0
class Mat {
private:
	double **p;										//二维指针
	int rows, cols;									//矩阵的尺寸
	void initMatrix();								//为矩阵分配空间

	/*************************函数声明****************************/
public:
	/*初始化构造*/
	explicit Mat(int r, int c);						//默认
	explicit Mat(int r, int c, double val);			//浮点
	explicit Mat(int r, int c, int val);			//整型
	Mat(const Mat & m);								//矩阵
	Mat(int r, int c, double *arr);					//数组
	virtual ~Mat();
public:
	/*重载操作运算符*/
	Mat operator=(const Mat& arr);								//矩阵赋值矩阵
	Mat operator=(double *arr);									//数组赋值矩阵(最好提前知道矩阵尺寸)
	bool 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;
	Mat operator-(const Mat& arr)const;
	double* operator[](int index)const;							//[] 访问
	friend istream &operator >> (istream &in, Mat &arr);		//cin输入
	friend ostream &operator << (ostream &out,const Mat &arr);	//cout输出
public:
	void Solve(Mat& outArr)const;				//求解n元一次方程
	Mat T()const;								//转置矩阵
	Mat Inv()const;								//矩阵的逆
	double det()const;							//行列式
	Mat matPow(int n)const;						//矩阵的n次幂
	void Normalizing(int r, int c);				//单位化矩阵
	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(int r, int c, int 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(int r, int c, double *arr) {
	rows = r, cols = c;
	initMatrix();
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			p[i][j] = *(arr + i*cols + 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;
}
bool Mat::operator == (const Mat& arr)const {
	if (!(rows == arr.rows&&cols == arr.cols)) return false;
	for (int i = 0; i < rows; i++) {
		for (int j = 0; j < cols; j++) {
			if (p[i][j] != arr.p[i][j]) return false;
		}
	}
	return true;
}
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;
}
double* Mat::operator[](int index)const {
	try{
		if (index >= cols - 1)
			throw "访问越界!\n";
	}
	catch (char *a) {
		cout << a;
	}
	return *(p + index);
}

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);
	if (A.p[A.rows - 1][A.rows - 1] == 0 && A.p[A.rows - 1][A.cols - 1] != 0) assert(0);//有一行系数矩阵全为0且常数项不为0则无解
	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;
	E.Normalizing(rows, cols);
	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++) {
		if (A.p[i][i] == 0) assert(0);//无逆矩阵
		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;
}
Mat Mat::matPow(int n)const {
	/*快速幂乘*/
	if (cols != rows) assert(0);
	Mat t(rows, cols), a = *this;
	t.Normalizing(rows, cols);
	while (n) {
		if (n & 1)
			t *= a;
		n >>= 1;
		a *= a;
	}
	return t;
}
void Mat::Normalizing(int r, int c) {
	if (r != c) assert(0);
	*this = Mat(r, c, 0);
	for (int i = 0; i < r; i++)
		p[i][i] = 1;
}
int Mat::row()const {
	return rows;
}
int Mat::col()const {
	return cols;
}
void Mat::swapLine(int r1, int r2) {
	double *t = p[r1];
	p[r1] = p[r2];
	p[r2] = 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';
	}
}
istream &operator >> (istream &in, Mat &arr) {
	for (int i = 0; i < arr.rows; i++) {
		for (int j = 0; j < arr.cols; j++) {
			in >> arr.p[i][j];
		}
	}
	return in;
}
ostream &operator << (ostream &out,const Mat &arr) {

	for (int i = 0; i < arr.rows; i++) {
		for (int j = 0; j < arr.cols; j++)
			out << arr.p[i][j] << ' ';
		out << '\n';
	}
	return out;
}

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;
}

更新 8.9

  • 6
    点赞
  • 34
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值