英文原文(What happens inside Eigen, on a simple example)
考虑以下示例程序:
#include<Eigen/Core>
int main()
{
int size = 50;
// VectorXf is a vector of floats, with dynamic size.
Eigen::VectorXf u(size), v(size), w(size);
u = v + w;
}
本页的目标是了解 Eigen 如何编译,假设启用了 SSE2
矢量化(GCC 选项 -msse2
)。
为什么讨论这个问题
也许你认为上面的示例程序很简单,编译它应该不涉及任何非常有趣的事情。在开始之前,让我们解释一下编译它时的一些不普通的部分——也就是生成优化代码——以便我们在这里解释Eigen的复杂性真正有用。
看一下这行代码:
u = v + w; // (*)
编译它的第一件重要的事情是数组应该只遍历一次,就像:
for(int i = 0; i < size; i++) u[i] = v[i] + w[i];
问题是,如果我们创建一个简单的 C++ 库,其中 VectorXf
类有一个返回 VectorXf
的运算符+
,那么代码行 (*)
将相当于:
VectorXf tmp = v + w;
VectorXf u = tmp;
显然,这里引入临时tmp
是没有用的。它对性能有非常糟糕的影响,首先是因为 tmp
的创建需要在此上下文中进行动态内存分配,其次是因为现在有两个 for
循环:
for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i];
for(int i = 0; i < size; i++) u[i] = tmp[i];
如果我们需要进行两次数组遍历而不是一次,这对性能来说非常糟糕,因为这意味着我们需要做很多冗余的内存访问。
编译上述程序的第二个重要方面是正确使用SSE2
指令。请注意,Eigen还支持AltiVec
,我们在这里讨论的所有内容也适用于AltiVec
。
像AltiVec
一样,SSE2
是一组指令,允许一次处理128位数据。由于float
占32位,这意味着SSE2
指令可以一次处理4个float
。这意味着,如果正确使用,它们可以使我们的计算速度提高4倍。
然而,在上面的程序中,我们选择了size=50
,因此我们的向量由50个float
组成,而50不是4的倍数。这意味着我们不能指望使用SSE2
指令来执行所有的计算。我们应该着眼于次优解,即使用SSE2
指令处理前48个系数,因为48是小于50的最大4的倍数,然后单独处理第49个和第50个系数。类似这样:
for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i);
for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];
因此,让我们逐行查看示例程序,然后跟随 Eigen 进行编译。
构建向量
我们分析一下第一行代码:
Eigen::VectorXf u(size), v(size), w(size);
首先,VectorXf
是 typedef
定义类型:
typedef Matrix<float, Dynamic, 1> VectorXf;
类模板Matrix
在src/Core/util/ForwardDeclarations.h
中声明,有6个模板参数,但最后3个参数由前3个自动确定,所以暂时不解释它们。在这里,Matrix<float,Dynamic,1>
表示一个浮点数矩阵,具有1列和动态的行数。
Matrix
类继承了基类 MatrixBase
。可以说 MatrixBase
统一了矩阵/向量和所有表达式类型,下面将详细介绍。
当只执行 Eigen::VectorXf u(size);
时,调用的构造函数是 Matrix::Matrix(int)
,位于 src/Core/Matrix.h 中。除了一些断言之外,它所做的只是构造 m_storage
成员,其类型为 DenseStorage<float, Dynamic, Dynamic, 1>
。
你可能会想,把存储放在一个单独的类中是不是设计过度了?原因是Matrix类模板涵盖了各种矩阵和向量:固定大小和动态大小。这两种情况下的存储方法不同。对于固定大小的矩阵,矩阵系数被存储为一个普通的成员数组。对于动态大小的矩阵,系数将被存储为指向动态分配数组的指针。因此,需要将存储从Matrix类中抽象出来。这就是DenseStorage
的作用。
让我们看看在 src/Core/DenseStorage.h 中的这个构造函数。这里有许多DenseStorages
的部分模板特化,单独处理维度为Dynamic
或在编译时固定的情况。我们要看的部分特化是:
template<typename T, int Cols_> class DenseStorage<T, Dynamic, Dynamic, Cols_>
这里,调用的构造函数是 DenseStorage::DenseStorage(int size, int rows, int columns)
,其中 size=50
、rows=50
、columns=1
。
这个构造函数如下:
inline DenseStorage(int size, int rows, int) : m_data(internal::aligned_new<T>(size)), m_rows(rows) {
}
这里,m_data
成员是矩阵系数的实际数组,它是动态分配的。与其调用new[]
或malloc()
,Eigen 定义了自己的internal::aligned_new
,它位于 src/Core/util/Memory.h中。它的作用是,如果启用了向量化,那么它会使用平台特定的调用来分配一个128位对齐的数组,因为对于使用SSE2
和AltiVec
进行向量化非常有用。如果未启用向量化,则相当于标准的new[]
。
构造函数还将m_rows
成员设置为size
。请注意,这里没有m_columns
成员:实际上,在DenseStorage
的这个部分特化中,列数是编译时确定的,因为Cols_
模板参数与Dynamic
不同。在这种情况下,Cols_
为1
,这意味着我们的向量只是一个只有1
列的矩阵。因此,没有必要将列数存储为运行时变量。
当调用VectorXf::data()
获取系数数组的指针时,它会返回DenseStorage::data()
,后者返回m_data
成员。
当调用VectorXf::size()
获取向量的大小时,实际上是调用了基类MatrixBase
中的一个方法。它确定向量是一个列向量,因为ColsAtCompileTime==1
(这来自于typedef VectorXf
中的模板参数)。它推断出大小就是行数,所以返回VectorXf::rows()
,后者返回DenseStorage::rows()
,后者返回构造函数设置的m_rows
成员。