樹莓派 Raspberry Pi GPU 加速運算教學：QPULib 函式庫

最大公因數範例

這是 Euclidean 演算法的 C 語言實作版本，可以找出兩個整數的最大公因數（greatest common divisor）：

void gcd(int* p, int* q, int* r) {
  int a = *p;
  int b = *q;
  while (a != b) {
    if (a > b)
      a = a-b;
    else
      b = b-a;
  }
  *r = a;
}

這裡我們以指標（pointer）的方式來傳遞參數，如果是一般的 C 語言程式是可以不需要這麼麻煩，直接以傳值呼叫即可，這裡這樣寫是為了後續我們要把這個函數改為平行化的向量版本。

若要讓這段城市可以在 GPU 上面執行，我們要使用 QPULib 的語法改寫這一個函數，以下是改寫後的版本：

#include <QPULib.h>
void gcd(Ptr<Int> p, Ptr<Int> q, Ptr<Int> r) {

  // 取讀輸入向量，長度為 16
  Int a = *p;
  Int b = *q;

  // 平行檢查 a 向量與 b 向量，
  // 若 a 與 b 中有任何元素不相等，
  // 則繼續執行 While 迴圈
  While (any(a != b))

    // 個別比較 a 向量與 b 向量中的元素，
    // 針對 a 大於 b 的那些元素進行運算
    Where (a > b)
      a = a - b;
    End  // Where 條件判斷式結尾

    // 個別比較 a 向量與 b 向量中的元素，
    // 針對 a 小於 b 的那些元素進行運算
    Where (a < b)
      b = b - a;
    End  // Where 條件判斷式結尾

  End  // While 迴圈結尾

  // 傳回計算結果
  *r = a;
}

QPULib 是以 C++ 物件導向的方式建立自己的程式語言語法，所以表面上看起來是 C++ 的程式，事實上在使用起來有很大的差異，用這種語法所撰寫的程式碼是放在 QPU 中執行的，所有的運算都是以長度為 16 的向量為基本單位，以平行化的方式在執行。

• Int 是一個 QPULib 所定義的類別，代表長度為 16 的 32 位元整數向量。
• Ptr<Int> 則是代表 Int 向量的記憶體位址向量。
• *p 代表記憶體中以 p 開頭的 Int 向量（跟 C 語言的陣列相同，第一個向量元素位址為 p）。
• a != b 代表兩個向量個別元素的比較，結果是長度為 16 的布林向量。
• any(a != b) 是判斷 a != b 這個布林向量中，是否有任何值為 true。
• Where (a > b) 與 End 兩個是一組的條件判斷式，向量中只有符合條件判斷的對應元素會執行情中的程式碼。

以上是 GPU 部分的核心函數，以下是 CPU 部分的 C/C++ 程式碼：

int main() {
  // 將 gcd 函數編譯成 QPU 核心函數
  auto k = compile(gcd);

  // 配置 CPU 與 QPU 共享的陣列
  SharedArray<int> a(16), b(16), r(16);

  // 設定亂數種子
  srand(0);

  // 以亂數產生測試用資料
  for (int i = 0; i < 16; i++) {
    a[i] = 100 + rand()%100;
    b[i] = 100 + rand()%100;
  }

  // 設定要使用的 QPU 數量
  k.setNumQPUs(1);

  // 使用 QPU 進行平行運算
  k(&a, &b, &r);

  // 輸出結果
  for (int i = 0; i < 16; i++)
    printf("gcd(%i, %i) = %i\n", a[i], b[i], r[i]);

  return 0;
}

這個範例就是前面我們在安裝 QPULib 函式庫測試時所用的 GCD 程式，編譯與執行方式請參考上面的說明。

回圈展開

回圈展開（loop unrolling）是程式碼最佳化常用的技巧，可以有效降低分支指令（branch instructions）的執行次數，加速整體程式的速度。

QPU 的分支指令會需要 3 個延遲間隙（delay slots），相當於 12 個時脈週期，目前 QPULib 在這個間隙中並沒有做其他用途（也就是浪費掉了），所以若要加速程式的執行，必須儘可能減少條件判斷出現的頻率。

QPULib 允許程式設計者以一般 C 語言的 for 迴圈，對特定的程式碼做回圈展開，例如：

#include <QPULib.h>
void gcd(Ptr<Int> p, Ptr<Int> q, Ptr<Int> r) {
  Int a = *p;
  Int b = *q;
  While (any(a != b))
    for (int i = 0; i < 32; i++) {  // 迴圈展開
      Where (a > b)
        a = a - b;
      End
      Where (a < b)
        b = b - a;
      End
    }  // for 迴圈展開結尾
  End
  *r = a;
}

這裡加入的 for 迴圈並不是真正的迴圈，而是做迴圈展開，也就是讓編譯器產生多個重複的 QPU 指令，改善程式執行速度。

迴圈範例

這是一個旋轉座標的 C 語言函數：

void rot3D(int n, float cosTheta, float sinTheta, float* x, float* y) {
  for (int i = 0; i < n; i++) {
    float xOld = x[i];
    float yOld = y[i];
    x[i] = xOld * cosTheta - yOld * sinTheta;
    y[i] = yOld * cosTheta + xOld * sinTheta;
  }
}

若以 QPULib 改寫後，會變成這樣：

void rot3D(Int n, Float cosTheta, Float sinTheta, Ptr<Float> x, Ptr<Float> y) {
  For (Int i = 0, i < n, i = i + 16)
    Float xOld = x[i];
    Float yOld = y[i];
    x[i] = xOld * cosTheta - yOld * sinTheta;
    y[i] = yOld * cosTheta + xOld * sinTheta;
  End
}

Non-Blocking 記憶體存取

上面的程式碼只是最初淺的 GPU 版本，這段程式碼在執行時會因為記憶體搬移資料，造成程式產生多餘的等待時間，若要改善效能可以將計算與記憶體存取的時間區間重疊，讓程式一邊計算、一邊存取記憶體，以下是修改過後的程式碼：

void rot3D(Int n, Float cosTheta, Float sinTheta, Ptr<Float> x, Ptr<Float> y) {
  // index() 函數可產生 0, 1, 2, ..., 15 的向量
  Ptr<Float> p = x + index();
  Ptr<Float> q = y + index();

  // 預先載入資料
  gather(p); gather(q);

  Float xOld, yOld;
  For (Int i = 0, i < n, i = i+16)

    // 預先載入下一次迭代要用的資料
    gather(p+16); gather(q+16);

    // 取得當次迭代要用的資料
    receive(xOld); receive(yOld);

    // 將計算結果儲存回記憶體中
    store(xOld * cosTheta - yOld * sinTheta, p);
    store(yOld * cosTheta + xOld * sinTheta, q);

    p = p+16; q = q+16;
  End
}

QPULib 提供了 non-blocking 的記憶體載入方式：

• gather(p) 可以預先載入 p 中的每一個值。
• receive(x) 可以將 gather 載入的值存入 x 中，程式在這裡會等待所有的直都載入至 x 中之後，才會繼續執行（也就是會 block 住）。

gather(p) 的語法與 x = *p 不同，它是針對每一個記憶體位址操作，所以要使用 index() 做一些調整。

記憶體的資料會在 gather 與 receive 之間進行載入，所以這段時間我們可以加入其他的計算，讓計算與資料載入同時進行。

QPU 中有 8 個 FIFO 緩衝區可以儲存 gather 載入的資料，而資料在載入 FIFO 之後，再以 receive 取出，也就是說同時間可以有 8 組 gather 與 receive 一起運作。

另外一個 non-blocking 操作是使用 store 將資料寫回記憶體，store(x, p) 可以將 x 向量寫入 p 這個記憶體位置，與 *p = x 不同的是，store(x, p) 是 non-blocking 的操作，不過若遇到下一個 store 時，就要等待前一個 store 執行完成才能接續執行。

多 QPU 平行運算

QPULib 可以透過 setNumQPUs(n) 函數來設定 QPU 數量，同時使用多顆 QPU 來做平行計算，而在使用多顆 QPU 時，程式設計者必須要能辨識其所使用的 GPU，這樣才能將計算分配到每一顆 QPU 中。

QPULib 的 me() 函數可以傳回 QPU 的 id（型態為 Int），而 numQPUs() 則會傳回可以使用的 QPU 總數，QPU 的 id 值會介於 0 到 numQPUs - 1 之間。

以下是使用多顆 QPU 平行計算的版本：

void rot3D(Int n, Float cosTheta, Float sinTheta, Ptr<Float> x, Ptr<Float> y) {
  Int inc = numQPUs() << 4;  // numQPUs() * 16
  Ptr<Float> p = x + index() + (me() << 4);
  Ptr<Float> q = y + index() + (me() << 4);
  gather(p); gather(q);

  Float xOld, yOld;
  For (Int i = 0, i < n, i = i + inc)
    gather(p+inc); gather(q+inc);
    receive(xOld); receive(yOld);
    store(xOld * cosTheta - yOld * sinTheta, p);
    store(yOld * cosTheta + xOld * sinTheta, q);
    p = p + inc; q = q + inc;
  End
}

以上就是簡單的 QPULib 介紹與教學，更詳細的語法與範例請參考 QPULib GibHub 網站上的說明文件以及其原始碼。