部落格轉自：從Eigen向量化談記憶體對齊
Eigen是一個非常常用的矩陣運算庫，至少對於SLAM的研究者來說不可或缺。然而，有時候會由於Eigen向量化的記憶體對齊問題使程式執行異常。

事情起源:我的程式原本在NVIDIA TX2上跑的好好的，直到有一天，我打算把它放到伺服器上，看看傳說中的RTX 2080GPU能不能加速一把。結果悲劇發生了，編譯正常，但是一執行就立即double free。我很是吃驚，怎麼能一行程式碼都沒執行就崩了呢。但崩了就是崩了，一定是哪裡有bug，我用valgrind檢查記憶體問題，發現種種線索都指向g2o。g2o是一個SLAM後端優化庫，裡面封裝了大量SLAM相關的優化演算法，內部使用了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位元組對齊。