一个通用、可重用且可扩展的矩阵类
4.77/5 (36投票s)
2003年2月4日
41分钟阅读
293441
6652
一篇描述通用矩阵类及其实现方式的文章。
引言
有时，为了简化编程，需要将数据存储在矩阵中。这在科学编程中尤其如此，因为矩阵常用于计算。然而，C++ 中没有内置的 Matrix 类型。
在 C++ 中有几种实现矩阵的方法。最简单的方法是使用数组或指针。例如，要使用数组创建一个简单的矩阵，只需执行 double matrix[10][10]; 来创建一个具有 10 行 10 列的变量“matrix”。编程书籍会告诉你，数据并非真正以二维形式存储，但这与我们无关。重要的是我们可以轻松访问矩阵中的元素；例如，要访问第 2 行第 5 列的元素，我们执行 matrix[1][4]，考虑到数组从 0 开始。为便于解释，现在起，所有矩阵都从零开始，因此如果我提到第 2 行第 5 列，就是 matrix[2][5]。使用数组作为矩阵的缺点是行和列的大小必须在编译时固定。
要动态调整矩阵的大小，我们可以使用指针创建矩阵。我们可以使用
double* matrix = new double[100]
或
double** matrix = new double*[10]; for(int i=0; i<10;++i) matrix[i] = new double[10];
第一种形式是一种特殊的矩阵，要访问第 i 行第 j 列的元素，需要执行 matrix[i*10+j]。这显然不像 matrix[i][j] 那么直观。但是，它有自己的优点，即矩阵创建更简单。第二种形式创建矩阵稍微复杂一些，但访问元素与数组中的方式相同，即 matrix[i][j]。使用指针的缺点是，使用完矩阵后必须记住调用 delete[]，否则会导致内存泄漏。
创建矩阵的第三种方法是使用 STL 容器，如 vector。我们可以将两个 vector 嵌套在一起，例如 vector<vector<double> > matrix。要访问元素，我们仍然可以使用 matrix[i][j]。此外，我们可以通过 vector 容器的 size 函数相对轻松地获取矩阵的大小。我们还可以通过 resize 函数轻松地增加矩阵的大小。使用 STL 容器的另一个优点是它不像指针那样容易导致内存泄漏。
使用数组、指针或 STL 容器实现的矩阵进行数学运算比较困难。例如，要将两个矩阵相加，不能写成“matrixA + matrixB”。相反，您必须编写两个循环来逐个元素相加。为了方便计算并提供矩阵操作的有用函数，我们可以创建一个矩阵类。许多程序员都创建了自己的矩阵类。其中比较著名的是 C++ 科学计算（Numerical Recipes in C++）1 中实现的简单矩阵类，印第安纳大学开放系统实验室开发的矩阵模板库（MTL）2，以及 Masakatsu Ito3 实现的矩阵表达式模板。
既然已经有优秀的矩阵类可用，为什么还需要创建另一个呢？创建新矩阵类有几个动机，这也是本文的重点。首先，现有的矩阵类没有我需要的所有功能。向这些矩阵类添加额外功能需要理解类的实现，这很耗时。其次，作为一名业余程序员，创建新矩阵类将提高我的编程技能，因为我需要彻底考虑实现问题。第三个原因是，我不必花时间学习如何使用其他矩阵类。虽然有些很简单，但有些则不然。实现我自己的矩阵将消除这种学习曲线。
我实现的矩阵类并不完美。它并非旨在成为最好的矩阵类。在继续深入之前，先列出矩阵类的不足之处：
- 仅实现了密集矩阵存储类型。未实现稀疏矩阵、三角矩阵、对角矩阵等其他存储类型。
- 矩阵的性能不是现有矩阵类中最快的。尽管我没有进行广泛的比较测试，但我相信我的矩阵类在计算速度方面不是最快的。
- 用于操作矩阵的函数数量少于某些可用矩阵类。
- 最严重的不足是，此类只能由符合 C++ 标准的编译器编译。我只在 Visual C++ 2005 上测试过当前版本。
除了创建此类矩阵的动机之外，我的矩阵类相比其他矩阵类具有优势。否则，我不会发布此类的。其优点是：
- 该矩阵可用于存储数值类型和非数值类型。存储非数值类型时，数学函数不可用，如果意外使用，将导致编译器错误，而不是在运行时导致错误。
- 用于执行矩阵计算的函数易于使用，不像某些矩阵类那样。
- 尽管其他存储类型未实现，但用户可以通过遵循提供的
MatrixStorageTemplate来相对轻松地创建它们。 - 熟悉 STL 容器的用户可以使用类似 STL 的函数和 typedef。也可以在矩阵上执行某些 STL 算法。
对于那些认为此类可能适合他们，或者希望学习某些编程技术（如使用策略模板和表达式模板）的人来说，本文的其余部分分为以下几个部分。第一部分概述了可用的函数。第二部分举例说明了如何使用矩阵类。第三部分讨论了矩阵类的实现，详细介绍了如何使用策略进行泛型编程以及如何使用表达式模板来加速某些矩阵计算。
可用的矩阵函数
用户需要使用的唯一类是模板 Matrix<ValueT,MatrixTypeP,MathP> 类。ValueT 是矩阵中存储的数据类型的类型名称。它可以是 int、double 或 complex 等数值类型，也可以是 char* 等非数值类型，甚至是用户定义的类。MatrixTypeP 是数据存储类型的模板类策略。目前，只有 DenseMatrix 可用。MathP 是数学函数的模板类策略。有三种选择：MathMatrix、MathMETMatrix 和 NonMathMatrix。
无论创建何种类型的矩阵，以下内容对于所有类型的矩阵都是通用的：
类似 STL 的 typedef
|
|
存储在矩阵中的数据类型。 |
|
|
对存储在矩阵中的数据的引用。 |
|
|
对存储在矩阵中的数据的常量引用。 |
|
|
指向存储在矩阵中的数据的指针。 |
|
|
指向存储在矩阵中的数据的常量指针。 |
|
|
指针计算的差值类型。 |
|
|
大小类型。 |
构造函数、析构函数、复制构造函数和赋值运算符
|
|
默认构造函数 |
|
|
需要行数、列数和元素初始值的构造函数。 |
|
|
需要行数、列数以及包含元素初始值的 1D 数组的构造函数。 |
|
|
需要行数以及包含元素初始值的 2D 数组的构造函数。 |
|
|
析构函数。 |
|
|
复制构造函数。 |
|
|
接受任何可互换 |
|
|
接受同类型矩阵的赋值运算符。 |
|
|
接受任何可互换 |
|
|
将 1D 数组复制到矩阵。数组元素的数量必须等于或大于矩阵元素的数量，并且数组元素的顺序必须是连接起来以放入矩阵的行。 |
|
|
将 2D 数组复制到矩阵。数组的大小必须与矩阵相同。 |
矩阵成员函数
|
|
清空矩阵。 |
|
|
检查矩阵是否为空。 |
|
|
调整矩阵大小，用提供的或默认值填充新元素并删除多余元素。 |
|
|
交换当前矩阵与另一个矩阵的内容。要交换的矩阵必须是相同类型。 |
|
|
交换矩阵中的两行。 |
|
|
交换矩阵中的两列。 |
|
|
比较运算符。 |
|
|
更新矩阵。 |
|
|
成员访问函数 |
行和列迭代器函数
size_type size() const |
返回行/列的总数。 |
template<typename ForwardIterator> void insert(size_type rowNo, const ForwardIterator& first) |
在指定行/列处插入行/列。 |
void erase(size_type rowNo) |
删除指定行/列处的行/列。 |
template<typename ForwardIterator> void push_back(const ForwardIterator& first) |
在矩阵末尾添加行/列。 |
template<typename ForwardIterator> void push_front(const ForwardIterator& first) |
在矩阵开头添加行/列。 |
void pop_back() |
移除矩阵末尾的行/列。 |
void pop_front() |
移除矩阵开头的行/列。 |
|
|
获取矩阵数据迭代器的函数。 |
作用于矩阵的函数对象
|
|
获取矩阵的转置，即交换行和列。第一种形式返回转置后的矩阵。第二种形式将转置后的矩阵保存到传入的变量中。 |
|
|
获取矩阵的主对角线，或将向量置于矩阵的主对角线上。第一种形式返回对角矩阵。第二种形式将对角矩阵保存到传入的变量中。 |
|
|
获取矩阵的协方差。第一种形式返回协方差矩阵。第二种形式将协方差矩阵保存到传入的变量中。 |
|
|
获取一个矩阵，其中所有元素都提高到指定幂。第一种形式返回幂次矩阵。第二种形式将幂次矩阵保存到传入的变量中。 |
|
|
获取一个行向量，其中包含矩阵每列的平均值。第一种形式返回行向量。第二种形式将行向量保存到传入的变量中。 |
|
|
获取一个行向量，其中包含矩阵每列的中位数。第一种形式返回行向量。第二种形式将行向量保存到传入的变量中。 |
|
|
获取一个行向量，其中包含矩阵每列的总和。第一种形式返回行向量。第二种形式将行向量保存到传入的变量中。 |
|
|
获取矩阵元素的累积和。第一种形式返回累积和矩阵。第二种形式将累积和矩阵保存到传入的变量中。 |
矩阵类示例
##include <string> #include "Matrix.hpp" using namespace std; using namespace YMatrix; int main(int argc, char* argv[]) { // typedef for non-mathematical matrix with string as datatype. typedef Matrix<string,DenseMatrix,NonMathMatrix> MatrixNM; // typedef for normal mathematical matrix with double as datatype. typedef Matrix<double,DenseMatrix,MathMatrix> MatrixM; // typedef for mathematical matrix implementing expression template. typedef Matrix<double,DenseMatrix,MathMETMatrix> MatrixEM; // typedef for normal mathematical row vector with double as datatype. typedef RowVector<double,DenseMatrix,MathMatrix> RVectorM; // typedef for normal mathematical col vector with double as datatype. typedef ColVector<double,DenseMatrix,MathMatrix> CVectorM; int i, j; // construct a 3 x 4 matrix and fill it with 'Hello' for all elements. MatrixNM stringMatrix(3, 4, "Hello"); // resize matrix to 2,5 and fill new elements with 'World' stringMatrix.resize(2, 5, "World"); // print matrix. cout << "First way to print matrix\n"; for (i=0; i<stringMatrix.row.size(); ++i) { for (j=0; j<stringMatrix.col.size(); ++j) { cout << stringMatrix(i,j) << "\t"; } cout << "\n"; } cout << endl; vector<string> s(5, "New"); stringMatrix.row.insert(1, s.begin()); stringMatrix.col.insert(2, s.begin()); // another way to print matrix cout << "Second way to print matrix\n"; for (i=0; i<stringMatrix.row.size(); ++i) { copy (stringMatrix.row(i).begin(), stringMatrix.row(i).end(), ostream_iterator<string>(cout, "\t")); cout << "\n"; } cout << endl; stringMatrix.col.erase(0); // another way to print matrix cout << "Third way to print matrix\n"; cout << stringMatrix; // erase matrix. cout << "Deleting contents of matrix\n"; stringMatrix.clear(); // check if matrix is empty cout << (stringMatrix.empty() ? "Matrix is empty\n" : "Matrix is not empty\n"); // prepare two 3 x 3 arrays double a1[9] = {1,2,3, 4,5,6, 7,8,9}; double** a2 = new double*[3]; for (i=0; i<3; ++i) a2[i] = new double[3]; for (i=0; i<3; ++i) { for (j=0; j<3; ++j) { a2[i][j] = 9 - (i*3+j); } } // construct three 3 x 3 matrix with no initial value MatrixM mathMatrix1(3,3), mathMatrix2(3,3), mathMatrix3; // fill matrices with array mathMatrix1 = a1; mathMatrix2 = a2; cout << "First matrix:\n" << mathMatrix1; cout << "Second matrix:\n" << mathMatrix2; cout << "Element(1,2) of first matrix: " << mathMatrix1.at(1,2) << endl; // perform mathematical operations on matrix and print the results cout << "Some mathematical operations\n"; mathMatrix3 = mathMatrix1 + mathMatrix2; cout << mathMatrix3; mathMatrix3 = mathMatrix1 - mathMatrix2; cout << mathMatrix3; mathMatrix3 = mathMatrix1 * mathMatrix2; cout << mathMatrix3; mathMatrix3 = mathMatrix1 * 2; cout << mathMatrix3; mathMatrix3 = mathMatrix1 + mathMatrix2 * mathMatrix2; cout << mathMatrix3; mathMatrix3 += mathMatrix1; cout << mathMatrix3; mathMatrix3 -= mathMatrix1; cout << mathMatrix3; mathMatrix3 *= mathMatrix1; cout << mathMatrix3; mathMatrix3 *= mathMatrix1 + mathMatrix2; cout << mathMatrix3; mathMatrix3 *= 2; cout << mathMatrix3; // swap mathMatrix1 and mathMatrix2 mathMatrix1.swap(mathMatrix2); cout << "Swapping first and second matrix:\n"; cout << "First matrix:\n" << mathMatrix1; cout << "Second matrix:\n" << mathMatrix2; // swap row 0 and 1 of mathMatrix1 cout << "Swapping first and second row of first matrix:\n"; mathMatrix1.swaprows(0,1); cout << mathMatrix1; // swap column 0 and 2 of mathMatrix1 cout << "Swapping first and third column of first matrix:\n"; mathMatrix1.swapcols(0,2); cout << mathMatrix1; // print transpose of mathMatrix1 cout << "Transpose\n" << Transpose<MatrixM>()(mathMatrix1); // print diagonal of mathMatrix1 cout << "Diagonal\n" << Diagonal<MatrixM>()(mathMatrix1); // print covariance of mathMatrix1 cout << "Covariance\n" << Covariance<MatrixM>()(mathMatrix1); // print power 2 of mathMatrix1 cout << "Power 2\n" << Power<MatrixM>()(mathMatrix1, 2); // print mean of mathMatrix1 cout << "Mean\n" << Mean<MatrixM>()(mathMatrix1); // print median of mathMatrix1 cout << "Median\n" << Median<MatrixM>()(mathMatrix1); // print sum of mathMatrix1 cout << "Sum\n" << Sum<MatrixM>()(mathMatrix1); // print cumulative sum of mathMatrix1 cout << "CumulativeSum\n" << CumulativeSum<MatrixM>()(mathMatrix1); // construct three 3 x 3 matrix using the two arrays and mathMatrix1 MatrixEM mathEmatrix1(3, 3, a1); MatrixEM mathEmatrix2(3, 3, a2); MatrixEM mathEmatrix3(mathMatrix1); cout << "First matrix:\n" << mathEmatrix1; cout << "Second matrix:\n" << mathEmatrix2; // perform mathematical operations on matrix and print the results cout << "Some mathematical operations\n"; mathEmatrix3 = mathEmatrix1 + mathEmatrix2; cout << mathEmatrix3; mathEmatrix3 = mathEmatrix1 - mathEmatrix2; cout << mathEmatrix3; mathEmatrix3 = mathEmatrix1 * mathEmatrix2; cout << mathEmatrix3; mathEmatrix3 = mathEmatrix1 * 2; cout << mathEmatrix3; mathEmatrix3 = mathEmatrix1 + mathEmatrix2 * mathEmatrix2; cout << mathEmatrix3; mathEmatrix3 += mathEmatrix1; cout << mathEmatrix3; mathEmatrix3 -= mathEmatrix1; cout << mathEmatrix3; mathEmatrix3 *= mathEmatrix1; cout << mathEmatrix3; mathEmatrix3 *= 2; cout << mathEmatrix3; // invalid mathematical operation (compare with above) // mathEmatrix3 *= mathEmatrix1 + mathEmatrix2; cout << mathMatrix3; // demonstrate speed of expression template // (can't tell unless timing is done) mathEmatrix3 = mathEmatrix1 + mathEmatrix2 - mathEmatrix1 + mathEmatrix2; cout << mathMatrix3; // demonstrate interchangably of matrices mathEmatrix3 = mathMatrix1 + mathEmatrix2; cout << mathEmatrix3; // get transpose of mathEmatrix3 MatrixEM mathEmatrix3Transposed; Transpose<MatrixEM>()(mathEmatrix3, mathEmatrix3Transposed); cout << "Transposed:\n" << mathEmatrix3Transposed; RVectorM rv(5, 3.0); cout << rv; rv.push_back(8); cout << "push_back 8\n" << rv; rv.push_front(1); cout << "push_front 1\n" << rv; rv.pop_back(); cout << "pop_back\n" << rv; rv.pop_front(); cout << "pop_front\n" << rv; cout << endl; CVectorM cv(5, 2.6); cout << cv; cv.push_back(3.5); cout << "push_back 3.5\n" << cv; cv.push_front(8.5); cout << "push_front 8.5\n" << cv; cv.pop_back(); cout << "pop_back\n" << cv; cv.pop_front(); cout << "pop_front\n" << cv; cout << endl; }
实现细节
引言
此类矩阵的设计目标是通用性和可重用性。此设计的灵感来自 Andrei Alexandrescu4。由于矩阵是通用的，因此模板对于实现至关重要。矩阵必须能够存储不同类型的数据，而不仅仅是数值数据。因此，必须有一个 typename ValueT 来指示存储的数据类型。
有不同类型的矩阵：密集矩阵、稀疏矩阵、三角矩阵、对角矩阵等。每种类型的组织方式和数据存储方式都不同。因此，应该有一个策略来确定矩阵的存储方式和数据访问方式。因此，应为每种类型的矩阵创建单独的类。此类将处理数据的存储和访问。因此，矩阵类模板参数中的第二个元素是 template<typename> class MatrixTypeP。目前，MatrixTypeP 只能是 DenseMatrix。
设计矩阵类时必须考虑的第三件事是，为了使矩阵在科学编程中有用，它应该具有数学运算。但是，如果矩阵包含非数值类型，则不应提供数学运算。这个问题可以通过使用数学策略来解决。通过实现数学策略，包含非数值类型的矩阵无法执行数学运算，并且如果尝试执行，将导致编译器错误，而对于数值类型，数学运算将可用。最初，创建了两个类来实现此策略：MathMatrix 和 NonMathMatrix，矩阵类模板参数中的第三个元素是 template<typename,typename> 类 MathP。在矩阵类的初始版本中，MathP 从 MatrixTypeP 派生，然后矩阵类从任一数学类派生。后来更改为当前版本，矩阵类同时从 MatrixTypeP 和 MathP 派生，并且 MathP 不再从 MatrixTypeP 派生。此外，还创建了第三个数学类 MathMETMatrix，它实现了表达式模板以提高执行某些类型矩阵计算的性能。表达式模板的使用将在后面讨论。
因此，形成了矩阵类的骨架。矩阵类将接受三个模板参数，分别表示要存储的数据类型、存储和访问数据的策略，以及确定是否为矩阵提供数学运算的策略。唯一未考虑的是矩阵的维度。假设矩阵是二维的。多维矩阵未被考虑，因为这会增加实现的难度。
矩阵类的操作设计为类似于 STL 容器。矩阵类实现了各种类似 STL 的 typedef，如 value_type、size_type 等。还实现了用于遍历矩阵中元素的迭代器。实现了一些常见的 STL 容器函数，如 clear、empty、resize 和 swap。其他函数，如 assign、erase、insert、front、back、push_back、push_front、pop_back、pop_front、begin 和 end，在行和列迭代器中实现。
MatrixTypeP 的实现
任何 MatrixTypeP 矩阵都必须具有几个基本公共函数。必须有一个不接受参数的默认构造函数。应该有接受矩阵初始大小和初始值的构造函数。初始值可以是所有元素的默认值，也可以来自使用数组或指针实现的矩阵。由于存在原始指针作为成员变量，因此复制构造函数是必不可少的。应该有一个用于复制另一个类似 MatrixTypeP 矩阵的赋值运算符，以及一个用于复制任何矩阵类的赋值运算符。MatrixTypeP 矩阵还必须实现 resize、swap 和比较运算符。此外，还添加了一个 Update 函数，允许用户在数据更改时通知矩阵进行更新。DenseMatrix 不使用此 Update 函数，但对于其他 MatrixTypeP 来说，它可能是必需的。
还有其他基本公共函数，理论上可以将它们放在最终的矩阵类中，而不是放在 MatrixTypeP 矩阵中。它们是 clear、empty、各种成员访问函数、operator() 和 at 等。理论上可以将它们放在最终的矩阵类中，因为对于所有 MatrixTypeP 矩阵，它们的代码都无需更改。它们的实现基于所有 MatrixTypeP 矩阵中都会存在的通用成员变量和函数。然而，目前它们仍然放在 MatrixTypeP 矩阵中，因为它们提供了 MatrixTypeP 矩阵在实现其他函数时可能很有用的函数。
MatrixTypeP 矩阵中还应存在两个私有函数：ResizeIterators 和 UpdateIterators。ResizeIterators 用于调整四个私有成员变量：rowIteratorStart、rowIteratorFinish、colIteratorStart 和 colIteratorFinish 的大小。UpdateIterators 用于更新这四个成员变量存储的数据。最初，只有 UpdateIterators。四个私有成员变量的大小调整是在 UpdateIterators 内部完成的。之后，将大小调整与 UpdateIterators 分开，形成 ResizeIterators，以提高矩阵的性能。这四个私有成员变量是原始指针，用于保存指向每行和每列的开始和结束的迭代器的信息。通过存储这些迭代器，成员访问函数和各种迭代器函数可以被加速并对所有 MatrixTypeP 变得通用。最初，这四个私有成员变量是使用 STL vector 实现的。然后，为了提高性能，更改为使用原始指针。
最后，每个 MatrixTypeP 矩阵都必须实现一个行迭代器类和一个列迭代器类。行迭代器负责遍历矩阵中的一行，列迭代器负责遍历一列。这两个迭代器类是 MatrixTypeP 矩阵实现的最重要方面。我将尝试详细解释如何实现迭代器类。
由于迭代器类应类似于 STL 容器，因此它应具有各种类似 STL 的 typedef。因此，定义了 iterator_category、value_type、reference、pointer、size_type 和 difference_type。
迭代器应具有适当的构造函数。它必须有一个不接受参数的默认构造函数。它应该有一个接受迭代器（VectorTypeIterator）的构造函数，该迭代器遍历存储矩阵数据的类型的行或列。例如，在 DenseMatrix 中，矩阵数据由 value_type* 类型的原始指针存储。因此，VectorTypeIterator 的类型应为 value_type*。如果矩阵数据由 STL vector 存储，则 VectorTypeIterator 的类型为 std::vector::iterator。如果矩阵数据由 STL deque<deque> 存储，则 VectorTypeIterator 的类型为 std::deque::iterator。此 VectorTypeIterator 存储在成员变量中，供各种函数使用。除了 VectorTypeIterator 之外，构造函数可能还会接受任何有助于迭代器操作的变量。例如，DenseMatrix 的列迭代器中的构造函数接受矩阵的总列数。它在其 Subtract、Increment、Decrement 和 Advance 函数中使用此信息。
由于可以使用默认构造函数来构造迭代器，因此应该有一个接受 VectorTypeIterator 的赋值运算符。这是为了允许更新存储 VectorTypeIterator 的成员变量。还应提供其他函数来允许更新迭代器操作所需的任何变量。例如，有一个 SetCols 函数用于更新 DenseMatrix 的列迭代器中的总列数变量。
迭代器中应存在的函数类型将取决于迭代器的类型，即输入、输出、双向或随机。所有类型的迭代器都应实现 operator* 以允许解引用，以及 operator-> 以允许访问矩阵元素成员变量或函数。还必须为所有类型的迭代器实现友元比较函数。
对于随机访问迭代器，它应该具有前缀和后缀的递增和递减运算符。它还应该有一个 operator- 来确定两个类似迭代器之间的距离。此外，它应该实现 operator+=、operator-=、operator+ 和 operator- 以允许迭代器增加或减少某个距离。这些函数可以通过将任务委托给其他函数（如 Subtract、Increment、Decrement 和 Advance）来使其通用。随机访问迭代器还应实现 operator[] 以使其行为类似数组。
template<typename ValueT> class DenseMatrix
密集矩阵的实现经历了几个阶段。最初，使用 deque<deque<ValueT> > 来存储矩阵元素。然后，使用 deque<ValueT>。最后，使用 ValueT* 来提高矩阵的性能，因为使用 STL 容器会大大减慢矩阵计算速度。但是，在实现 DenseMatrix 类时，使用原始指针会稍微增加难度。实施过程中必须格外小心，以避免内存泄漏。此外，在发生异常时也可能发生内存泄漏。由于性能在科学编程中至关重要，因此在此情况下忽略了使用原始指针的这种缺点。
迭代器的实现与上面模板的实现略有不同。列迭代器使用上面的模板进行实现，但行迭代器使用原始指针 ValueT*。在初始版本中，使用一个遵循模板格式的类来创建行迭代器。后来发现，在矩阵计算中，重复使用行迭代器类中的函数实际上会减慢计算速度。由于数据存储允许将原始指针用于行迭代器，因此使用原始指针来提高矩阵计算速度。列迭代器不能这样做，因为遍历列需要了解矩阵的总列数，并且需要存储和使用此信息。因此，需要一个合适的类来用于列迭代器。由于使用行迭代器比使用列迭代器更快，因此在使用 DenseMatrix 存储策略时应使用行迭代器，除非绝对必要使用列迭代器以提高矩阵的效率。
使用 DenseMatrix 存储策略进行矩阵计算速度的近似分析表明，其速度远不如涉及原始指针或数组的手写循环。这在意料之中，因为底层的矩阵元素存储在多层代码之下，这会增加显著的开销。为了使此类矩阵对需要原始指针快速性能的科学程序员有吸引力，将三个成员变量和一个函数（通常应该是私有或保护的）设为公共。这违背了封装原则，但结果证明是正当的。三个变量是 rows_ 和 cols_，分别包含矩阵的总行数和列数，以及 matrix_，它是一个 ValueT* 类型的原始指针，包含矩阵元素。公开的函数是 ResizeIterators。追求速度但又希望方便的科学程序员可以利用矩阵类的各种函数来帮助维护矩阵，并在进行计算时，可以直接访问矩阵存储变量 matrix_，通过将其用作普通原始指针来执行快速查找、赋值或计算。公开 ResizeIterators 是为了让程序员在手动调整 matrix_ 大小或将其指向另一个矩阵时，可以更新内部迭代器，以保持矩阵结构。
MathP 的实现
必须在 MathP 类中实现的函数必须支持当数据为数值类型时矩阵的加法、减法、点乘和标量乘法。当数据为非数值类型时，所有这些函数都应不可用。有两种 MathP 类提供数学运算，一种则不提供。提供数学运算的两个类实现了 operator+=、operator+、operator-=、operator-、operator*= 和 operator*。
MathP 类的实现带来了一个独特的问题。MathP 类是最终矩阵类从中继承的基类。它与存储类是分开的，因此它不知道矩阵元素是如何存储的。由于存储类还负责提供访问器函数，因此 MathP 类实际上无法访问矩阵元素。但是，要执行矩阵计算，它必须能够访问矩阵元素。克服此问题的方法是将一个指向最终矩阵类的指针传递给 MathP 类，以便它可以访问矩阵元素。MathP 类必须维护一个指向其子矩阵的指针变量。因此，MathP 类的最终声明是 template<typename ValueT, typename MatrixT>，其中 MatrixT 是最终矩阵类的类型。由于 MathP 类在没有指向其子矩阵的指针的情况下无法存在，因此 MathP 类没有不接受参数的默认构造函数。只有一个需要指向子矩阵的指针的构造函数。
template<typename ValueT, typename MatrixT> class MathMatrix
MathMatrix 类相对简单。此类的主要功能是为最终矩阵类提供数学运算符，以便可以使用直观的数学运算符进行矩阵计算。
数学运算符的实现方式是，operator+ 和 operator- 分别调用 operator+= 和 operator-= 函数来执行其功能。此技术也用于实现标量乘法的 operator*= 和 operator*。但是，对于点乘，情况相反。operator*= 调用 operator* 函数来执行其功能。之所以反转，是因为它更有效。在矩阵点乘过程中，会创建一个临时矩阵来存储乘法结果。如果此代码位于 operator*= 而不是 operator* 中，则在使用 operator* 时，将创建三个临时矩阵，而不是通常的两个（在 operator+ 和 operator- 中）。这降低了矩阵点乘的效率。将代码放在 operator* 中，只创建两个临时矩阵，因此更有效。近似分析表明，这不会降低 operator*= 的性能，因此无需在 operator*= 中重复乘法代码。
MathMatrix 类实现进行各种矩阵计算前的尺寸检查。如果两个矩阵由于尺寸不匹配而无法相加、相减或点乘，它将返回一个异常。此检查不会显着降低矩阵计算速度，但如果确实需要性能，可以将其移除。
template<typename ValueT, typename MatrixT> class MathMETMatrix
MathMETMatrix 使用称为表达式模板的技术来提高某些类型矩阵加法和减法的性能。表达式模板由 Todd Veldhuizen5 引入。为了解释表达式模板，请考虑以下问题。
A、B、C、D、E、F 是矩阵。考虑代码 A = B + C + D + E + F;。为了执行矩阵计算，编译器将首先将 B 和 C 相加，然后将结果添加到 D，然后将结果添加到 E，然后将结果添加到 F，最后将最终结果放入 A。因此，在加法过程中会创建几个临时矩阵，这会大大降低性能。表达式模板的目标是通过创建等效于以下代码的代码来消除所有这些不必要的临时矩阵：
for (size_t i=0; i<rows; ++i) { for (size j=0; j<cols; ++j) { A(i,j) = B(i,j) + C(i,j) + D(i,j) +E(i,j) + F(i,j); } }
为了实现上述目标，使用了模板，并且需要几个支持类。我不会深入探讨表达式模板的理论。相反，我将解释各种支持模板类的创建。
第一个所需的模板类是操作类型（Operations type）的包装类。它命名为 MET。在 MathMETMatrix 类中，只有一种操作类型，即 METElementBinaryOp，它对矩阵中的元素执行二元运算。MET 的主要工作是包装 METElementBinaryOp，提供对 METElementBinaryOp 中存在的函数的访问。为此，它维护一个用于 METElementBinaryOp 模板类的变量。MET 是必需的，因为它有助于绑定不同的 METElementBinaryOp 模板类。没有 MET，表达式模板的整个概念就无法工作。
下一个模板类是操作类型。目前，只有 METElementBinaryOp。METElementBinaryOp 在其构造函数中接受两个模板变量，并将它们保存在指向这些变量数据的指针中。这两个模板变量可以是矩阵类或 MET 模板类。METElementBinaryOp 模板参数还接受一个运算符类型类。此运算符类型类决定对这两个模板变量执行的操作类型。
MathMETMatrix 中存在三种运算符类型：METAdd、METSubtract 和 METMultiply。METMultiply 尚未启用，因为尚未找到实现用于执行点乘的表达式模板的良好算法。METAdd 和 METSubtract 是提供静态函数 Evaluate 的结构。Evaluate 接受两个矩阵元素并返回它们相加或相减的结果。
为了实现表达式模板，重载了 operator+ 和 operator-，使其返回 MET 对象的实例。有一个 operator+ 和 operator- 用于执行 MathMETMatrix 派生矩阵与其他任何矩阵之间的运算。还有友元 operator+ 和 operator-，用于执行任何矩阵与任何 MET 对象之间，以及两个 MET 对象之间的运算。每个运算符重载都将返回一个 MET 对象，其中包含有关要执行的操作类型（目前仅为 METElementBinaryOp）、要执行的运算符类型（METAdd 用于 operator+，METSubtract 用于 operator-）以及用于执行运算符类型的变量（矩阵或表达式）的信息。
因此，当编译器遇到上述代码时：A = B + C + D + E + F;，它将生成以下代码：
A = MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >, Matrix>, Matrix>(& MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >, Matrix>, Matrix>(& MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >, Matrix>, Matrix>(& MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >, Matrix>(& MET< METElementBinaryOp< METAdd<value_type>, Matrix, Matrix> >(METElementBinaryOp< METAdd<value_type>, Matrix, Matrix>(&B, &C)), &D)), &E)), &F)); &C)), &D)), &E)), &F))
不要担心上面的代码。它的基本作用是生成一个 MET 对象，该对象封装了执行 B + C + D + E + F 所需的所有信息。现在需要做的就是实现结果的计算，并将其放入矩阵 A。
为了实现计算，在最终矩阵类中重载了赋值运算符（operator=）。赋值运算符将接受 MET 对象作为其参数。由于存在许多不同类型的 MET 对象，因此 operator= 必须是模板函数。结果代码如下：
template<typename ExprT> Self& operator=(MET<ExprT> expr) { size_t sizeRow = expr.rowsize; size_t sizeCol = expr.colsize; Self tempMatrix(sizeRow,sizeCol); for (size_t i=0; i<sizeRow; ++i) { for (size_t j=0; j<sizeCol; ++j) { tempMatrix(i,j) = expr(i,j); } } swap(tempMatrix); Update; return *this; }
代码相对容易理解。它的作用是从 MET 对象获取所需的总行数和列数，并创建一个临时矩阵来保存矩阵计算的结果。然后，它循环遍历每个单独的元素，并找出计算结果。以上面的例子为例，计算结果是 B(i,j) + C(i,j) + D(i,j) +E(i,j) + F(i,j) 的结果。然后将临时矩阵与当前矩阵（在示例中是矩阵 A）进行交换，以便计算结果现在分配给当前矩阵。
稍作离题，此函数比矩阵表达式模板3 慢，因为创建了临时矩阵。但是，此临时矩阵对于使操作更直观是必需的，即，在计算之前，矩阵 A 的大小无关紧要，因为它将被调整为 B + C + D + E + F 的矩阵计算所需的正确大小。矩阵表达式模板类的情况并非如此。在该类中，矩阵 A 需要与 B、C、D、E 和 F 大小相同，否则操作将失败。通过采用相同的立场，当前矩阵类可以实现比矩阵表达式模板更好的效率。只需移除临时矩阵，并将 tempMatrix(i,j) 替换为 (*this)(i,j)。但这样做是因为在这种情况下，操作的直观性比速度更重要。
现在，让我们回到表达式模板的讨论。从 operator= 代码可以看出，需要实现各种表达式模板支持类中的一些函数。需要实现两个函数 row.size 和 col.size 来找出接收结果的矩阵所需的总行数和列数。还需要实现一个接受矩阵元素位置并输出数学运算结果的运算符。
对于 MET 类，实现这三个函数很简单。这三个函数只是传递给 METElementBinaryOp 进行处理。
METElementBinaryOp 通过将其传递给其两个模板变量之一来实现 row.size 和 col.size。选择变量 left_。没有理由不能选择变量 right_。如果 left_ 是矩阵，它将具有 row.size 和 col.size 函数，因此将分别返回行数和列数。如果 left_ 是 MET 对象，它将将其传递给 METElementBinaryOp，后者将递归地将其传递给 left_，依此类推，直到到达矩阵。
METElementBinaryOp 通过在其 left_ 和 right_ 变量上应用运算符来实现运算符。如果 left_ 和 right_ 是矩阵，它们将返回位置 (i,j) 处的元素。如果它们是 MET 对象，它们将将其传递给 METElementBinaryOp，直到最终到达矩阵并返回位置 (i,j) 处的元素。然后将返回的元素发送到运算符类型的 Evaluate 函数，该函数将对其进行加法或减法运算并返回结果。因此，看起来无害的 expr(i,j) 实际上执行了所有必需的对单个元素的计算（例如，B(i,j) + C(i,j) + D(i,j) + E(i,j) + F(i,j)）。
到此为止，关于矩阵类的表达式模板实现就结束了。可以看出一个问题是，在进行矩阵计算之前，没有 proper 的方法来确保所有矩阵的大小都相同。因此，在执行计算之前，由用户确保矩阵大小相同。否则，可能会出现意外错误。另一个问题是，表达式模板可能不会加速所有矩阵的加法或减法。对于像 A = B + C; 这样的短加法，使用表达式模板实际上可能会减慢计算速度。因此，用户应尝试 MathMatrix 和 MathMETMatrix，以确定每种情况的最佳数学矩阵。两者都包含将一种形式转换为另一种形式的代码，并且两种类型的矩阵都可以用于同一次计算，因此在最终矩阵类中使用这两种类型没有损失。目前还没有实现使用表达式模板执行点乘的 proper 算法。因此，目前 MathMETMatrix 的点乘使用与 MathMatrix 相同的代码。
template<typename ValueT, typename MatrixT> NonMathMatrix
这是三个 MathP 类中最简单的。它包含最基本的功能，以便提供与 MathMatrix 和 MathMETMatrix 相同的接口。除此之外，它不包含任何数学运算符，因此从它派生的矩阵无法执行任何数学运算，并且在尝试执行时会产生编译器错误。
矩阵的实现
template<typename ValueT, template<typename> class MatrixTypeP, template<typename, typename> class MathP> class Matrix
我们已经到达了最重要的类，即 Matrix 类。它是用户与之交互的类。此类中存在的最终功能取决于它从哪个存储策略和哪个数学策略继承。
即使 Matrix 大部分函数继承自 MatrixTypeP 和 MathP，构造函数、析构函数、复制构造函数和赋值运算符也不能继承。因此，它需要实现这些。由于它是一个类似 STL 的容器，因此它必须定义各种类似 STL 的 typedef。这些 typedef 从 MatrixTypeP 获取，通过重新定义它们。MathP 策略没有任何构造函数或赋值运算符必须由 Matrix 镜像，因此 Matrix 只需定义 MatrixTypeP 中存在的构造函数、复制构造函数和赋值运算符，并将操作传递给它。因此，必须有一个不接受参数的默认构造函数。应该有接受矩阵初始大小和初始值的构造函数。初始值可以是所有元素的默认值，也可以来自使用数组或指针实现的矩阵。还应存在复制构造函数。应该有一个用于复制另一个类似矩阵的赋值运算符，以及一个用于复制任何矩阵类的赋值运算符。
此外，Matrix 实现了一个额外的复制构造函数，该构造函数在 MatrixTypeP 中不存在。此复制构造函数接受任何矩阵类。代码通过将要复制的矩阵传递给 MatrixTypeP 中的模板 operator= 来实现，而不是传递给 MatrixTypeP 中的类似复制构造函数。这样做的原因是，如果 MatrixTypeP 中有一个可以接受所有矩阵类的复制构造函数，它将阻止 Matrix 中的常规复制构造函数正常工作，因为传递到 MatrixTypeP 的所有矩阵都将由模板版本的复制构造函数而不是常规版本进行处理。这样做的原因很简单。Matrix 中常规复制构造函数的结构是 Matrix(const Self& m) : MathBase(this), StorageBase(m)。可以看出，StorageBase 接收 m 的常量引用。如果 StorageBase 中没有模板版本的复制构造函数，则 m 将由 StorageBase 中的常规复制构造函数处理。但是，如果 StorageBase 中存在模板版本，则 StorageBase 将使用模板版本来处理 m，因为 m 的类型是 Matrix，所以它不适合常规复制构造函数，后者需要 StorageBase 类型。因此，StorageBase 中不应存在模板版本的复制构造函数，否则所有复制都可能导致缓慢的操作或错误。然而，模板版本的复制构造函数将很有用，因为可能存在不同的 Matrix 类型，每种类型仅在它们的 MatrixTypeP 或 MathP 上有所不同，并且应该可以将一种类型复制到另一种类型。因此，模板版本的复制构造函数仅在 Matrix 类中创建，并将操作传递给 MatrixTypeP 中的模板 operator=，后者通常应通过所有矩阵可用的通用函数来执行矩阵复制。
Matrix 还实现了一个额外的赋值运算符（operator=），它接受一个 MET 对象。原因已在 MathMETMatrix 中讨论。
Matrix 中还实现了两个对所有矩阵都有用但对 MatrixTypeP 操作不是必需的函数：swaprows 和 swapcols。这两个函数仅使用迭代器进行操作，以便它们适用于所有类型的 MatrixTypeP 矩阵。
RowVector 和 ColVector 的实现
template<typename ValueT,template<typename> class MatrixTypeP,template<typename,typename> class MathP> class RowVector
template<typename ValueT,template<typename> class MatrixTypeP,template<typename,typename> class MathP> class ColVector
RowVector 和 ColVector 都继承自 Matrix 类。它们的存在几乎是不必要的，因为 Matrix 类可以通过将行数或列数分别设置为 1 来充当行向量或列向量。此外，STL 已经实现了三个有用且强大的 1 维容器。那么，为什么需要 RowVector 和 ColVector 呢？最初创建它们的想法是为了性能和易用性。
当 RowVector 和 ColVector 使用 DenseMatrix 类作为存储策略时，它们可以比 STL 容器更快。如前所述，DenseMatrix 向追求性能的用户公开了其存储机制。因此，用户可以通过公开的变量直接访问矩阵中的元素，而无需通过访问器函数。STL 容器不公开其存储机制，因此访问其元素的唯一方法是通过访问器函数。这会增加开销并降低性能。
对于在其生命周期内将成为向量的变量，RowVector 和 ColVector 比 Matrix 类更容易使用。对于这些变量，其行数或列数始终为 1。因此，Matrix 类在访问元素或获取迭代器方面提供的灵活性是不必要的。为了方便使用，RowVector 和 ColVector 实现了一些存在于 STL 容器中的函数： operator[]、at、push_back、resize、size、begin、end、rbegin 和 rend。RowVector 和 ColVector 还实现了只需要一个大小变量而不是 Matrix 类中的两个（一个用于行数，一个用于列数）的构造函数。
尽管 RowVector 和 ColVector 为访问和修改其元素和迭代器提供了简化的函数，但它们保留了 Matrix 类中可用的所有函数，以便与 Matrix 类兼容。它们也可以与 Matrix 类进行相互转换。但是，用户必须确保在将 Matrix 类复制或赋值给 RowVector 或 ColVector 类时，矩阵类必须分别是行向量或列向量（即，其行数或列数必须为 1）。否则，RowVector 或 ColVector 最终可能会变成一个矩阵。这不会导致使用此类时出现任何错误，用户甚至可能不会注意到它，但这在逻辑上是错误的，因为行向量或列向量不应存储矩阵。
作用于矩阵的函数
template<typename MatrixT> struct Transpose
Transpose 是一个查找矩阵转置的函数对象（转置矩阵是行和列交换的矩阵）。它有两个重载的运算符，第一个接受一个要转置的矩阵并返回转置后的矩阵，第二个接受一个要转置的矩阵和一个用于保存转置矩阵的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，或适用于任何合适的 STL 算法。第二个运算符更快，当需要找到矩阵的转置但不需要立即用于计算时，应该使用它。
template<typename MatrixT> struct Diagonal
Diagonal 是一个查找矩阵主对角线或将向量置于矩阵主对角线上的函数对象。它有两个重载的运算符，第一个接受一个矩阵并返回对角矩阵，第二个接受一个矩阵和一个用于保存对角矩阵的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，或适用于任何合适的 STL 算法。第二个运算符更快，当需要找到矩阵的对角线但不需要立即用于计算时，应该使用它。
template<typename MatrixT> struct Covariance
Covariance 是一个查找矩阵协方差的函数对象。它有两个重载的运算符，第一个接受一个矩阵并返回协方差矩阵，第二个接受一个矩阵和一个用于保存协方差矩阵的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，或适用于任何合适的 STL 算法。第二个运算符更快，当需要找到矩阵的协方差但不需要立即用于计算时，应该使用它。
template<typename MatrixT> struct Power
Power 是一个函子，用于计算将矩阵所有元素计算为指定次幂后的矩阵。它有两个重载的运算符，第一个接受一个矩阵并返回幂运算后的矩阵，第二个接受一个矩阵和一个用于存放幂运算后矩阵的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，也可以用于任何合适的 STL 算法。第二个运算符速度更快，当需要找到幂运算后的矩阵但不需要立即在计算中使用时，应使用它。
template<typename MatrixT> struct Mean
Mean 是一个函子，用于查找一个行向量，该行向量包含矩阵每列的均值。它有两个重载的运算符，第一个接受一个矩阵并返回行向量，第二个接受一个矩阵和一个用于存放行向量的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，也可以用于任何合适的 STL 算法。第二个运算符速度更快，当需要找到行向量但不需要立即在计算中使用时，应使用它。
template<typename MatrixT> struct Median
Median 是一个函子，用于查找一个行向量，该行向量包含矩阵每列的中位数。它有两个重载的运算符，第一个接受一个矩阵并返回行向量，第二个接受一个矩阵和一个用于存放行向量的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，也可以用于任何合适的 STL 算法。第二个运算符速度更快，当需要找到行向量但不需要立即在计算中使用时，应使用它。
template<typename MatrixT> struct Sum
Sum 是一个函子，用于查找一个行向量，该行向量包含矩阵每列的总和。它有两个重载的运算符，第一个接受一个矩阵并返回行向量，第二个接受一个矩阵和一个用于存放行向量的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，也可以用于任何合适的 STL 算法。第二个运算符速度更快，当需要找到行向量但不需要立即在计算中使用时，应使用它。
template<typename MatrixT> struct CumulativeSum
CumulativeSum 是一个函子，用于查找一个矩阵，该矩阵包含矩阵元素的累积和。它有两个重载的运算符，第一个接受一个矩阵并返回累积和矩阵，第二个接受一个矩阵和一个用于存放累积和矩阵的矩阵。第一个运算符适用于矩阵计算，可以直观地使用，也可以用于任何合适的 STL 算法。第二个运算符速度更快，当需要找到累积和矩阵但不需要立即在计算中使用时，应使用它。
结论
到此，关于 Matrix 类的讨论结束。Matrix 类在其实现中使用了策略和表达式模板等概念。它被设计成通用的、可重用的和易于扩展的。尽管其中一些概念可能没有被正确实现或以最佳方式实现，但它试图将各种 C++ 技术整合在一起。本文档用作如何使用 Matrix 类的示例，为所有各种类和函数提供文档，同时也为设计和实现这些类时遇到的各种阶段和思考提供文档。希望本文档能为所有希望使用、修改或扩展此类的人提供足够的信息。还希望对策略、迭代器和表达式模板的简要解释能激励那些对此感兴趣的人阅读更多关于这些技术的内容，并在自己的工作中加以使用。
非常欢迎对本类提出任何意见，任何新功能的请求都将被考虑，并在适当的时候尽快添加。
参考文献
- William H. Press, Numerical recipes in C++: the art of scientific computing (Cambridge University Press, 2002)。
- 矩阵模板库.
- 矩阵表达式模板.
- Andrei Alexandrescu, Modern C++ Design (Addison-Wesley, 2001)。
- Todd Veldhuizen, Expression Templates: C++ Report, Vol 7, No. 5, June 1995, pg 26-31。
更改日志
- 版本 1.1
- 移除了单位矩阵的实现，因为它没有得到很好的实现。
- 版本 1.2
- 添加了
RowVector和ColVector类。 - 修复了与反向迭代器相关的错误。
- 更改了
empty()函数。
- 添加了
- 版本 1.3
- 为 StoragePolicy 和 Expression Templates 添加了成员空间
row和col。 - 将
rowsize、colsize、rowbegin、rowend、colbegin和colend替换为row::size、col::size、row(i).begin、row(i).end、col(i).begin和col(i).end。 - 添加了
row::insert、row::erase、row::push_back、row::pop_back、col::insert、col::erase、col::push_back和col::pop_back。 - 为
RowVector和ColVector添加了pop_back。
- 为 StoragePolicy 和 Expression Templates 添加了成员空间
- 版本 1.4
- 将成员空间
row和col更改为Row和Col，以防止在声明RowVector或ColVector时引起row(i)无法工作的冲突。
- 将成员空间
- 版本 1.5 (2006 年 4 月 16 日)
- 修复了代码，使其不再依赖 STLport，并能在 Visual Studio 2005 和 GCC 3.3.2 上工作。
- 由于此修复，行和列迭代器的常量性可能无法真正起作用。
- 添加了
row::push_front、row::pop_front、col::push_front和col::pop_front。 - 为
RowVector和ColVector添加了push_front和pop_front。