cmake怎么编译eigenc++_从Eigen向量化谈内存对齐

作者：好人森森_195 | 来源：互联网 | 2023-09-03 03:56

缘起Eigen是一个非常常用的矩阵运算库，至少对于SLAM的研究者来说不可或缺。然而，向来乖巧的Eigen近来却频频闹脾气，把我的程序折腾

缘起

Eigen是一个非常常用的矩阵运算库&＃xff0c;至少对于SLAM的研究者来说不可或缺。然而&＃xff0c;向来乖巧的Eigen近来却频频闹脾气&＃xff0c;把我的程序折腾得死去活来&＃xff0c;我却是丈二和尚摸不着头脑。

简单说说我经历的灵异事件。我的程序原本在NVIDIA TX2上跑的好好的&＃xff0c;直到有一天&＃xff0c;我打算把它放到服务器上&＃xff0c;看看传说中的RTX 2080GPU能不能加速一把。结果悲剧发生了&＃xff0c;编译正常&＃xff0c;但是一运行就立即double free。我很是吃惊&＃xff0c;怎么能一行代码都没执行就崩了呢。但崩了就是崩了&＃xff0c;一定是哪里有bug&＃xff0c;我用valgrind检查内存问题&＃xff0c;发现种种线索都指向g2o。g2o是一个SLAM后端优化库&＃xff0c;里面封装了大量SLAM相关的优化算法&＃xff0c;内部使用了Eigen进行矩阵运算。阴差阳错之间&＃xff0c;我发现关闭-march&＃61;native这个编译选项后就能正常运行&＃xff0c;而这个编译选项其实是告诉编译器当前的处理器支持哪些SIMD指令集&＃xff0c;Eigen中又恰好使用了SSE、AVX等指令集进行向量化加速。此时&＃xff0c;机智的我发现Eigen文档中有一章叫做Alignment issues&＃xff0c;里面提到了某些情况下Eigen对象可能没有内存对齐&＃xff0c;从而导致程序崩溃。现在&＃xff0c;证据到齐&＃xff0c;基本可以确定我遇到的真实问题了&＃xff1a;编译安装g2o时&＃xff0c;默认没有使用-march&＃61;native&＃xff0c;因此里面的Eigen代码没有使用向量化加速&＃xff0c;所以它们并没有内存对齐。而在我的程序中&＃xff0c;启用了向量化加速&＃xff0c;所有的Eigen对象都是内存对齐的。两个程序链接起来之后&＃xff0c;g2o中未对齐的Eigen对象一旦传递到我的代码中&＃xff0c;向量化运算的指令就会触发异常。解决方案很简单&＃xff0c;要么都用-march&＃61;native&＃xff0c;要么都不用。

这件事就这么过去了&＃xff0c;但我不能轻易放过它&＃xff0c;毕竟花费了那么多时间找bug。后来我又做了一些深入的探究&＃xff0c;这篇文章就来谈谈向量化和内存对齐里面的门道。

什么是向量化运算&＃xff1f;

向量化运算就是用SSE、AVX等SIMD&＃xff08;Single Instruction Multiple Data&＃xff09;指令集&＃xff0c;实现一条指令对多个操作数的运算&＃xff0c;从而提高代码的吞吐量&＃xff0c;实现加速效果。SSE是一个系列&＃xff0c;包括从最初的SSE到最新的SSE4.2&＃xff0c;支持同时操作16 bytes的数据&＃xff0c;即4个float或者2个double。AVX也是一个系列&＃xff0c;它是SSE的升级版&＃xff0c;支持同时操作32 bytes的数据&＃xff0c;即8个float或者4个double。

但向量化运算是有前提的&＃xff0c;那就是内存对齐。SSE的操作数&＃xff0c;必须16 bytes对齐&＃xff0c;而AVX的操作数&＃xff0c;必须32 bytes对齐。也就是说&＃xff0c;如果我们有4个float数&＃xff0c;必须把它们放在连续的且首地址为16的倍数的内存空间中&＃xff0c;才能调用SSE的指令进行运算。

A Simple Example

为了给没接触过向量化编程的同学一些直观的感受&＃xff0c;我写了一个简单的示例程序&＃xff1a;

#include #include int main() {double input1[4] &＃61; {1, 1, 1, 1};double input2[4] &＃61; {1, 2, 3, 4};double result[4];std::cout <<"address of input1: " <}

这段代码使用AVX中的向量化加法指令&＃xff0c;同时计算4对double的和。这4对数保存在input1和input2中。 _mm256_load_pd指令用来加载操作数&＃xff0c;_mm256_add_pd指令进行向量化运算&＃xff0c;最后&＃xff0c; _mm256_store_pd指令读取运算结果到result中。可惜的是&＃xff0c;程序运行到第一个_mm256_load_pd处就崩溃了。崩溃的原因正是因为输入变量没有内存对齐。我特意打印出了两个输入变量的地址&＃xff0c;结果如下

address of input1: 0x7ffeef431ef0 address of input2: 0x7ffeef431f10

上一节提到了AVX要求32字节对齐&＃xff0c;我们可以把这两个输入变量的地址除以32&＃xff0c;看是否能够整除。结果发现0x7ffeef431ef0和 0x7ffeef431f10都不能整除。当然&＃xff0c;其实直接看倒数第二位是否是偶数即可&＃xff0c;是偶数就可以被32整除&＃xff0c;是奇数则不能被32整除。

如何让输入变量内存对齐呢&＃xff1f;我们知道&＃xff0c;对于局部变量来说&＃xff0c;它们的内存地址是在编译期确定的&＃xff0c;也就是由编译器决定。所以我们只需要告诉编译器&＃xff0c;给input1和input2申请空间时请让首地址32字节对齐&＃xff0c;这需要通过预编译指令来实现。不同编译器的预编译指令是不一样的&＃xff0c;比如gcc的语法为__attribute__((aligned(32)))&＃xff0c;MSVC的语法为 __declspec(align(32)) 。以gcc语法为例&＃xff0c;做少量修改&＃xff0c;就可以得到正确的代码

#include #include int main() {__attribute__ ((aligned (32))) double input1[4] &＃61; {1, 1, 1, 1};__attribute__ ((aligned (32))) double input2[4] &＃61; {1, 2, 3, 4};__attribute__ ((aligned (32))) double result[4];std::cout <<"address of input1: " <}

输出结果为

address of input1: 0x7ffc5ca2e640 address of input2: 0x7ffc5ca2e660 2 3 4 5

可以看到&＃xff0c;这次的两个地址都是32的倍数&＃xff0c;而且最终的运算结果也完全正确。

虽然上面的代码正确实现了向量化运算&＃xff0c;但实现方式未免过于粗糙。每个变量声明前面都加上一长串预编译指令看起来就不舒服。我们尝试重构一下这段代码。

重构

首先&＃xff0c;最容易想到的是&＃xff0c;把内存对齐的double数组声明成一种自定义数据类型&＃xff0c;如下所示

using aligned_double4 &＃61; __attribute__ ((aligned (32))) double[4];aligned_double4 input1 &＃61; {1, 1, 1, 1};aligned_double4 input2 &＃61; {1, 2, 3, 4};aligned_double4 result;

这样看起来清爽多了。更进一步&＃xff0c;如果4个double是一种经常使用的数据类型的话&＃xff0c;我们就可以把它封装为一个Vector4d类&＃xff0c;这样&＃xff0c;用户就完全看不到内存对齐的具体实现了&＃xff0c;像下面这样。

#include #include class Vector4d {using aligned_double4 &＃61; __attribute__ ((aligned (32))) double[4]; public:Vector4d() {}Vector4d(double d1, double d2, double d3, double d4) {data[0] &＃61; d1;data[1] &＃61; d2;data[2] &＃61; d3;data[3] &＃61; d4;}aligned_double4 data; };Vector4d operator&＃43; (const Vector4d& v1, const Vector4d& v2) {__m256d data1 &＃61; _mm256_load_pd(v1.data);__m256d data2 &＃61; _mm256_load_pd(v2.data);__m256d data3 &＃61; _mm256_add_pd(data1, data2);Vector4d result;_mm256_store_pd(result.data, data3);return result; }std::ostream& operator<<(std::ostream& o, const Vector4d& v) {o <<"(" <}int main() {Vector4d input1 &＃61; {1, 1, 1, 1};Vector4d input2 &＃61; {1, 2, 3, 4};Vector4d result &＃61; input1 &＃43; input2;std::cout <}

这段代码实现了Vector4d类&＃xff0c;并把向量化运算放在了operator&＃43;中&＃xff0c;主函数变得非常简单。

但不要高兴得太早&＃xff0c;这个Vector4d其实有着严重的漏洞&＃xff0c;如果我们动态创建对象&＃xff0c;程序仍然会崩溃&＃xff0c;比如这段代码

int main() {Vector4d* input1 &＃61; new Vector4d{1, 1, 1, 1};Vector4d* input2 &＃61; new Vector4d{1, 2, 3, 4};std::cout <<"address of input1: " <data <data <}

崩溃前的输出为

address of input1: 0x1ceae70 address of input2: 0x1ceaea0

很诡异吧&＃xff0c;似乎刚才我们设置的内存对齐都失效了&＃xff0c;这两个输入变量的内存首地址又不是32的倍数了。

Heap vs Stack

问题的根源在于不同的对象创建方式。直接声明的对象是存储在栈上的&＃xff0c;其内存地址由编译器在编译时确定&＃xff0c;因此预编译指令会生效。但用new动态创建的对象则存储在堆中&＃xff0c;其地址在运行时确定。C&＃43;&＃43;的运行时库并不会关心预编译指令声明的对齐方式&＃xff0c;我们需要更强有力的手段来确保内存对齐。

C&＃43;&＃43;提供的new关键字是个好东西&＃xff0c;它避免了C语言中丑陋的malloc操作&＃xff0c;但同时也隐藏了实现细节。如果我们翻看C&＃43;&＃43;官方文档&＃xff0c;可以发现new Vector4d实际上做了两件事情&＃xff0c;第一步申请sizeof(Vector4d)大小的空间&＃xff0c;第二步调用Vector4d的构造函数。要想实现内存对齐&＃xff0c;我们必须修改第一步申请空间的方式才行。好在第一步其实调用了operator new这个函数&＃xff0c;我们只需要重写这个函数&＃xff0c;就可以实现自定义的内存申请&＃xff0c;下面是添加了该函数后的Vector4d类。

class Vector4d {using aligned_double4 &＃61; __attribute__ ((aligned (32))) double[4]; public:Vector4d() {}Vector4d(double d1, double d2, double d3, double d4) {data[0] &＃61; d1;data[1] &＃61; d2;data[2] &＃61; d3;data[3] &＃61; d4;}void* operator new (std::size_t count) {void* original &＃61; ::operator new(count &＃43; 32);void* aligned &＃61; reinterpret_cast((reinterpret_cast(original) & ~size_t(32 - 1)) &＃43; 32);*(reinterpret_cast(aligned) - 1) &＃61; original;return aligned;}void operator delete (void* ptr) {::operator delete(*(reinterpret_cast(ptr) - 1));}aligned_double4 data; };

operator new的实现还是有些技巧的&＃xff0c;我们来详细解释一下。首先&＃xff0c;根据C&＃43;&＃43;标准的规定&＃xff0c;operator new的参数count是要开辟的空间的大小。为了保证一定可以得到count大小且32字节对齐的内存空间&＃xff0c;我们把实际申请的内存空间扩大到count &＃43; 32。可以想象&＃xff0c;在这count &＃43; 32字节空间中&＃xff0c; 一定存在首地址为32的倍数的连续count字节的空间。所以&＃xff0c;第二行代码&＃xff0c;我们通过对申请到的原始地址original做一些位运算&＃xff0c;先找到比original小且是32的倍数的地址&＃xff0c;然后加上32&＃xff0c;就得到了我们想要的对齐后的地址&＃xff0c;记作aligned。接下来&＃xff0c;第三行代码很关键&＃xff0c;它把原始地址的值保存在了aligned地址的前一个位置中&＃xff0c;之所以要这样做&＃xff0c;是因为我们还需要自定义释放内存的函数operator delete。毕竟aligned地址并非真实申请到的地址&＃xff0c;所以在该地址上调用默认的delete 是会出错的。可以看到&＃xff0c;我们在代码中也定义了一个operator delete&＃xff0c;传入的参数正是前面operator new返回的对齐的地址。这时候&＃xff0c;保存在aligned前一个位置的原始地址就非常有用了&＃xff0c;我们只需要把它取出来&＃xff0c;然后用标准的delete释放该内存即可。

为了方便大家理解这段代码&＃xff0c;有几个细节需要特地强调一下。::operator new中的::代表全局命名空间&＃xff0c;因此可以调用到标准的operator new。第三行需要先把aligned强制转换为void**类型&＃xff0c;这是因为我们希望在aligned的前一个位置保存一个void*类型的地址&＃xff0c;既然保存的元素是地址&＃xff0c;那么该位置对应的地址就是地址的地址&＃xff0c;也就是void**。

这是一个不大不小的trick&＃xff0c;C&＃43;&＃43;的很多内存管理方面的处理经常会有这样的操作。但不知道细心的你是否发现了这里的一个问题&＃xff1a;reinterpret_cast(aligned) - 1这个地址是否一定在我们申请的空间中呢&＃xff1f;换句话说&＃xff0c; 它是否一定大于original呢&＃xff1f; 之所以存在这个质疑&＃xff0c;是因为这里的-1其实是对指针减一。要知道&＃xff0c;在64位计算机中&＃xff0c;指针的长度是8字节&＃xff0c;所以这里得到的地址其实是reinterpret_cast(aligned) - 8。看出这里的区别了吧&＃xff0c;对指针减1相当于对地址的值减8。所以仔细想想&＃xff0c;如果original到aligned的距离小于8字节的话&＃xff0c;这段代码就会对申请的空间以外的内存赋值&＃xff0c;可怕吧。

其实没什么可怕的&＃xff0c;为什么我敢这样讲&＃xff0c;因为Eigen就是这样实现的。这样做依赖于现代编译器的一个共识&＃xff1a;所有的内存分配都默认16字节对齐。这个事实可以解释很多问题&＃xff0c;首先&＃xff0c;永远不用担心original到aligned的距离会不会小于8了&＃xff0c;它会稳定在16&＃xff0c;这足够保存一个指针。其次&＃xff0c;为什么我们用AVX指令集举例&＃xff0c;而不是SSE&＃xff1f;因为SSE要求16字节对齐&＃xff0c;而现代编译器已经默认16字节对齐了&＃xff0c;那这篇文章就没办法展开了。最后&＃xff0c;为什么我的代码在NVIDIA TX2上运行正常而在服务器上挂掉了&＃xff1f;因为TX2中是ARM处理器&＃xff0c;里面的向量化指令集NEON也只要求16字节对齐。

还有坑&＃xff1f;

如果你以为到这里就圆满结束了&＃xff0c;那可是大错特错。还有个天坑没展示给大家&＃xff0c;下面的代码中&＃xff0c;我的自定义类Point包含了一个Vector4d的成员&＃xff0c;这时候&＃xff0c;噩梦又出现了。

class Point { public:Point(Vector4d position) : position(position) {}Vector4d position; };int main() {Vector4d* input1 &＃61; new Vector4d{1, 1, 1, 1};Vector4d* input2 &＃61; new Vector4d{1, 2, 3, 4};Point* point1 &＃61; new Point{*input1};Point* point2 &＃61; new Point{*input2};std::cout <<"address of point1: " <position.data <position.data <position &＃43; point2->position;std::cout <}

输出的地址又不再是32的倍数了&＃xff0c;程序戛然而止。我们分析一下为什么会这样。在主函数中&＃xff0c;new Point动态创建了一个Point对象。前面提到过&＃xff0c;这个过程分为两步&＃xff0c;第一步申请Point对象所需的空间&＃xff0c;即sizeof(Point)大小的空间&＃xff0c;第二步调用Point的构造函数。我们寄希望于第一步申请到的空间恰好让内部的position对象对齐&＃xff0c;这是不现实的。因为整个过程中并不会调用Vector4d的operator new&＃xff0c;调用的只有Point的operator new&＃xff0c;而这个函数我们并没有重写。

可惜的是&＃xff0c;此处并没有足够优雅的解决方案&＃xff0c;唯一的方案是在Point类中也添加自定义operator new&＃xff0c;这就需要用户的协助&＃xff0c;类库的作者已经无能为力了。不过类库的作者能做的&＃xff0c;是尽量让用户更方便地添加operator new&＃xff0c;比如封装为一个宏定义&＃xff0c;用户只需要在Point类中添加一句宏即可。最后&＃xff0c;完整的代码如下。

#include #include #define ALIGNED_OPERATOR_NEW void* operator new (std::size_t count) { void* original &＃61; ::operator new(count &＃43; 32); void* aligned &＃61; reinterpret_cast((reinterpret_cast(original) & ~size_t(32 - 1)) &＃43; 32); *(reinterpret_cast(aligned) - 1) &＃61; original; return aligned;} void operator delete (void* ptr) { ::operator delete(*(reinterpret_cast(ptr) - 1)); }class Vector4d {using aligned_double4 &＃61; __attribute__ ((aligned (32))) double[4]; public:Vector4d() {}Vector4d(double d1, double d2, double d3, double d4) {data[0] &＃61; d1;data[1] &＃61; d2;data[2] &＃61; d3;data[3] &＃61; d4;}ALIGNED_OPERATOR_NEWaligned_double4 data; };Vector4d operator&＃43; (const Vector4d& v1, const Vector4d& v2) {__m256d data1 &＃61; _mm256_load_pd(v1.data);__m256d data2 &＃61; _mm256_load_pd(v2.data);__m256d data3 &＃61; _mm256_add_pd(data1, data2);Vector4d result;_mm256_store_pd(result.data, data3);return result; }std::ostream& operator<<(std::ostream& o, const Vector4d& v) {o <<"(" <}class Point { public:Point(Vector4d position) : position(position) {}ALIGNED_OPERATOR_NEWVector4d position; };int main() {Vector4d* input1 &＃61; new Vector4d{1, 1, 1, 1};Vector4d* input2 &＃61; new Vector4d{1, 2, 3, 4};Point* point1 &＃61; new Point{*input1};Point* point2 &＃61; new Point{*input2};std::cout <<"address of point1: " <position.data <position.data <position &＃43; point2->position;std::cout <}

这段代码中&＃xff0c;宏定义ALIGNED_OPERATOR_NEW 包含了operator new和operator delete&＃xff0c;它们对所有需要内存对齐的类都适用。因此&＃xff0c;无论是需要内存对齐的类&＃xff0c;还是包含了这些类的类&＃xff0c;都需要添加这个宏。

再谈Eigen

在Eigen官方文档中有这么一页内容

有没有觉得似曾相识&＃xff1f;Eigen对该问题的解决方案与我们不谋而合。这当然不是巧合&＃xff0c;事实上&＃xff0c;本文的灵感正是来源于Eigen。但Eigen只告诉了我们应该怎么做&＃xff0c;没有详细讲解其原理。本文则从问题的提出&＃xff0c;到具体的解决方案&＃xff0c;一一剖析&＃xff0c;希望可以给大家一些更深的理解。

总结

最后做一个简短的总结。对于基本数据类型和自定义类型&＃xff0c;我们需要用预编译指令来保证栈内存的对齐&＃xff0c;用重写operator new的方式保证堆内存对齐。对于嵌套的自定义类型&＃xff0c;申请栈内存时会自动保证其内部数据类型的对齐&＃xff0c;而申请堆内存时仍然需要重写operator new。

有一种特殊情况本文并未提到&＃xff0c;如果使用std::vector &＃xff0c;需要传入自定义内存申请器&＃xff0c;即std::vector&＃xff0c;其中AlignedAllocator是我们自定义的内存申请器。这是因为&＃xff0c;std::vector中使用了动态申请的空间保存数据&＃xff0c;因此默认的operator new是无法让其内存对齐的。在无法重写std::vector类的operator new的情况下&＃xff0c;标准库提供了自定义内存申请器的机制&＃xff0c;让用户可以以自己的方式申请内存。具体做法本文就不再展开了&＃xff0c;理解了前面的内容&＃xff0c;这个问题应该很容易解决。

本文用到的所有示例代码已上传GitHub&＃xff1a;jingedawang/AlignmentExample。