這裡介紹如何在樹莓派上使用 QPULib 這套 QPU（GPU） 平行運算函式庫，加速各種運算的執行，解決樹莓派 CPU 運算速度不足的問題。
QPULib 是一套可以運用樹莓派的 QPU 進行平行運算的 C++ 函式庫，它包含特殊的程式語言語法以及編譯器，可在執行時產生適用於 QPU 執行的程式，讓 QPU 負責較為大量的運算，大幅增加運算的速度，對於有實時（real-time）需求的應用應該很有幫助。以下是 QPULib 的使用教學與範例程式碼。
這裡我所使用的測試環境是 Raspberry Pi 3 Model B，作業系統與核心環境如下：
lsb_release -a
No LSB modules are available. Distributor ID: Raspbian Description: Raspbian GNU/Linux 8.0 (jessie) Release: 8.0 Codename: jessie
uname -a
Linux raspberrypi 4.4.43-v7+ #948 SMP Sun Jan 15 22:20:07 GMT 2017 armv7l GNU/Linux
GPU 運算基礎觀念
QPU（Quad Processing Unit）是 Broadcom 所研發的一種向量處理器，它可以一次處理 16 個 32 位元整數或浮點數的運算。
舉例來說，如果有兩個長度為 16 的數值向量：
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
與
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
QPU 的 integer-add 指令可以將這兩個向量相加，得到一個相加後的結果向量：
30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60
這個新的向量長度也是 16，其中每個元素就是由前面兩個向量中對應元素相加所得到的。
每一個長度為 16 的向量包含 4 個 quads，一個 QPU 每一個時脈週期（clock cycle）可以處理一個 quad（這也是它稱為 QPU 的原因），所以計算一個長度為 16 的向量需要 4 個時脈週期。
樹莓派中總共含有 12 個 QPU，時脈速度（clock rate）為 250MHz，也就是說每秒可以進行 250M / 4 * 12 = 750M 個長度為 16 的向量運算，等同於每秒 750M * 16 = 12G 個浮點運算，而某些特殊的 QPU 指令還可以一次執行兩個運算，所以樹莓派的 QPU 最高可達到 24GFLOPS 的運算能力。
QPU 屬於樹莓派圖形管線（pipeline）的一部分，如果您的應用會牽涉到繪圖，可以使用 OpenGL ES，而如果只是要加速非繪圖的程式運算速度，那麼就適合使用 QPULib。
安裝 QPULib
使用 git 下載 QPULib 專案的原始碼：
# 安裝 git
sudo apt-get install git
# 下載 QPULib 原始碼
git clone https://github.com/mn416/QPULib
編譯 GCD 這個範例程式，編譯時要加上 QPU=1 參數，讓程式可以使用樹莓派實體的 GPU（QPU）：
cd QPULib/Tests
make QPU=1 GCD
使用 sudo 執行 GCD 測試：
sudo ./GCD
輸出為：
gcd(183, 186) = 3 gcd(177, 115) = 1 gcd(193, 135) = 1 gcd(186, 192) = 6 gcd(149, 121) = 1 gcd(162, 127) = 1 gcd(190, 159) = 1 gcd(163, 126) = 1 gcd(140, 126) = 14 gcd(172, 136) = 4 gcd(111, 168) = 3 gcd(167, 129) = 1 gcd(182, 130) = 26 gcd(162, 123) = 3 gcd(167, 135) = 1 gcd(129, 102) = 3
若程式可以正常執行，就表示沒問題。
Hello World
在撰寫使用 GPU 運算的程式時，程式碼主要可分為兩大部份：
- CPU 部份
- 在 CPU 上面執行的程式，使用一般 C/C++ 語言的語法來撰寫（等同於一般的 C/C++ 程式），主要處理資料的輸入與輸出，還有一些次要的少量運算。
- GPU 部份
- 在 GPU 上面以平行化方式所執行的程式，使用 QPULib 以 C++ 類別為基礎所建立的特殊語法撰寫而成，負責最核心的平行運算，通常這部份的程式碼會包裝成一個獨立的函數，跟 CPU 的程式碼區隔開來，而這種在 GPU 上面執行的函數就稱為核心函數（kernel function）。
這是 QPULib 所提供的 Hello.cpp 範例程式碼，示範最簡單的 CPU 與 GPU 資料傳遞與平行運算的使用方式：
#include <QPULib.h>
// 定義在 QPU 上執行的核心函數
void hello(Ptr<Int> p) {
*p = 1;
}
int main() {
// 建立核心函數
auto k = compile(hello);
// 配置 CPU 與 QPU 共用的記憶體陣列
SharedArray<int> array(16);
for (int i = 0; i < 16; i++) {
array[i] = 100;
}
// 使用 QPU 運算
k(&array);
// 輸出結果
for (int i = 0; i < 16; i++) {
printf("%i: %i\n", i, array[i]);
}
return 0;
}
這裡我們一開始定義一個 hello 核心函數，這個函數就是要放在 GPU 執行的部份，在程式執行時，先使用 compile 函數將這個 hello 編譯成可以被 GPU 執行的程式，接著用 SharedArray 配置一個可以同時被 CPU 與 GPU 存取的陣列，用於 GPU 運算上的資料傳遞，然後呼叫編建立好的核心函數進行平行運算，最後輸出結果。
這個 Hello 程式的 makefile 已經事先寫好了，所以我們可以直接執行 make 編譯程式，這裡同樣要加上 QPU=1 讓程式使用 GPU 運算：
make QPU=1 Hello
使用 sudo 執行：
sudo ./Hello
輸出為：
0: 1 1: 1 2: 1 3: 1 4: 1 5: 1 6: 1 7: 1 8: 1 9: 1 10: 1 11: 1 12: 1 13: 1 14: 1 15: 1
這裡我們放在 GPU 中的運算內容很簡單，只是單純把每個數值設定為 1 而已，所以輸出是 16 個 1。
以上就是一個最簡單的樹莓派 GPU hello world 程式，接下來我們要繼續介紹 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 網站上的說明文件以及其原始碼。