MATLAB的矩阵可以有任意多维，在矩阵只有两维时可以进行乘法运算，以及进行奇异值分解svd等运算。用C或C++的数组实现矩阵运算，很难模拟矩阵的块赋值运算，即只对矩阵中某一方块中的元素（可能隔行或隔列）进行赋值。此外，MATLAB矩阵在进行sum及区块等运算时会自动降维，其结果仍然是一个矩阵，MATLAB其它矩阵运算还包括矩阵的拼接等。
如果用一系列函数去实现上述运算，就不如MATLAB的表达简洁和美观。当然，利用C++的类和运算符重载，可以使矩阵运算显得简洁美观。但是，模仿MATLAB的矩阵运算仍然存在许多障碍，例如，用矩阵A进行奇异值分解可得三个矩阵：[S, V, D] = svd(A)，以往的C++很难对赋值进行模仿，因为C和C++的函数返回值只能有一个值。另外，MATLAB矩阵中的块可以被赋值，只是改变了原矩阵的一部分值；而用C++的类重新构造一个矩阵块，则矩阵块与被赋值的原矩阵无关。
因此，用C++的类模拟MATLAB的矩阵存在许多问题，许多问题与C++语言本身的支持程度相关，而有些问题是可以通过一定的设计技巧克服的。模仿MATLAB存在的主要问题如下：
(1) 如何定义任意多维数组？
(2) 如果用类模仿，需要为二维、三维等定义多个类型吗？如何用一个类表示任意多维数组？
(3) 如何使块赋值做出的修改影响原有矩阵的部分值？
(4) 如何解决奇异值分解函数不能返回多个值赋给3个变量的问题？
(5) 如何表达某一维的上下界以及步长，例如MATLAB的A[1:6, 1:6]以支持取块等操作？
还可以列出一些模仿MATLAB矩阵存在的问题。先谈第(1)个问题的解决方案，如果采用变参类模板，随着模板使用实参类型进行展开，实际上会产生许多重载的构造函数。因此，在定义类模板时不建议采用递归展开的形式，否则会降低编译速度且使程序冗长。
第(2)个问题解决方案是采用一维数组模拟多维数组，为此需要同时记录每一维的大小。第(3)个问题的解决首先涉及维的上下界描述，甚至涉及维的步长，与第(5)个问题密切相关；其次涉及取块得到的矩阵(构造的新对象）元素如何同原矩阵关联，解决方案是设置一个原始矩阵标志，取块得到的矩阵元素（指针类型）与原始矩阵值相同，如此可以在修改矩阵元素的值时修改原始矩阵的元素；析构时若原始矩阵标志为假，不释放块矩阵中指针指向的元素内存。
第(4)个问题的解决则依赖于C++引入的结构化类型推导，例如对于矩阵a使用auto [s, v, d] = a.svd()。第(5)个问题的解决则依赖于C++的类模板initializer_list，如此则MATLAB的A[1:6, 1:6]可用A[{1,6}][{1,6}]模拟，其中{1,6}的类型为initializer_list<int>实例类，若{...}有三个元素则可以描述维的步长。对于MATLAB的维用begin等符号，可用initializer_list的实例类initializer_list<const char*>表示。
如此，利用《C++程序设计精要教程》的C++最新标准，可定义超级矩阵类MAT，并模拟实现MATLAB矩阵运算的功能。头文件MAT.h用于说明类型信息，MAT.cpp用于定义相关超级矩阵类的函数成员。如下为MAT.h，可根据需要扩展该定义更多的矩阵运算函数。
#pragma once
#include <list>
#include <typeinfo>
#include <iostream>
using namespace std;
class DIM { //定义维类
int* d, n;
public:
DIM(int x);
DIM(const std::list<int>&); //定义DIM(const std::list<const char *>&)，可支持begin等表示
DIM(const DIM&);
DIM(DIM&&);
DIM& operator=(const DIM&);
DIM& operator=(DIM&&);
int& operator[](int x);
operator int()const;
~DIM();
};
template <typename T> struct SVD;
template <typename T>
class MAT {
T** const e;//存放元素
long z; //e的元素个数
int *d; //矩阵各维的界
const int n; //维数
bool o; //是否原始矩阵
public:
MAT();
template<typename...Args>
MAT(int, Args...args)requires conjunction_v<is_same<int, Args>...>;
MAT(const MAT&);
MAT(MAT&&);
operator T& ();
MAT operator[](int x);
MAT operator[](std::list<DIM> x);
MAT operator+(const MAT&)const;
MAT operator-(const MAT&)const;
MAT operator*(const MAT&)const;
MAT& operator=(const T&);
MAT& operator=(T&&);
MAT& operator=(const MAT&);
MAT& operator=(MAT&&);
MAT sum(int x)const;
SVD<T> svd( )const;
void print();
~MAT();
};
template <typename T>
struct SVD { MAT<T> u, s, v; };
//template <typename T, typename...Args>
//SVD<T, Args...> svd(const MAT<T, Args...>&); //可以加上auto [U,S,V]=svd(A)
以下是MAT.cpp文件定义的相关运算或函数。注意，构造函数的可变形参没有采用递归展开的方式进行定义：
#include "MAT.h"
#include <type_traits>
DIM::DIM(int x):d(new int(x)),n(d?1:0) {
if (d == nullptr) throw "Memory allocation error when constructing DIM!";
}
DIM::DIM(const std::list<int>& x): d(nullptr), n(0) {
int c = x.size();
if (c != 2 && c != 3) throw "Boundary format error!";
if (c == 2) {
int f = x.front();
int l = x.back();
if(f>=l) throw "Boundary format error!";
n = l - f + 1;
d = new int[n];
for (int s = 0; s < n; f++, s++) d[s] = f;
return;
}
int t[3], m=0;
for (int s: x) {
t[m] = s; m++;
}
if (t[1] == 0) throw "Boundary format error!";
if (t[1] > 0 && t[0] >= t[2]) throw "Boundary format error!";
if (t[1] < 0 && t[0] <= t[2]) throw "Boundary format error!";
if (t[1] > 0) {
n = (t[2] - t[0]) / t[1]+1;
d = new int[n];
for (int s = 0; s < n; t[0]+=t[1], s++) d[s] = t[0];
return;
}
n = (t[2] - t[0]) / t[1] + 1;
d = new int[n];
for (int s = 0; s < n; t[0] += t[1], s++) d[s] = t[0];
}
DIM::DIM(const DIM& a) : d(new int[a.n]), n(d ? a.n : 0) {
if(d == nullptr) throw "Memory allocation error when constructing DIM!";
for (int x = 0; x < n; x++) d[x] = a.d[x];
}
DIM::DIM(DIM&& a) : d(a.d), n(a.n) {
a.d = nullptr;
a.n = 0;
}
DIM& DIM::operator=(const DIM&a) {
if(this==&a) return *this;
if (d) delete d;
d = new int[a.n];
if (d == nullptr) throw "Memory allocation error when constructing DIM!";
for (n = 0; n < a.n; n++) d[n] = a.d[n];
return *this;
}
DIM& DIM::operator=(DIM&&a) {
if (this == &a) return *this;
if (d) delete d;
d = a.d;
n = a.n;
a.d = nullptr;
a.n = 0;
return *this;
}
int& DIM::operator[](int x) {
if(x<0 || x>=n) throw "Subscriptin error for DIM!";
return d[x];
}
DIM::operator int()const {
return this->n;
}
DIM::~DIM() {
if (d) {
delete d; d = nullptr; n = 0;
}
}
template <typename T>
MAT<T>::MAT(): e(nullptr), z(0), d(nullptr), n(0), o(false) { }
template <typename T>
template <typename ...Args>
MAT<T>::MAT(int f, Args...args)requires conjunction_v<is_same<int, Args>...> : e(nullptr), n(0), o(true) {
static_assert(conjunction_v<is_same<int, Args>...>, "dimensions must be integers");
int s = sizeof...(Args);
//bool b = conjunction_v<is_same<int, Args>...>;
int* p = 1 + &f; //第2维地址
(int&)n = 1 + s; //维数
d = new int[n]; //为各维的界分配空间
if (d == nullptr) throw "Memory allocation error when constructing MAT!";
z = d[0] = f;
for (int m = 0; m < s; m++)
z *= d[m + 1] = p[m];
(T**&)e = new T * [z] {nullptr};
if (e == nullptr) throw "Memory allocation error when constructing MAT!";
for (long x = 0; x < z; x++) {
e[x] = new T{ 0 };
if (e[x] == nullptr) throw "Memory allocation error when constructing MAT!";
}
}
template <typename T>
MAT<T>::MAT(const MAT& a): e(new T* [a.z]), z(e ? a.z : 0), d(new int[a.n]), n(d ? a.n : 0), o(true) {
if (e == nullptr) throw "Memory allocation error when constructing MAT!";
if (d == nullptr) throw "Memory allocation error when constructing MAT!";
for (int m = 0; m < n; m++) d[m] = a.d[m];
for (long x = 0; x < z; x++) {
e[x] = new T{ 0 };
if (e[x] == nullptr) throw "Memory allocation error when constructing MAT!";
*e[x] = *a.e[x];
}
}
template <typename T>
MAT<T>::MAT(MAT&& a) : e(a.e), n(a.n), d(a.d), z(a.z), o(a.o) {
(T**&)(a.e) = nullptr;
a.z = 0;
a.d = nullptr;
(int&)(a.n) = 0;
a.o = false;
}
template <typename T>
MAT<T>::~MAT() {
if (e == nullptr) return;
if (o == true) {
for (long x = 0; x < z; x++)
if (e[x] != nullptr) delete e[x];
}
delete[]e;
if (d != nullptr) delete[]d;
(T**&)e = nullptr;
z = 0;
d = nullptr;
(int &)n = 0;
o = false;
}
template <typename T>
MAT<T> MAT<T>::operator[](int x) {
if (n < 1) throw "decreasing dimension error!";
if (x < 0 || x >= d[0]) throw "subscription overflow!";
MAT<T> r;
r.z = z / d[0];
(int&)(r.n) = n - 1;
(T**&)(r.e) = new T * [r.z]{ nullptr };
if (r.e == nullptr) throw "Memory allocation error when decresing dimension!";
long p = x * r.z;
for (long h = 0; h < r.z; h++) r.e[h] = e[p + h];
r.d = new int[r.n];
if (r.d == nullptr) throw "Memory allocation error when decresing dimension!";
for (int m = 1; m < n; m++) r.d[m-1] = d[m];
r.o = false;
return r;
}
template <typename T>
MAT<T> MAT<T>::operator[](std::list<DIM> x) {
if (n < 1) throw "decreasing dimension error!";
int t = 0;
for (DIM m : x) {
t+=m.operator int();
for(int f=0, k=m; f<k; f++)
if (m[f] < 0 || m[f] >= d[0]) throw "subscription overflow!";
}
MAT<T> r;
r.z = z / d[0] * t;
(int&)(r.n) = n;
(T**&)(r.e) = new T * [r.z]{ nullptr };
r.d = new int[r.n];
if (r.d == nullptr) throw "Memory allocation error when decresing dimension!";
r.d[0] = t;
for (int m = 1; m < n; m++) r.d[m] = d[m];
long rz = z / d[0];
t = 0;
for (DIM m : x) {
for (int f = 0, k= m; f < k; f++) {
long p = m[f] * rz;
long q = t * rz;
for (long h = 0; h < rz; h++) r.e[q + h] = e[p + h];
t++;
}
}
r.o = false;
return r;
}
template <typename T>
MAT<T>::operator T& () {
if(z!=1) throw "The number of elements must be 1 when fetching element!";
return *e[0];
}
template <typename T>
MAT<T> MAT<T>::operator+(const MAT& a)const {
if (z != a.z) throw "The size is not equal to each other when adding two MATs!";
if (n != a.n) throw "The number of dimensions is not equal to each other when adding two MATs!";
for (int m = 0; m < n; m++) {
if (d[m] != a.d[m]) throw "The size of each dimension is not equal to each other when adding two MATs!";
}
if (typeid(**e)!=typeid(**a.e)) "The type of element is not equal to each other when adding two MATs!";
MAT<T> r(a); //r.o=true即是分配元素内存的
for (long h = 0; h < z; h++) (*r.e[h]) += *e[h];
return r;
}
template <typename T>
MAT<T> MAT<T>::operator-(const MAT& a)const {
if (z != a.z) throw "The size is not equal to each other when substracting two MATs!";
if (n != a.n) throw "The number of dimensions is not equal to each other when substracting two MATs!";
for (int m = 0; m < n; m++) {
if (d[m] != a.d[m]) throw "The size of each dimension is not equal to each other when substracting two MATs!";
}
if (typeid(**e) != typeid(**a.e)) "The type of element is not equal to each other when substracting two MATs!";
MAT<T> r(*this); //r.o=true即是分配元素内存的
for (long h = 0; h < z; h++) (*r.e[h]) -= *a.e[h];
return r;
}
//可实现一个较“通用”的矩阵乘法：检查n和a.n==2，对返回结果可定义MAT<T, Args...> r()，重新为r分配内存
MAT<int> MAT<int>::operator*(const MAT& a)const {
if (n != 2 || a.n != 2) throw "The number of dimensions must be 2 when multiplying two MATs!";
if (d[1] != a.d[0]) throw "The size of second dimension of first matrix is not equal to the size of first dimension of second matrix when multiplying two MATs!";
if (typeid(**e) != typeid(**a.e)) "The type of element is not equal to each other when multiplying two MATs!";
MAT<int> r(d[0], a.d[1]); //r.o=true即是分配元素内存的
for (int x = 0; x < d[0]; x++)
for (int y = 0; y < a.d[1]; y++) {
long p = x * a.d[1] + y;
*r.e[p] = 0;
for (int z = 0; z < d[1]; z++) {
*r.e[p] += *e[x * d[1] + z] * *a.e[z * a.d[1] + y];
}
}
return r;
}
MAT<double> MAT<double>::operator*(const MAT& a)const {
if (n != 2 || a.n != 2) throw "The number of dimensions must be 2 when multiplying two MATs!";
if (d[1] != a.d[0]) throw "The size of second dimension of first matrix is not equal to the size of first dimension of second matrix when multiplying two MATs!";
if (typeid(**e) != typeid(**a.e)) "The type of element is not equal to each other when multiplying two MATs!";
MAT<double> r(d[0], a.d[1]); //r.o=true即是分配元素内存的
for (int x = 0; x < d[0]; x++)
for (int y = 0; y < a.d[1]; y++) {
long p = x * a.d[1] + y;
*r.e[p] = 0;
for (int z = 0; z < d[1]; z++) {
*r.e[p] += *e[x * d[1] + z] * *a.e[z * a.d[1] + y];
}
}
return r;
}
template <typename T>
MAT<T>& MAT<T>::operator=(const T& t) {
if (z != 1) throw "The number of elements must be 1 when fetching element!";
*e[0] = t;
return *this;
}
template <typename T>
MAT<T>& MAT<T>::operator=(T&& t) {
if (z != 1) throw "The number of elements must be 1 when fetching element!";
*e[0] = t;
return *this;
}
template <typename T>
MAT<T>& MAT<T>::operator=(const MAT& a) {
if (this == &a) return *this;
if (o == false) {//赋值给非原生MAT：两个MAT的维的界必须相同
if (z != a.z) throw "The size is not equal to each other when assigning an matrix!";
if (n != a.n) throw "The number of dimensions is not equal to each other when assigning an matrix!";
for (int m = 0; m < n; m++) {
if (d[m] != a.d[m]) throw "The size of each dimension is not equal to each other when assigning an matrix!";
}
if (typeid(**e) != typeid(**a.e)) "The type of element is not equal to each other when assigning an matrix!";
for (long h = 0; h < z; h++) (*e[h]) = *a.e[h];
}
else { //赋值给原生MAT：两个MAT的维的界可以相同
if (e != nullptr && a.e != nullptr && typeid(**e) != typeid(**a.e)) "The type of element is not equal to each other when assigning an matrix!";;
if (e) for (long x = 0; x < z; x++) if (e[x] != nullptr) delete e[x];
delete[]e;
delete[]d;
(T**&)e = new T * [z = a.z];
if (e == nullptr) throw "Memory allocation error when assigning an matrix!";
for (long x = 0; x < z; x++) {
e[x] = new T{ 0 };
if (e[x] == nullptr) throw "Memory allocation error when assigning an matrix!";
*e[x] = *a.e[x];
}
d = new int[(int&)n = a.n];
if (d == nullptr) throw "Memory allocation error when assigning an matrix!";
for (int m = 0; m < n; m++) (int&)(d[m]) = int(a.d[m]);
}
return *this;
}
template <typename T>
MAT<T>& MAT<T>::operator=(MAT&& a) {
if (this == &a) return *this;
if (o && e) for (long x = 0; x < z; x++) if (e[x] != nullptr) delete e[x];
if(e) delete[]e;
if(d) delete[]d;
(T**&)e = a.e;
z = a.z;
d = a.d;
(int &)n = a.n;
o = a.o;
(T**&)(a.e) = nullptr;
a.z = 0;
a.d = 0;
(int&)(a.n) = 0;
a.o = false;
return *this;
}
template <typename T>
SVD<T> MAT<T>::svd( )const {
SVD<T> s;
return s;
};
template <typename T>
MAT<T> MAT<T>::sum(int x)const {
MAT<T> m;
return m;
}
template <typename T>
void MAT<T>::print( ){
if (n != 2) throw "Not a 2 dimension matrix, can not print out!";
for (long x = 0; x < d[0]; x++) {
for (long y = 0; y < d[1]; y++) {
cout << "\t" << *e[x * d[1] + y];
}
cout << endl;
}
}
template MAT<int>;
template MAT<int>::MAT(int);
template MAT<int>::MAT(int, int);
template MAT<int>::MAT(int, int, int);
template MAT<double>;
template MAT<double>::MAT(int);
template MAT<double>::MAT(int, int);
template MAT<double>::MAT(int, int, int);
可以根据自己的需要，扩展定义更多的矩阵运算重载函数。在定义好上述类模板以后，可以添加如下主函数进行测试:
#include "MAT.h"
int main()
{
MAT<int> m(3);
MAT<int> a(3, 3);
MAT<int> e(2, 3);
MAT<double> n(3, 3, 3);
string st = typeid(MAT<int>).name();
st = typeid(MAT<double>).name();
a[0][0] = 1;
a[0][1] = 2;
a[0][2] = 3;
a[1][0] = 4;
a[1][1] = 5;
a[1][2] = 6;
a[2][0] = 7;
a[2][1] = 8;
a[2][2] = 9;
auto [s, v, d] = a.svd();
std::list<int> b = { 0,2 };
e = a[{0, DIM({ 1,2 })}];
e[0][0] = 11;
int z = e[0][0];
e = a[{ 1, 2 }]; //也可以对块a[{1,2}]赋值，修改原始类的部分元素。
e = e * a;
e.print();
}
1858
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
