CUDA 學習筆記 (二) 【Chapter4 CUDA並行程式設計】
Capter4 CUDA並行程式設計
前面我們看到將一個標準C函式放到GPU裝置上執行是很容易的。只需要在函式定義前面加上 __globle__ 修飾符，並通過一種特殊的尖括號語法來呼叫它，就可以在GPU上執行這個函式。然而，前面的示例只調用了一個和函式，並且該函式在GPU上以序列方式執行。本章將學習如何啟動一個並行執行的裝置函式。
4.1 向量求和運算
我們先通過一個簡單的示例說明執行緒的概念，以及如何使用 CUDA C來實現執行緒。假設有兩組資料，我們需要將這兩組資料中對應的元素兩兩相加，並將結果儲存在第三陣列中。
4.1.1 基於CPU的向量求和
首先看看如何通過傳統的C程式碼來實現這個求和運算：
#include "../common/book.h" #define N 10 void add( int *a, int *b, int *c ) { int tid = 0; // 這是第0個CPU,因此索引從0開始 while (tid < N) { c[tid] = a[tid] + b[tid]; tid += 1; // 由於只有一個CPU,所以每次遞增1 } } int main( void ) { int a[N], b[N], c[N]; // 在CPU上為陣列a,b賦初值 for (int i=0; i<N; i++) { a[i] = -i; b[i] = i * i; } add( a, b, c ); // display the results for (int i=0; i<N; i++) { printf( "%d + %d = %d\n", a[i], b[i], c[i] ); } return 0; }
我們對add()函式做簡要分析，看看為什麼這個函式有些過於複雜：
函式在while迴圈中計算總和，其中索引 tid 的取值範圍為0到N-1。我們將a[ ] 和 b[ ] 的對應元素相加起來，並將結果儲存在c[ ] 的相應元素中。通常，可以用更簡單的方式來編寫這段程式碼，例如：
void add(int *a, *b, *c){
for(i=0; i<N; i++)
c[i]=a[i]+b[i];
}然而，上面採用while迴圈的方式雖然有些複雜，但這是為了使程式碼能夠在擁有多個CPU或者CP核的系統上並行執行。例如，在雙核處理器上可以將每次遞增的大小改為2，這樣其中一個核從 tid=0 開始執行迴圈 ，而另一個核從 tid=1開始執行迴圈。第一個核將偶數索引的元素相加，而第二個核將奇數索引的元素相加。著相當於在每個CPU核上執行以下程式碼：
// 第1個CPU核心
void add(int *a, int *b, int *c){
int tid = 0; // 第1個核從tid=0開始迴圈
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += 2; // 有兩個核, 所以每次累加2
}
}
// 第2個CPU核心
void add(int *a, int *b, int *c){
int tid = 1; // 第2個核從tid=1開始迴圈
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += 2; // 有兩個核, 所以每次累加2
}
}當然，要在CPU上世紀執行這個運算，還需要增加更多的工作，例如編寫一定數量的程式碼來建立工作執行緒，每個執行緒都執行函式add()，並假設每個執行緒都將並行執行。但是，這種假設是理想而又不實際的，執行緒排程機制的實際執行情況往往並非如此。
4.1.2 基於GPU的向量求和
我們可以在GPU上實現相同的加法運算，這需要將 add() 編寫為一個裝置函式。下面首先給出main()函式：
/* add_loop_gpu.cu */
#include "../common/book.h"
#define N 10
int main( void ) {
int a[N], b[N], c[N];
int *dev_a, *dev_b, *dev_c;
// 在GPU上分配記憶體
HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_c, N * sizeof(int) ) );
// 在CPU上為陣列a, b賦初值
for (int i=0; i<N; i++) {
a[i] = -i;
b[i] = i * i;
}
// 將陣列a,b從CPU複製到GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(int),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(int),
cudaMemcpyHostToDevice ) );
add<<<N,1>>>( dev_a, dev_b, dev_c );
// 將陣列c從GPU複製到CPU
HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(int),
cudaMemcpyDeviceToHost ) );
// display the results
for (int i=0; i<N; i++) {
printf( "%d + %d = %d\n", a[i], b[i], c[i] );
}
// 釋放在GPU上分配的記憶體
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_c ) );
return 0;
}上述程式碼使用了一些通用模式：
- 呼叫cudaMalloc() 在裝置上為三個陣列分配記憶體： 在其中兩個陣列（dev_a和dev_b）中包含了輸入值，而在陣列dev_c中包含了計算結果；
- 為了避免記憶體洩漏，在使用完GPU記憶體後通過cudaFree() 釋放他們；
- 通過cudaMemcpy()將輸入資料複製到裝置中，同時指定引數 cudaMemcpyHostToDevice，在計算完成後，將計算結果通過引數cudaMemcpyDeviceToHost是複製回主機。
- 通過尖括號語法， 在主機main() 函式中執行add()裝置程式碼。
上面程式碼中在CPU中對陣列 a, b 賦初值只是為了說明如何在GPU上實現兩個向量的加法運算，事實上，如果在GPU上為陣列賦值，這個步驟會執行得更快。
接下來是add()函式，這個函式看上去非常類似於基於CPU實現的add() :
__global__ void add( int *a, int *b, int *c ) {
int tid = blockIdx.x; // 計算該索引處的資料
if (tid < N)
c[tid] = a[tid] + b[tid];
}上面add()函式使用了一種通用模式：
- 編寫一個在裝置上執行的函式add。採用C來編寫程式碼，並在函式名前新增修飾符__global__ .
前面提到了尖括號中的兩個引數將傳遞給執行時，作用是告訴執行時如何啟動核函式，明確對執行時引數進行解釋見CUDA和函式執行引數。
4.2 一個有趣的示例
接下來的示例將介紹如何繪製Julia集的曲線，對於不熟悉Julia集的讀者，可以簡單地將Julia集認為是滿足某個複數計算函式的所有點構成的邊界。
生成Julia集的演算法非常簡單。Julia集的基本演算法是，通過一個簡單的迭代等式對複平面中的點求值。如果在計算某個點時，迭代等式的計算結果是發散的，那麼這個點就不屬於Julia集合。相反，如果在迭代等式中計算得到的一系列值都位於某個邊界範圍之內，那麼這個點就屬於Julia集合。迭代等式如下：
迭代過程包括：首先計算當前值的平方，然後再加一個常數以得到下一個值。
基於CPU的Julia集
/* julia_cpu.cu */
#include "../common/book.h"
#include "../common/cpu_bitmap.h"
#define DIM 1000
struct cuComplex {
float r;
float i;
cuComplex( float a, float b ) : r(a), i(b) {}
float magnitude2( void ) { return r * r + i * i; }
cuComplex operator*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
cuComplex operator+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
};
int julia( int x, int y ) {
const float scale = 1.5;
float jx = scale * (float)(DIM/2 - x)/(DIM/2);
float jy = scale * (float)(DIM/2 - y)/(DIM/2);
cuComplex c(-0.8, 0.156);
cuComplex a(jx, jy);
int i = 0;
for (i=0; i<200; i++) {
a = a * a + c;
if (a.magnitude2() > 1000)
return 0;
}
return 1;
}
void kernel( unsigned char *ptr ){
for (int y=0; y<DIM; y++) {
for (int x=0; x<DIM; x++) {
int offset = x + y * DIM;
int juliaValue = julia( x, y );
ptr[offset*4 + 0] = 255 * juliaValue;
ptr[offset*4 + 1] = 0;
ptr[offset*4 + 2] = 0;
ptr[offset*4 + 3] = 255;
}
}
}
int main( void ) {
CPUBitmap bitmap( DIM, DIM );
unsigned char *ptr = bitmap.get_ptr();
kernel( ptr );
bitmap.display_and_exit();
}
基於GPU的Julia集
/* julia_gpu.cu */
#include "../common/book.h"
#include "../common/cpu_bitmap.h"
#define DIM 1000
struct cuComplex {
float r;
float i;
cuComplex( float a, float b ) : r(a), i(b) {}
__device__ float magnitude2( void ) {
return r * r + i * i;
}
__device__ cuComplex operator*(const cuComplex& a) {
return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
}
__device__ cuComplex operator+(const cuComplex& a) {
return cuComplex(r+a.r, i+a.i);
}
};
__device__ int julia( int x, int y ) {
const float scale = 1.5;
float jx = scale * (float)(DIM/2 - x)/(DIM/2);
float jy = scale * (float)(DIM/2 - y)/(DIM/2);
cuComplex c(-0.8, 0.156);
cuComplex a(jx, jy);
int i = 0;
for (i=0; i<200; i++) {
a = a * a + c;
if (a.magnitude2() > 1000)
return 0;
}
return 1;
}
__global__ void kernel( unsigned char *ptr ) {
// map from blockIdx to pixel position
int x = blockIdx.x;
int y = blockIdx.y;
int offset = x + y * gridDim.x;
// now calculate the value at that position
int juliaValue = julia( x, y );
ptr[offset*4 + 0] = 255 * juliaValue;
ptr[offset*4 + 1] = 0;
ptr[offset*4 + 2] = 0;
ptr[offset*4 + 3] = 255;
}
// globals needed by the update routine
struct DataBlock {
unsigned char *dev_bitmap;
};
int main( void ) {
DataBlock data;
CPUBitmap bitmap( DIM, DIM, &data );
unsigned char *dev_bitmap;
HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );
data.dev_bitmap = dev_bitmap;
dim3 grid(DIM,DIM);
kernel<<<grid,1>>>( dev_bitmap );
HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
bitmap.image_size(),
cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaFree( dev_bitmap ) );
bitmap.display_and_exit();
}
4.3 本章小結
到目前為止，我們已經看到了如何告訴CUDA執行時線上程塊上並行執行程式。我們把在GPU上啟動的執行緒塊集合稱為一個執行緒格。從名字的含義可以看出，執行緒格既可以是一位的執行緒塊集合，也可以是二維的執行緒塊集合。核函式的每個副本都可以通過內建變數blockIdx來判斷那個執行緒塊正在執行它。同樣，它還可以通過內建變數gridDim來獲得執行緒格的大小。這兩個內建變數在核函式中都是非常有用的，可以用來計算每個執行緒塊需要的資料索引。