博客轉自：從Eigen向量化談內存對齊
Eigen是一個非常常用的矩陣運算庫，至少對於SLAM的研究者來說不可或缺。然而，有時候會由於Eigen向量化的內存對齊問題使程序運行異常。
事情起源:我的程序原本在NVIDIA TX2上跑的好好的，直到有一天，我打算把它放到服務器上，看看傳說中的RTX 2080GPU能不能加速一把。結果悲劇發生了，編譯正常，但是一運行就立即double free。我很是吃驚，怎么能一行代碼都沒執行就崩了呢。但崩了就是崩了，一定是哪里有bug，我用valgrind檢查內存問題，發現種種線索都指向g2o。g2o是一個SLAM后端優化庫，里面封裝了大量SLAM相關的優化算法，內部使用了Eigen進行矩陣運算。陰差陽錯之間，我發現關閉-march=native這個編譯選項后就能正常運行，而這個編譯選項其實是告訴編譯器當前的處理器支持哪些SIMD指令集，Eigen中又恰好使用了SSE、AVX等指令集進行向量化加速。此時，機智的我發現Eigen文檔中有一章叫做Alignment issues，里面提到了某些情況下Eigen對象可能沒有內存對齊，從而導致程序崩潰。現在，證據到齊，基本可以確定我遇到的真實問題了：編譯安裝g2o時，默認沒有使用-march=native，因此里面的Eigen代碼沒有使用向量化加速，所以它們並沒有內存對齊。而在我的程序中，啟用了向量化加速，所有的Eigen對象都是內存對齊的。兩個程序鏈接起來之后，g2o中未對齊的Eigen對象一旦傳遞到我的代碼中，向量化運算的指令就會觸發異常。解決方案很簡單，要么都用-march=native，要么都不用。
這件事就這么過去了，但我不能輕易放過它，畢竟花費了那么多時間找bug。后來我又做了一些深入的探究，這篇文章就來談談向量化和內存對齊里面的門道。
什么是向量化運算
向量化運算就是用SSE、AVX等SIMD（Single Instruction Multiple Data）指令集，實現一條指令對多個操作數的運算，從而提高代碼的吞吐量，實現加速效果。SSE是一個系列，包括從最初的SSE到最新的SSE4.2，支持同時操作16 bytes的數據，即4個float或者2個double。AVX也是一個系列，它是SSE的升級版，支持同時操作32 bytes的數據，即8個float或者4個double。
但向量化運算是有前提的，那就是內存對齊。SSE的操作數，必須16 bytes對齊，而AVX的操作數，必須32 bytes對齊。也就是說，如果我們有4個float數，必須把它們放在連續的且首地址為16的倍數的內存空間中，才能調用SSE的指令進行運算。
A Simple Example
為了給沒接觸過向量化編程的同學一些直觀的感受，我寫了一個簡單的示例程序：
#include <immintrin.h>
#include <iostream>
int main() {
double input1[4] = {1, 1, 1, 1};
double input2[4] = {1, 2, 3, 4};
double result[4];
std::cout << "address of input1: " << input1 << std::endl;
std::cout << "address of input2: " << input2 << std::endl;
__m256d a = _mm256_load_pd(input1);
__m256d b = _mm256_load_pd(input2);
__m256d c = _mm256_add_pd(a, b);
_mm256_store_pd(result, c);
std::cout << result[0] << " " << result[1] << " " << result[2] << " " << result[3] << std::endl;
return 0;
}
這段代碼使用AVX中的向量化加法指令，同時計算4對double的和。這4對數保存在input1和input2中。 _mm256_load_pd指令用來加載操作數，_mm256_add_pd指令進行向量化運算，最后， _mm256_store_pd指令讀取運算結果到result中。可惜的是，程序運行到第一個_mm256_load_pd處就崩潰了。崩潰的原因正是因為輸入變量沒有內存對齊。我特意打印出了兩個輸入變量的地址，結果如下
address of input1: 0x7ffeef431ef0
address of input2: 0x7ffeef431f10
上一節提到了AVX要求32字節對齊，我們可以把這兩個輸入變量的地址除以32，看是否能夠整除。結果發現0x7ffeef431ef0和 0x7ffeef431f10都不能整除。當然，其實直接看倒數第二位是否是偶數即可，是偶數就可以被32整除，是奇數則不能被32整除。
如何讓輸入變量內存對齊呢？我們知道，對於局部變量來說，它們的內存地址是在編譯期確定的，也就是由編譯器決定。所以我們只需要告訴編譯器，給input1和input2申請空間時請讓首地址32字節對齊，這需要通過預編譯指令來實現。不同編譯器的預編譯指令是不一樣的，比如gcc的語法為__attribute__((aligned(32)))，MSVC的語法為 __declspec(align(32)) 。以gcc語法為例，做少量修改，就可以得到正確的代碼
#include <immintrin.h>
#include <iostream>
int main() {
__attribute__ ((aligned (32))) double input1[4] = {1, 1, 1, 1};
__attribute__ ((aligned (32))) double input2[4] = {1, 2, 3, 4};
__attribute__ ((aligned (32))) double result[4];
std::cout << "address of input1: " << input1 << std::endl;
std::cout << "address of input2: " << input2 << std::endl;
__m256d a = _mm256_load_pd(input1);
__m256d b = _mm256_load_pd(input2);
__m256d c = _mm256_add_pd(a, b);
_mm256_store_pd(result, c);
std::cout << result[0] << " " << result[1] << " " << result[2] << " " << result[3] << std::endl;
return 0;
}
輸出結果為
address of input1: 0x7ffc5ca2e640
address of input2: 0x7ffc5ca2e660
2 3 4 5
可以看到，這次的兩個地址都是32的倍數，而且最終的運算結果也完全正確。
雖然上面的代碼正確實現了向量化運算，但實現方式未免過於粗糙。每個變量聲明前面都加上一長串預編譯指令看起來就不舒服。我們嘗試重構一下這段代碼。
重構
首先，最容易想到的是，把內存對齊的double數組聲明成一種自定義數據類型，如下所示
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
aligned_double4 input1 = {1, 1, 1, 1};
aligned_double4 input2 = {1, 2, 3, 4};
aligned_double4 result;
這樣看起來清爽多了。更進一步，如果4個double是一種經常使用的數據類型的話，我們就可以把它封裝為一個Vector4d類，這樣，用戶就完全看不到內存對齊的具體實現了，像下面這樣。
#include <immintrin.h>
#include <iostream>
class Vector4d {
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
public:
Vector4d() {
}
Vector4d(double d1, double d2, double d3, double d4) {
data[0] = d1;
data[1] = d2;
data[2] = d3;
data[3] = d4;
}
aligned_double4 data;
};
Vector4d operator+ (const Vector4d& v1, const Vector4d& v2) {
__m256d data1 = _mm256_load_pd(v1.data);
__m256d data2 = _mm256_load_pd(v2.data);
__m256d data3 = _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 << "(" << v.data[0] << ", " << v.data[1] << ", " << v.data[2] << ", " << v.data[3] << ")";
return o;
}
int main() {
Vector4d input1 = {1, 1, 1, 1};
Vector4d input2 = {1, 2, 3, 4};
Vector4d result = input1 + input2;
std::cout << result << std::endl;
return 0;
}
這段代碼實現了Vector4d類，並把向量化運算放在了operator+中，主函數變得非常簡單。但不要高興得太早，這個Vector4d其實有着嚴重的漏洞，如果我們動態創建對象，程序仍然會崩潰，比如這段代碼
int main() {
Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
std::cout << "address of input1: " << input1->data << std::endl;
std::cout << "address of input2: " << input2->data << std::endl;
Vector4d result = *input1 + *input2;
std::cout << result << std::endl;
delete input1;
delete input2;
return 0;
}
崩潰前的輸出為
address of input1: 0x1ceae70
address of input2: 0x1ceaea0
很詭異吧，似乎剛才我們設置的內存對齊都失效了，這兩個輸入變量的內存首地址又不是32的倍數了。
Heap vs Stack
問題的根源在於不同的對象創建方式。直接聲明的對象是存儲在棧上的，其內存地址由編譯器在編譯時確定，因此預編譯指令會生效。但用new動態創建的對象則存儲在堆中，其地址在運行時確定。C++的運行時庫並不會關心預編譯指令聲明的對齊方式，我們需要更強有力的手段來確保內存對齊。
C++提供的new關鍵字是個好東西，它避免了C語言中丑陋的malloc操作，但同時也隱藏了實現細節。如果我們翻看C++官方文檔，可以發現new Vector4d實際上做了兩件事情，第一步申請sizeof(Vector4d)大小的空間，第二步調用Vector4d的構造函數。要想實現內存對齊，我們必須修改第一步申請空間的方式才行。好在第一步其實調用了operator new這個函數，我們只需要重寫這個函數，就可以實現自定義的內存申請，下面是添加了該函數后的Vector4d類。
class Vector4d {
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
public:
Vector4d() {
}
Vector4d(double d1, double d2, double d3, double d4) {
data[0] = d1;
data[1] = d2;
data[2] = d3;
data[3] = d4;
}
void* operator new (std::size_t count) {
void* original = ::operator new(count + 32);
void* aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~size_t(32 - 1)) + 32);
*(reinterpret_cast<void**>(aligned) - 1) = original;
return aligned;
}
void operator delete (void* ptr) {
::operator delete(*(reinterpret_cast<void**>(ptr) - 1));
}
aligned_double4 data;
};
operator new的實現還是有些技巧的，我們來詳細解釋一下。 首先，根據C++標准的規定，operator new的參數count是要開辟的空間的大小。 為了保證一定可以得到count大小且32字節對齊的內存空間，我們把實際申請的內存空間擴大到count + 32。可以想象，在這count + 32字節空間中， 一定存在首地址為32的倍數的連續count字節的空間。 所以，第二行代碼，我們通過對申請到的原始地址original做一些位運算，先找到比original小且是32的倍數的地址，然后加上32，就得到了我們想要的對齊后的地址，記作aligned。 接下來，第三行代碼很關鍵，它把原始地址的值保存在了aligned地址的前一個位置中，之所以要這樣做，是因為我們還需要自定義釋放內存的函數operator delete。畢竟aligned地址並非真實申請到的地址，所以在該地址上調用默認的delete是會出錯的。可以看到，我們在代碼中也定義了一個operator delete，傳入的參數正是前面operator new返回的對齊的地址。這時候，保存在aligned前一個位置的原始地址就非常有用了，我們只需要把它取出來，然后用標准的delete釋放該內存即可。
為了方便大家理解這段代碼，有幾個細節需要特地強調一下。::operator new中的::代表全局命名空間，因此可以調用到標准的operator new。第三行需要先把aligned強制轉換為void**類型，這是因為我們希望在aligned的前一個位置保存一個void*類型的地址，既然保存的元素是地址，那么該位置對應的地址就是地址的地址，也就是void**。
這是一個不大不小的trick，C++的很多內存管理方面的處理經常會有這樣的操作。但不知道細心的你是否發現了這里的一個問題：reinterpret_cast<void**>(aligned) - 1這個地址是否一定在我們申請的空間中呢？換句話說， 它是否一定大於original呢？ 之所以存在這個質疑，是因為這里的-1其實是對指針減一。要知道，在64位計算機中，指針的長度是8字節，所以這里得到的地址其實是reinterpret_cast<size_t>(aligned) - 8。看出這里的區別了吧，對指針減1相當於對地址的值減8。所以仔細想想，如果original到aligned的距離小於8字節的話，這段代碼就會對申請的空間以外的內存賦值，可怕吧。
其實沒什么可怕的，為什么我敢這樣講，因為Eigen就是這樣實現的。這樣做依賴於現代編譯器的一個共識：所有的內存分配都默認16字節對齊。這個事實可以解釋很多問題，首先，永遠不用擔心original到aligned的距離會不會小於8了，它會穩定在16，這足夠保存一個指針。其次，為什么我們用AVX指令集舉例，而不是SSE？因為SSE要求16字節對齊，而現代編譯器已經默認16字節對齊了，那這篇文章就沒辦法展開了。 最后，為什么我的代碼在NVIDIA TX2上運行正常而在服務器上掛掉了？因為TX2中是ARM處理器，里面的向量化指令集NEON也只要求16字節對齊。
仍然存在BUG
如果你以為到這里就圓滿結束了，那可是大錯特錯。還有個天坑沒展示給大家，下面的代碼中，我的自定義類Point包含了一個Vector4d的成員，這時候，噩夢又出現了。
class Point {
public:
Point(Vector4d position) : position(position) {
}
Vector4d position;
};
int main() {
Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
Point* point1 = new Point{*input1};
Point* point2 = new Point{*input2};
std::cout << "address of point1: " << point1->position.data << std::endl;
std::cout << "address of point2: " << point2->position.data << std::endl;
Vector4d result = point1->position + point2->position;
std::cout << result << std::endl;
delete input1;
delete input2;
delete point1;
delete point2;
return 0;
}
輸出的地址又不再是32的倍數了，程序戛然而止。我們分析一下為什么會這樣。在主函數中，new Point動態創建了一個Point對象。前面提到過，這個過程分為兩步，第一步申請Point對象所需的空間，即sizeof(Point)大小的空間，第二步調用Point的構造函數。我們寄希望於第一步申請到的空間恰好讓內部的position對象對齊，這是不現實的。因為整個過程中並不會調用Vector4d的operator new，調用的只有Point的operator new，而這個函數我們並沒有重寫。
可惜的是，此處並沒有足夠優雅的解決方案，唯一的方案是在Point類中也添加自定義operator new，這就需要用戶的協助，類庫的作者已經無能為力了。 不過類庫的作者能做的，是盡量讓用戶更方便地添加operator new，比如封裝為一個宏定義，用戶只需要在Point類中添加一句宏即可。最后，完整的代碼如下。
#include <immintrin.h>
#include <iostream>
#define ALIGNED_OPERATOR_NEW \
void* operator new (std::size_t count) { \
void* original = ::operator new(count + 32); \
void* aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~size_t(32 - 1)) + 32); \
*(reinterpret_cast<void**>(aligned) - 1) = original; \
return aligned;\
} \
void operator delete (void* ptr) { \
::operator delete(*(reinterpret_cast<void**>(ptr) - 1)); \
}
class Vector4d {
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
public:
Vector4d() {
}
Vector4d(double d1, double d2, double d3, double d4) {
data[0] = d1;
data[1] = d2;
data[2] = d3;
data[3] = d4;
}
ALIGNED_OPERATOR_NEW
aligned_double4 data;
};
Vector4d operator+ (const Vector4d& v1, const Vector4d& v2) {
__m256d data1 = _mm256_load_pd(v1.data);
__m256d data2 = _mm256_load_pd(v2.data);
__m256d data3 = _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 << "(" << v.data[0] << ", " << v.data[1] << ", " << v.data[2] << ", " << v.data[3] << ")";
return o;
}
class Point {
public:
Point(Vector4d position) : position(position) {
}
ALIGNED_OPERATOR_NEW
Vector4d position;
};
int main() {
Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
Point* point1 = new Point{*input1};
Point* point2 = new Point{*input2};
std::cout << "address of point1: " << point1->position.data << std::endl;
std::cout << "address of point2: " << point2->position.data << std::endl;
Vector4d result = point1->position + point2->position;
std::cout << result << std::endl;
delete input1;
delete input2;
delete point1;
delete point2;
return 0;
}
這段代碼中，宏定義ALIGNED_OPERATOR_NEW包含了operator new和operator delete，它們對所有需要內存對齊的類都適用。因此，無論是需要內存對齊的類，還是包含了這些類的類，都需要添加這個宏。
再談Eigen
在Eigen官方文檔中有這么[一頁內容
有沒有覺得似曾相識？Eigen對該問題的解決方案與我們不謀而合。這當然不是巧合，事實上，本文的靈感正是來源於Eigen。但Eigen只告訴了我們應該怎么做，沒有詳細講解其原理。本文則從問題的提出，到具體的解決方案，一一剖析，希望可以給大家一些更深的理解。
總結
最后做一個簡短的總結。對於基本數據類型和自定義類型，我們需要用預編譯指令來保證棧內存的對齊，用重寫operator new的方式保證堆內存對齊。對於嵌套的自定義類型，申請棧內存時會自動保證其內部數據類型的對齊，而申請堆內存時仍然需要重寫operator new。
有一種特殊情況本文並未提到，如果使用std::vector<Vector4d>，需要傳入自定義內存申請器，即std::vector<Vector4d, AlignedAllocator>，其中AlignedAllocator是我們自定義的內存申請器。這是因為，std::vector中使用了動態申請的空間保存數據，因此默認的operator new是無法讓其內存對齊的。在無法重寫std::vector類的operator new的情況下，標准庫提供了自定義內存申請器的機制，讓用戶可以以自己的方式申請內存。具體做法本文就不再展開了，理解了前面的內容，這個問題應該很容易解決。
本文用到的所有示例代碼已上傳GitHub：jingedawang/AlignmentExample。
參考資料
Eigen Memory Issues ethz-asl/eigen_catkin wiki
Explanation of the assertion on unaligned arrays Eigen Doc
在C/C++代碼中使用SSE等指令集的指令(4)SSE指令集Intrinsic函數使用 gengshenghong
alignas cppreference
Data Alignment, Part 1 Noel Llopis
Data Alignment, Part 2: Objects on The Heap and The Stack Noel Llopis
GCC中的aligned和packed屬性 Shengbin
new expression cppreference
Why are all arrays aligned to 16 bytes on my implementation? stackoverflow
Why is dynamically allocated memory always 16 bytes aligned? stackoverflow