行主映射:将矩阵依此从每一行的索引从左至右连续编号，将二维矩阵的索引映射为一维数组。
列主映射:将矩阵依此从每一列的索引从上至下连续编号，将二维矩阵的索引映射为一维数组。
本文的多种矩阵均采用行主映射，使用一维数组来表达多种矩阵。例如按照常规对二维矩阵进行描述的C++二维数组中，矩阵元素M(3,5)对应二维数组Array[2][4]。而采用一维映射后，该元素对应一维数组的Array[3*矩阵的列+5-1]。这样映射的好处在于节省内存空间，提高一些操作(如转置，矩阵乘法)的运行效率，虽然对于一般矩阵的提升不够明显，但是对于类似一些特殊矩阵(如对角矩阵、三对角矩阵、稀疏矩阵)所带来的效果极其明显。
对于特殊矩阵，需要重新定义数据结构进行存储。常见的一种方法是用三元组表达(row,col,value)，以节省空间并提高效率。同时，因为与常规矩阵结构的不同，其加减乘除等操作需要重新编写算法。
本文中会用到以前编写的自定义异常类、数组线性表等程序。见线性表–数组/vecotr描述。
一般矩阵(m*n):
实现矩阵的加/减/乘运算和矩阵类的各种拷贝控制成员以及其它一些读写操作。具体思想以及其他一些思考见注释。
#pragma once
#include <iostream>
#include <algorithm>
#include <iomanip>
#include "../../../ch06/ChainList/ChainList/illegalParameterValue.h"
using std::cin; using std::cout; using std::ends; using std::endl;
//前置声明
template <typename T> class myMatrix;
//将operator<<声明为模板 使用myMatrix的模板形参作为自己的模板实参 限定友好关系
template <typename T> std::ostream& operator<<(std::ostream& os,const myMatrix<T>& matrix);
template <typename T>
class myMatrix{
friend std::ostream& operator<<<T>(std::ostream& os,const myMatrix<T>& matrix);
public:
myMatrix()= default;
myMatrix(size_t _row,size_t _col);
myMatrix(const myMatrix<T>& rhs);//copy constructor
myMatrix(myMatrix<T>&& rhs);//move construcotr
~myMatrix() { delete element; }
size_t get_row()const { return row; }
size_t get_col()const { return col; }
T& operator()(size_t _row,size_t _col) const;//重载()来表示矩阵行列的索引 并完成矩阵数学逻辑上的索引到内存数组索引的映射
void set(size_t _row,size_t _col,const T& _element);//改变矩阵(i,j)元素
void reset();//将矩阵全部元素复位为0
myMatrix<T>& tranpose();//将矩阵转置
myMatrix<T>& operator=(const myMatrix<T>& rhs);//copy assignment
myMatrix<T>& operator=(myMatrix<T>&& rhs);//move assignment
myMatrix<T> operator+(const myMatrix<T>& rhs) const;//矩阵加法
myMatrix<T> operator-(const myMatrix<T>& rhs) const;//矩阵减法
myMatrix<T> operator*(const myMatrix<T>& rhs) const;//矩阵乘法
myMatrix<T> operator+() const;//一元操作符+ 取正
myMatrix<T> operator-() const;//一元操作符- 取负
myMatrix<T>& operator+=(const myMatrix<T>& rhs);//复合赋值运算符
myMatrix<T>& operator+=(size_t _unit);//复合赋值运算符 该矩阵加上一个数_unit(即_unit个单位矩阵)
myMatrix<T>& operator-=(const myMatrix<T>& rhs);//复合赋值运算符
myMatrix<T>& operator-=(size_t _unit);//复合赋值运算符 该矩阵减去一个数_unit(即_unit个单位矩阵)
myMatrix<T>& operator*=(size_t _unit);//复合赋值运算符 该矩阵乘以一个数_unit(即_unit个单位矩阵)
private:
size_t row, col;//矩阵的行列
T* element;//数组element
void check_if_same_size(size_t rhs_col,size_t rhs_row)const;//检查2个矩阵是否行可加减(行数列数相同)
void check_if_multiplicable(size_t rhs_row) const;//检查2个矩阵是否可相乘(列数=行数)
};
template <typename T>
myMatrix<T>::myMatrix(size_t _row,size_t _col):row(_row),col(_col) {
//检验矩阵行列参数的有效性
if (row < 0 || col < 0) throw illegalParameterValue("Rows and columns must be>=0!");
//行列不能只有其中一个为0(可以同时为0 生成0*0的矩阵)
if ((row == 0 || col == 0) && (row != 0 || col != 0))
throw illegalParameterValue("Either both or neither rows and columns should be 0!");
element = new T[row*col];
}
template <typename T>
myMatrix<T>::myMatrix(const myMatrix<T>& rhs):row(rhs.row),col(rhs.col) { //拷贝构造函数
element = new T[row*col];
//利用std::copy拷贝rhs的元素
//注意 copy(start,end,target)的拷贝范围是[start,end) 所以end要比元素个数多1
std::copy(rhs.element, rhs.element + rhs.row*rhs.col, element);
}
template <typename T>
myMatrix<T>::myMatrix(myMatrix<T>&& rhs) :row(rhs.row),col(rhs.col),element(rhs.element){ //移动构造函数
//直接接管资源 然后将rhs置于可析构状态
cout << "move constructor" << endl;
rhs.element = nullptr;
}
template <typename T>
T& myMatrix<T>::operator()(size_t _row,size_t _col) const{
//检查序列(row,col)是否合法
if (_row<1 || _col<1 || _row>row || _col>col)
throw illegalParameterValue("i or j of matrix(i,j) out of range");
return element[(_row-1)*col+(_col-1)];//数学逻辑到物理地址的映射
}
template <typename T>
void myMatrix<T>::set(size_t _row,size_t _col,const T& _element) {
operator()(_row,_col) = _element;
}
template <typename T>
void myMatrix<T>::reset() { //将所有元素置为0
for (size_t i = 0; i < row*col; ++i)
element[i] = 0;
}
template <typename T>
void myMatrix<T>::check_if_same_size(size_t rhs_row,size_t rhs_col)const {
if (row != rhs_row || col != rhs_col)
throw illegalParameterValue("rows and columns of two matrixs must be the same!");
}
template <typename T>
void myMatrix<T>::check_if_multiplicable(size_t rhs_row)const {
if (col != rhs_row)
throw illegalParameterValue("columns of left matrix and rows of right matrix must be the same to multiply!");
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator=(const myMatrix<T>& rhs) { //拷贝赋值运算符
//先拷贝再释放内存避免自赋值
row = rhs.row;
col = rhs.col;
T* newdata = new T[row*col];
std::copy(rhs.element, rhs.element + row*col, newdata);
delete[] element;
element = newdata;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator=(myMatrix<T>&& rhs) { //移动赋值运算符
//直接接管资源
if (this!=&rhs) { //避免自赋值
row = rhs.row;
col = rhs.col;
element = rhs.element;
}
//将rhs置于可析构状态
rhs.col = rhs.row = 0;
rhs.element = nullptr;
cout << "move assignment" << endl;
return *this;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator+=(const myMatrix<T>& rhs){ //复合赋值运算符+=
check_if_same_size(rhs.row,rhs.col);
for (size_t i = 0; i < row*col;++i) {
element[i] += rhs.element[i];
}
return *this;
}
template <typename T>
myMatrix<T> myMatrix<T>::operator+(const myMatrix<T>& rhs) const { //矩阵加法
//检查是否两个矩阵行数列数相同
//Q:是否可以使用+=来实现+? A:会影响效率.因为涉及到对数组的2次拷贝 考虑以下一种实现情况:
//myMatrix<T> retMatrix(*this);
//retMatrix += rhs;
//return retMatrix;
//可以看到在拷贝构造*this是对element数组进行了一次for遍历 在+=时又对element数组进行了一次for遍历
//因此 其他一般用+=运算符定义+运算符的情况在这里不适用 -= *=下同
check_if_same_size(rhs.row,rhs.col);
myMatrix<T> retMatrix(row, col);//返回的矩阵
for (size_t i = 0; i < row*col; ++i) {
retMatrix.element[i] = element[i] + rhs.element[i];
}
return retMatrix;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator+=(size_t _unit) { //复合赋值运算符+= 矩阵加上一个数_unit(即_unit个单位矩阵)
for (size_t i = 0; i < row*col;++i) {
element[i] += _unit;
}
return *this;
}
template <typename T>
myMatrix<T> myMatrix<T>::operator-(const myMatrix<T>& rhs) const{ //矩阵减法
check_if_same_size(rhs.row,rhs.col);
myMatrix<T> retMatrix(row,col);
for (size_t i = 0; i < row*col;++i) {
retMatrix.element[i] = element[i] - rhs.element[i];
}
return retMatrix;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator-=(const myMatrix<T>& rhs) { //复合赋值运算符-=
check_if_same_size(rhs.col,rhs.row);
for (size_t i = 0; i < col*row;++i) {
element[i] -= rhs.element[i];
}
return *this;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator-=(size_t _unit) { //复合赋值运算符-= 矩阵减去一个数_unit(即_unit个单位矩阵)
for (size_t i = 0; i < row*col; ++i) {
element[i] -= _unit;
}
return *this;
}
template <typename T>
myMatrix<T> myMatrix<T>::operator+() const { //矩阵所有元素取正
myMatrix<T> retMatrix(*this);
for (size_t i = 0; i < row*col; ++i) {
retMatrix.element[i] = +retMatrix.element[i];
}
return retMatrix;
}
template <typename T>
myMatrix<T> myMatrix<T>::operator-() const { //矩阵所有元素取负
myMatrix<T> retMatrix(*this);
for (size_t i = 0; i < row*col; ++i) {
retMatrix.element[i] = -retMatrix.element[i];
}
return retMatrix;
}
template <typename T>
myMatrix<T> myMatrix<T>::operator*(const myMatrix<T>& rhs) const { //矩阵乘法
check_if_multiplicable(rhs.row);
myMatrix<T> retMatrix(row,rhs.col);//结果矩阵是row行rhs.col列
//注意 默认矩阵到一维数组是行主映射
size_t row_ptr = 0;//用于遍历乘号左边矩阵的当前行的每一个元素
size_t col_ptr = 0;//用于遍历乘号右边矩阵的当前列的每一个元素
size_t arr_ptr = 0;//用于遍历被映射的一维数组
for (size_t i = 0; i < row;++i) { //结果矩阵的每一行
for (size_t j = 0; j < rhs.col;++j) { //结果矩阵的每一列
T sum = 0;
for (size_t k = 0; k < col;++k) { //结果矩阵的每一个元素matrix(i,j)=左边矩阵的第i行乘以右边矩阵的第j列
sum += element[row_ptr] * rhs.element[col_ptr];
++row_ptr;//移动到左边矩阵当前行的下一个元素
col_ptr += rhs.col;//移动到右边矩阵当前列的下一个元素
}
retMatrix.element[arr_ptr++] = sum;
row_ptr -= col;//后退到左边矩阵当前行的行首元素
col_ptr = j+1;//移动到右边矩阵的下一列的首元素
}
row_ptr += col;//移动到左边矩阵的下一行的首元素
col_ptr = 0;//后退到右边矩阵的第一列的首元素
}
return retMatrix;
}
template <typename T>
myMatrix<T>& myMatrix<T>::operator*=(size_t _unit) { //复合赋值运算符*= 矩阵乘以一个数_unit(即_unit个单位矩阵)
for (size_t i = 0; i < row*col; ++i) {
element[i] *= _unit;
}
return *this;
}
template <typename T>
myMatrix<T>& myMatrix<T>::tranpose() {
for (size_t i = 1; i <= row;++i) {
for (size_t j = i; j <= row;++j) {
std::swap(operator()(i,j),operator()(j,i));
}
}
return *this;
}
//重载友元operator<<操作符
template <typename T>
std::ostream& operator<<(std::ostream& os,const myMatrix<T>& matrix) {
for (size_t i = 1; i <= matrix.get_row(); ++i) {
cout << "| ";
for (size_t j = 1; j <= matrix.get_col(); ++j) {
cout << std::setw(5) << std::left << matrix(i, j) << ends;
}
cout << " |" << endl;
}
return os;
}
对于矩阵类的测试:
//======================普通矩阵myMatrix=======================
#include "myMatrix.h"
int main(){
myMatrix<int> matrix1(5,5);//构造一个5*5的矩阵
myMatrix<int> matrix2(5,5);//构造一个5*5的矩阵
//初始化
std::default_random_engine e;//不适用种子方便调试
std::uniform_int_distribution<int> u(0,10);
for (size_t i = 1; i <= 5;++i) {
for (size_t j = 1; j <= 5;++j) {
matrix1(i, j) = u(e);//随机赋值
}
}
for (size_t i = 1; i <= 5; ++i) {
for (size_t j = 1; j <= 5; ++j) {
matrix2(i, j) = u(e);//随机赋值
}
}
cout << matrix1 << endl;
cout << "================tranpose==============" << endl;
cout << matrix1.tranpose() << endl;
cout << matrix2 << endl;
cout << "=================operator+/move constructor==================" << endl;
myMatrix<int> matrix3(matrix1+matrix2);//matrix3调用移动构造函数
cout << matrix3 << endl;
cout << "===================operator-/move construcotr/move assignment======================" << endl;
matrix2=matrix3 - matrix1;//编译器对matrix3-matrix1的返回值进行优化 调用移动构造函数 matrix2=又调用了移动赋值运算符
cout << matrix2 << endl;
cout << "======================opearator+=/-=/*=/======================" << endl;
matrix3 += matrix2;
cout << matrix3 << endl;
matrix3 += 3;
cout << matrix3 << endl;
matrix3 -= 3;
cout << matrix3 << endl;
matrix3 -= matrix2;
cout << matrix3 << endl;
matrix3 *= 2;
cout << matrix3 << endl;
cout << "======================get_col/get_row/set/reset======================" << endl;
cout << "matrix3 row=" << matrix3.get_row() << " col=" << matrix3.get_col() << endl;
matrix3.reset();
matrix3.set(4,3,4);
cout << matrix3 << endl;
cout << "========================operator*=========================" << endl;
myMatrix<int> matrix4(3,4);
myMatrix<int> matrix5(4,5);
for (size_t i = 1; i <= matrix4.get_row(); ++i) {
for (size_t j = 1; j <= matrix4.get_col(); ++j) {
matrix4(i, j) = u(e);//随机赋值
}
}
for (size_t i = 1; i <= matrix5.get_row(); ++i) {
for (size_t j = 1; j <= matrix5.get_col(); ++j) {
matrix5(i, j) = u(e);//随机赋值
}
}
cout << matrix4 << endl << matrix5 << endl;
myMatrix<int> result_matrix(matrix4 * matrix5);
cout << result_matrix << endl;
}
对角矩阵(n*n):
对角矩阵只有n个元素，然而使用一般矩阵类来描述不仅浪费空间，矩阵操作也是效率低下，因此只存储对角线上的元素，完成从矩阵的数学描述到计算机数组存储的物理映射。部分实现，其他思考见注释。
#pragma once
//===================对角矩阵的部分实现===================
//将对角矩阵设计为仅仅对n阶对角阵中对角线上n个元素的存储
#include "../../../ch06/ChainList/ChainList/illegalParameterValue.h"
template <typename T>
class diagonalMatrix {
public:
diagonalMatrix() = default;
diagonalMatrix(size_t _n):n(_n),element(new T[n]()){ }
~diagonalMatrix() { delete element; }
T& get_element(size_t i, size_t j)const {
check(i,j);
if (i == j) return element[i - 1];//对角线上的元素
else return 0;
}
void set(size_t i, size_t j,const T& new_value) {
check(i,j);
//非对角元素只能为0
if (i == j) element[i - 1] = new_value;
else throw illegalParameterValue("nondiagonal elements must be 0");
}
private:
size_t n;//n阶对角矩阵
T* element;//对角元素的存储数组
void check() const{ //(i,j)索引是否有效
if (i<1 || i>n || j<1 || j>n)
throw illegalParameterValue("i or j must be >=1 and <=n!");
}
};
三对角矩阵(n*n):
类似对角矩阵。部分实现和其他思考见代码和注释。
#pragma once
//=================三对角矩阵的不完全实现=======================
//对于n阶三对角矩阵，可以只存储三条对角线的元素，共3n-2个
//采用逐条对角线映射
//也就是在数组element[]中先存储n-1个下对角线的元素，再存储n个主对角线的元素，最后存储n-1个上对角线的元素
#include "../../../ch06/ChainList/ChainList/illegalParameterValue.h"
template <typename T>
class tridiagonalMatrix {
public:
tridiagonalMatrix() = default;
tridiagonalMatrix(size_t _n):n(_n),element(new T[3*n-2]){ }
~tridiagonalMatrix() { delete element; }
T& get_element(size_t i,size_t j) const{
check(i,j);
switch (i-j){
case 1://下对角线
return element[i-2];//Q:i-2 ? A:下对角线的i是从2开始的，即(2,j) 而下对角线首先存储在element[0,n-1)中
case 0://主对角线
return element[n+i-2];//Q:n+i-2 ? A:(n-1)+(i-1)
case -1://上对角线
return element[2*n+i-2];//Q:2*n+i-2 ? A:(n-1)+(n)+(i-1)
default://非三对角线元素
return 0;
}
}
void set(size_t i,size_t k,const T& new_val) {
check(i,j);
switch (i-j){
case 1://下对角线
element[i - 2] = new_val;
break;
case 0://主对角线
element[n + i - 2] = new_val;
break;
case -1://上对角线
element[2 * n + i - 2] = new_val;
break;
default:
throw illegalParameterValue("nontridiagonal element must be 0!");//非三对角元素必须为0
}
}
private:
size_t n;//n阶矩阵
T* element;//存储三对角的元素
void check()const { //检查(i,j)索引是否有效
if (i<1 || i>n || j<1 || j>n)
throw illegalParameterValue("i or j must be >=1 and <=n!");
}
};
下三角矩阵(n*n):
类似地，一个set方法的实现，思考见代码和注释。
//===================类似地，下三角矩阵的映射公式========================
//构造函数中 element=new T[n*(n+1)/2]; 三角矩阵一共有(n+1)n/2个元素
template <typename>
void lowerTriangularMatrix<T>::set(size_t i,size_t j,const T& new_value) {
check(i,j);
if (i >= j) //在下三角的位置
element[n*(n-1)/2+j-1] = new_value;//第i行上面一共有i(i-1)/2个元素
else throw illegalParameterValue("elements not in lower triangle must be 0!");//非下三角元素必须为0
}
稀疏矩阵(m*n):
稀疏矩阵和稠密矩阵之间没有明确的界限。这里我们规律稀疏矩阵的非0元素要小于n*n/2，甚至有时要求n*n/5。对角阵和三对角矩阵的非0元素是O(n)，因此算稀疏矩阵，而三角矩阵的非0元素最多可达(n+1)n/2，因此算入稠密矩阵。
实现了稀疏矩阵的加法、乘法、转置操作，其他思考见代码和注释。
#pragma once
//=====================稀疏矩阵=============================
//稀疏矩阵和稠密矩阵之间没有明确的界限 这里我们规律稀疏矩阵的非0元素要小于n*n/2 甚至有时要求n*n/5
//对角阵和三对角矩阵的非0元素是O(n) 因此算稀疏矩阵 而三角矩阵的非0元素最多可达(n+1)n/2 因此算入稠密矩阵
#include "../../../ch05/ArrayList/ArrayList/myArrayList.h"
//本实现中 对稀疏矩阵的每个元素进行行主次序映射
//======================稀疏矩阵每个元素的数据结构========================
template <typename T>
struct matrixTerm{
size_t row;
size_t col;
T value;
//因为arrayList<T>中的index_of/output等函数要求T实现==符号 <<符号
bool operator==(const matrixTerm<T>& rhs)const {
return row == rhs.row && col == rhs.col && value == rhs.value;
}
};
template <typename T>
std::ostream& operator<<(std::ostream& os,const matrixTerm<T>& term) {
os << "a(" << term.row << "," << term.col << ")=" << term.value << endl;
return os;
}
//======================稀疏矩阵的类spareMatrix========================
//前置声明 将<< >>运算符定义为模板T 与spareMatrix<T>形成一对一友好关系
template <typename T> class spareMatrix;
template <typename T> std::ostream& operator<<(std::ostream& os,const spareMatrix<T>& matrix);
template <typename T> std::istream& operator>>(std::istream& is,spareMatrix<T>& matrix);
template <typename T>
class spareMatrix{
friend std::ostream& operator<<<T>(std::ostream& os, const spareMatrix<T>& matrix);
friend std::istream& operator>><T>(std::istream& is, spareMatrix<T>& matrix);
public:
//实际上 我们的spareMatrix类不需要自定义任何拷贝控制成员 (默认/非默认/拷贝/移动构造函数/析构函数)
//每一个拷贝构造成员都可以使用编译器合成的版本 它们也是安全合理的
//原因在于我们的数据成员中自定义类型成员myArrayList<T>它已经定义了自己相应(所有)的拷贝控制成员
//编译器在合成spareMatrix类会自动调用对应的myArrayList拷贝控制成员 以一个移动构造函数为例进行说明 其它成员同理:
//用于说明 operator+中编译器的优化:对于返回临时对象的移动构造
//spareMatrix(spareMatrix<T>&& rhs) :rows(rhs.rows),cols(rhs.cols),terms(std::move(rhs.terms)){
移动构造函数直接接管资源 调用myArrayList的移动构造函数
对于size_t这些内置类型不用移动构造(也没必要移动构造)
//cout << "move construcotr" << endl;
//}
spareMatrix<T> tranpose() const;
spareMatrix<T> add(const spareMatrix<T>& rhs) const;
spareMatrix<T> multiply(const spareMatrix<T>& rhs)const;
spareMatrix<T> operator+(const spareMatrix<T>& rhs)const {
//按照常理来说 我们这样编写代码会产生多个临时的中间变量对象的拷贝 然而实际上会被编译器优化掉
//编译器会对add()返回的临时对象进行优化 并不会直接拷贝 而是移动构造 可以见移动构造函数中的调试语句
return add(rhs);
}
spareMatrix<T> operator*(const spareMatrix<T>& rhs) const {
return multiply(rhs);
}
private:
size_t rows;//矩阵的行数
size_t cols;//矩阵的列数
myArrayList<matrixTerm<T>> terms;//非0项表
};
template <typename T>
std::ostream& operator<<(std::ostream& os,const spareMatrix<T>& matrix) {
//输出稀疏矩阵的基本信息
os << "Matrix rows=" << matrix.rows << " columns=" << matrix.cols << endl;
os << "nonzero terms=" << matrix.terms.get_size() << endl;
//每个项单独一行输出
for (auto iter = matrix.terms.begin(); iter != matrix.terms.end();++iter) {
os << *iter << endl;
}
return os;
}
template <typename T>
std::istream& operator>>(std::istream& is,spareMatrix<T>& matrix) {
size_t numbers_of_terms;
cout << "please enter number of rows,columns and terms:" << endl;
is >> matrix.rows >> matrix.cols >> numbers_of_terms;
if (numbers_of_terms > matrix.rows*matrix.cols / 2) //非0元素多余1/2 则不算入稀疏矩阵
throw illegalParameterValue("spare matrix requires number of nonzero elements being less than half!");
matrix.terms.resize(numbers_of_terms);//使用arraylist的resize 保证set(size)合法并且有足够的容量(设定为2*number)
//输入项
matrixTerm<T> perTerm;
for (size_t i = 0; i < numbers_of_terms;++i) {
cout << "Enter row,column and value of term(" << (i + 1) << ") in matrix:" << endl;
is >> perTerm.row >> perTerm.col >> perTerm.value;
if (perTerm.row<1 || perTerm.row>matrix.rows || perTerm.col<1 || perTerm.col>matrix.cols)
throw illegalParameterValue("terms index (i,j) out of matrix(row,col)");
if (perTerm.value == 0) throw illegalParameterValue("spare matrix save nonzero element1");
matrix.terms.set(i, perTerm);
}
return is;
}
template <typename T>
spareMatrix<T> spareMatrix<T>::tranpose() const { //返回稀疏矩阵的转置矩阵
//设置转置矩阵的特征
spareMatrix<T> tranMatrix;
tranMatrix.rows = cols;
tranMatrix.cols = rows;
tranMatrix.terms.resize(terms.get_size());
//计算转置前后元素的位置i 用于set(i,term) (注意 是按照行主次序映射的)
size_t* colSize = new size_t[cols];//记录转置前矩阵每一列的元素个数
size_t* newStart = new size_t[cols];//转置后每一行(原每一列)的元素的起始位置
//初始化
for(size_t i = 0; i < cols;++i) {
colSize[i] = 0;
newStart[i] = 0;
}
//遍历非0元素点 根据它们的列信息递增colSize[i]
for (auto iter = terms.begin(); iter != terms.end();++iter) {
++colSize[(*iter).col-1];//数学逻辑row,col需要减去1才能映射到数组
}
//寻找转置后矩阵每一行(原每一列)中各个元素的起始点
newStart[0] = 0;//转置后第一行(原第一列)的元素一定是从0开始的
for (size_t i = 1; i < cols;++i) {
newStart[i] = colSize[i - 1] + newStart[i - 1];
}
//开始转置
matrixTerm<T> perTerm;
for (auto iter = terms.begin(); iter != terms.end();++iter) {
size_t pos = newStart[(*iter).col-1]++;//当前行(原矩阵列)已经存储了一个元素了 下一个元素的起始点索引应该+1
perTerm.col = (*iter).row;//行列转置
perTerm.row = (*iter).col;
perTerm.value = (*iter).value;
//cout << "tranpose pos:" << pos;//调试语句 注意pos从数组取值时要对(*iter).colcol-1
tranMatrix.terms.set(pos,perTerm);
}
return tranMatrix;
}
template <typename T>
spareMatrix<T> spareMatrix<T>::add(const spareMatrix<T>& rhs)const {
if (cols != rhs.cols || rows != rhs.rows)
throw illegalParameterValue("rows and columns of two matrixs being added must be the same!");
spareMatrix<T> resMatrix;//结果矩阵
resMatrix.rows = rows;
resMatrix.cols = cols;
size_t res_size=0;//结果矩阵中元素的数量 也是将term insert到结果矩阵时的索引_index
auto lhs_iter = terms.begin();
auto rhs_iter = rhs.terms.begin();
//遍历*this和rhs 相同的项相加 不同的项直接插入结果矩阵 直到其中一个完成遍历
while (lhs_iter != terms.end() && rhs_iter!=rhs.terms.end()) {
//按照行主次序得到2个迭代器指代元素的相对位置索引
size_t lhs_index = (*lhs_iter).row*cols + (*lhs_iter).col;
size_t rhs_index = (*rhs_iter).row*cols + (*rhs_iter).col;
if (lhs_index<rhs_index) { //*this矩阵的项在前
resMatrix.terms.insert(*lhs_iter,res_size++);
++lhs_iter;
}
else if (lhs_index==rhs_index) { //两个项在同一位置
T add_sum = (*lhs_iter).value + (*rhs_iter).value;
if (add_sum != 0) { //两项相加后结果不为0则插入到结果矩阵
matrixTerm<T> term(*lhs_iter);//随便解引用lhs/rhs其中之一都行
term.value = add_sum;
resMatrix.terms.insert(term,res_size++);
}
++lhs_iter;
++rhs_iter;
}
else { //*this的项在rhs的后面
resMatrix.terms.insert(*rhs_iter,res_size++);
++rhs_iter;
}
}//其中一个(可能两个)矩阵遍历完成
//把剩下的另外一个矩阵的元素直接插入到结果矩阵
for (; lhs_iter != terms.end(); ++lhs_iter)
resMatrix.terms.insert(*lhs_iter,res_size++);
for (; rhs_iter != rhs.terms.end(); ++rhs_iter)
resMatrix.terms.insert(*rhs_iter,res_size++);
return resMatrix;
}
template <typename T>
spareMatrix<T> spareMatrix<T>::multiply(const spareMatrix<T>& rhs)const {
//基本思路是 首先构建lhs/rhs的位置矩阵 记录每一行第一个元素开始的索引(多个元素记得将索引++得到后面的元素)
//然后根据这个矩阵 对lhs每一行的每一个元素进行遍历 每一次遍历找到这个元素的列号(同时也是rhs中的行号)
//然后对rhs的这一行(注意是行!不按照常规的a[i][k]*b[k][j]来乘 因为我们可以将结果临时保存 每次遍历进行累加)
//依此乘以rhs行中每个元素 将每个乘积累加到一个大小为lhs.row/rhs.col的矩阵对应的位置中 lhs每行循环时记得将这个临时数组清空
//最重要的一点就是相乘时 并不是常规的lhs的行乘以rhs的列 而是lhs的行乘以对应rhs中不同的行 进行累加最后存入结果矩阵
if (cols != rhs.rows)
throw illegalParameterValue("columns of left matrix must be the same to rows of right matrix to multiply!");
spareMatrix<T> resMatrix;//结果矩阵
//设置结果矩阵的特征
resMatrix.rows = rows;
resMatrix.cols = rhs.cols;
size_t res_size = 0;//结果矩阵的非零元素个数 同时也是插入arraylist时的索引
if (terms.get_size() == 0 || rhs.terms.get_size() == 0)//保证都不是非零矩阵
return resMatrix;//全部元素为0
//构造lhs/rhs的位置矩阵 保存其每行元素在矩阵行主映射下起始的索引
size_t* lhs_row_start = new size_t[rows+1](); //创建时初始化
size_t* rhs_row_start = new size_t[rhs.rows+1](); //创建时初始化
//遍历 先得到每一行的元素个数
for (auto iter = terms.begin(); iter != terms.end(); ++iter) {
//第一行起始位置一定为0
++lhs_row_start[(*iter).row];
}
for (auto iter = rhs.terms.begin();iter!=rhs.terms.end();++iter) {
++rhs_row_start[(*iter).row];
}
//将每行大小数组转为每行首元素开始位置的索引
for (size_t s = 1; s < rows + 1;++s) {
lhs_row_start[s] += lhs_row_start[s - 1];
}
for (size_t s = 1; s < rhs.rows + 1; ++s) {
rhs_row_start[s] += rhs_row_start[s - 1];
}
//进行矩阵乘法
T* temArr = new T[rhs.cols]();//临时数组 用来保存每次乘积的累积和
for (size_t i = 0; i < rows;++i) { //对lhs的每一行
//将上一行的临时矩阵清空
for (size_t tem_col = 0; tem_col < rhs.cols;++tem_col) {
temArr[tem_col] = 0;
}
//构造临时矩阵
for (size_t j = lhs_row_start[i]; j < lhs_row_start[i + 1];++j) { //对lhs当前i行的每一个元素j
size_t rhs_row = terms.get_element(j).col;//找到该元素j的列数 也就是对应的rhs的行数
//对rhs矩阵中rhs_row行的起始位置到下一行的起始位置间的元素 注意col到数组的映射需要-1
for (size_t k = rhs_row_start[rhs_row - 1]; k < rhs_row_start[rhs_row];++k) {
//lhs的元素j与rhs中rhs_row行的每一个元素k相乘 累加在元素k对应的临时数组的列中
temArr[rhs.terms.get_element(k).col-1] += terms.get_element(j).value * rhs.terms.get_element(k).value;
}
}
//将临时矩阵的非零乘积结果存入结果矩阵
for (size_t tem_col = 0; tem_col < rhs.cols; ++tem_col) {
if (temArr[tem_col] != 0) {
matrixTerm<T> term;
term.row = i + 1;
term.col = tem_col + 1;
term.value = temArr[tem_col];
resMatrix.terms.insert(term,res_size++);//存入
}
}
}
return resMatrix;
}
对于稀疏矩阵的测试:
//======================稀疏矩阵=======================
#include "spareMatrix.h"
int main() {
//手动输入可能比较麻烦 建议做成文件输入 手动或文件输入矩阵特征后 每项格式都为 row col value
//测试合成的拷贝控制成员/输入运算符>>/输出运算符<<
//spareMatrix<int> m1;
//cin >> m1;
//cout << "======================================" << endl << m1 << endl;
//测试矩阵的转置
//spareMatrix<int> m2(m1.tranpose());
//cout << "======================================" << endl << m2 << endl;
//测试两个矩阵相加
//spareMatrix<int> m3;
//cin >> m3;
//spareMatrix<int> m4;
//cin >> m4;
//spareMatrix<int> m5 = m3 + m4;
//cout << "======================================" << endl << m5 << endl;
//spareMatrix<int>m6;
//cin >> m6;
//spareMatrix<int>m7(m6+m6);
//cout << m7;
//测试两个矩阵相乘
spareMatrix<int> m8;
spareMatrix<int> m9;
cin >> m8;
cin >> m9;
spareMatrix<int> m10(m8*m9);
cout << m10;
return 0;
}
