最大公因數範例

這是 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() 做一些調整。

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

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

另外一個 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 值會介於 0numQPUs - 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 網站上的說明文件以及其原始碼。