線程協作
目标
- 了解CUDA C中的線程
- 了解不同線程之間的通信機制
- 了解并執行線程的同步機制
并行線程塊的分解
前面我們了解了如何在GPU上啟動并行代碼。具體的方法是，告訴CUDA運作時啟動核函數的多個并行副本。這裡将這些并行副本稱為線程塊(Block)。
CUDA運作時将這些線程塊分解為多個線程，當需要啟動多個并行線程塊時，隻需要将尖括号中的第一個參數由1改為想要啟動的線程塊的數量。例如在介紹矢量加法示例的時候，曾經就是為矢量中的每個元素(共有N個元素)都啟動一個線程塊。調用如下：
add<<<N,>>>(dev_a,dev_b,dev_c);
在尖括号中，第二個參數表示CUDA運作時在每個線程塊中建立的線程數量。上面這行代碼表示每個線程塊中隻啟動一個線程。總共啟動線程的數量可以按下面公式計算：
N個線程塊*1個線程/線程=N個并行線程
事實上我們也可以啟動N/2個線程塊，每個線程塊包含兩個線程，或者N/4個線程塊，每個線程塊包含4個線程。
矢量求和：重新回顧
這裡将實作一個和前面相同的任務，即将兩個輸入矢量相加并将得到的結果儲存在第三個輸出矢量中。這裡将使用線程而不是線程塊來實作這一工作。這兩種方法現階段沒有任何差別，但是，線程塊中的并行線程能夠完成并行線程塊無法完成的工作。
使用線程實作GPU上的矢量求和
在将實作方式由并行線程塊改為并行線程時，首先需要修改兩個地方。第一個地方是，之前核函數的調用中是啟動N個線程塊，并且每個線程塊對應一個線程，具體如下所示：
add<<<N,>>>(dev_a,dev_b,dev_c);
現在我們将啟動N個線程，所有的線程都在一個線程塊中：
add<<<,N>>>(dev_a,dev_b,dev_c);
第二個要修改的地方就是對資料的索引的方法。之前核函數中是通過線程塊索引對輸入資料和輸出資料進行索引。
這裡由于隻有一個線程塊，是以将采用線程索引對資料進行索引。
是以采用并行線程實作GPU上的矢量求和的例子：
#include <stdio.h>
#define N 10
__global__ void add(int *a,int *b,int *c)
{
int tid = threadIdx.x;
if(tid < N){
c[tid] = a[tid] + b[tid];
}
}
int main(void)
{
int a[N],b[N],c[N];
int *dev_a,*dev_b,*dev_c;
//在GPU上配置設定記憶體
cudaMalloc((void**)&dev_a,N*sizeof(int));
cudaMalloc((void**)&dev_b,N*sizeof(int));
cudaMalloc((void**)&dev_c,N*sizeof(int));
//在CPU上為數組‘a’和數組‘b’指派
for(int i=;i<N;i++){
a[i] = i;
b[i] = i*i;
}
//将數組‘a’和數組‘b’複制到GPU上
cudaMemcpy(dev_a,a,N*sizeof(int),cudaMemcpyHostToDevice);
cudaMemcpy(dev_b,b,N*sizeof(int),cudaMemcpyHostToDevice);
add<<<,N>>>(dev_a,dev_b,dev_c);
//将數組‘C’從GPU複制到CPU中
cudaMemcpy(c,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost);
//顯示結果
for(int i=;i<N;i++){
printf("%d+%d=%d\n",a[i],b[i],c[i]);
}
//釋放記憶體
cudaFree(&dev_a);
cudaFree(&dev_b);
cudaFree(&dev_c);
}
在GPU上求更長矢量的和
在前面我們了解到硬體将線程塊的數量限制為不超過65535個。同樣對于啟動核函數時每個線程塊中的線程數量，硬體也做了限制。具體的來說，最大線程數量不能超過裝置屬性中的maxThreadsPerBlock域的值，是以對于一些比較大的矩陣就需要将線程與線程塊結合起來實作計算。這種方法主要需要修改兩個地方：
* 核函數的索引計算方法
int tid = threadIdx.x + blockIdx.x*blockDim.x
這裡面使用了一個新的變量
blockDim
。對于所有的線程來說這個變量是一個常數，儲存的是線程塊中每一維的線程的數量。由于使用的是使用一維線程線程塊，是以隻需要用到
blockDim.x
。回顧第四章的内容，在
gridDim
中儲存了一個類似的值，即線上程格每一維的線程塊數量。此外
gridDim
是二維的，而
blockDim
實際上是三維的。也就是說CUDA運作時允許啟動一個二維的線程格，并且線程中每個線程塊都是一個三維的線程數組。
* 核函數的調用方式。這種方法需要修改核函數調用本身。雖然任然需要N個并行線程，但是我們希望在多個線程塊中啟動這些線程，這樣就不會出現每個線程塊最多不會超過512個線程的最大數量的限制。這裡将每個線程塊中包含的線程的數量設定為128。然後，可以啟動N/128個線程塊，這樣總共就啟動了N個線程同時運作。這裡的問題在于，N/128是一個整除，是以這裡我們希望能夠向上取整數。具體核函數調用方法如下：
add<<<(N+127)/128,128>>>(dev_a,dev_b,dev_c)
。這裡會發現将啟動過多的線程。然而，在核函數中已經解決了這一問題：
if(tid<N)
{
c[tid] = a[tid]+b[tid];
}
是以，當索引越過數組邊界時，例如當啟動的并行線程數量不是128的整數倍的時候就會出現這種情況，那麼核函數将自動停止計算。更重要的是，核函數不會對越過數組邊界的記憶體進行存儲或者讀寫操作。
在GPU上對任意長度的矢量求和
當第一次介紹在GPU上啟動并行線程塊的時候，我們并沒有考慮所有可能的情況，除了線上程數量商存在限制之外，線上程塊的數量上同樣存在着一個硬體的限制，線程格的每一個次元的大小不能超過65535。
是以，這就對矢量加法的實作帶來一個問題。如果啟動N/128個線程相加，那麼當矢量長度超過65535 X 128的時候核函數就會調用失敗。這裡需要對核函數做下面的修改：
__global_ void add(int *a,int *b,int *c)
{
int tid = threadIdx.x + blockIdx.x*blockDim.x;
while(tid<N){
c[tid] = a[tid] + b[tid];
tid += blockDim.x*gridDim.x;
}
}
這裡可以看出當使用的線程數量未達到要計算的量，該線程不會在進行下一個計算。在GPU系統中我們将并行線程的數量看作是處理器的數量。盡管GPU處理單元的數量可能小于這個值，但是我們認為每個線程在邏輯單元上都可以并行執行，并且硬體可以排程這個線程以便實際執行。通過并行化過程與硬體的實際執行過程解耦開來。
在了解了上述實作方式所包含的核心思想後，我們隻需要知道如何計算每個并行線程的初始索引，以确定遞增量的值。我們希望每個并行線程從不同的索引開始，是以就需要對線程索引和線程塊索引進行線性化。
int tid = threadIdx.x + blockIdx.x*blockDim.x
。在完成上面的計算後就需要對索引的增量進行計算，這個數值等于每個線程塊中線程的數量乘以線程格中線程塊的數量。
在核函數調用的時候為了保證不啟動過多的線程塊，我們将線程固定為某個較小的值。完整的代碼如下：（這裡計算的矢量長度取決與GPU上的記憶體容量）
#include <stdio.h>
#define N (33*1024)
__global_ void add(int *a,int *b,int *c)
{
int tid = threadIdx.x + blockIdx.x*blockDim.x;
while(tid < N){
c[tid] = a[tid] + b[tid];
tid += blockDim.x*gridDim.x;
}
}
int main(void)
{
int a[N],b[N],c[N];
int *dev_a,*dev_b,*dev_c;
//在GPU上配置設定記憶體
cudaMalloc((void**)&dev_a,N*sizeof(int));
cudaMalloc((void**)&dev_b,N*sizeof(int));
cudaMalloc((void**)&dev_c,N*sizeof(int));
//在CPU上為數組'a'和數組'b'指派
for(int i=;i<N;i++){
a[i] = i;
b[i] = i*i;
}
//将數組‘a’和數組‘b’複制到GPU記憶體中
cudaMemcpy(dev_a,a,N*sizeof(int),cudaMemcpyHostToDevice);
cudaMemcpy(dev_b,b,N*sizeof(int),cudaMemcpyHostToDevice);
add<<<,>>>(dev_a,dev_b,dev_c);
//将數組‘C’從GPU複制到CPU中
cudaMemcpy(c,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost);
//驗證GPU計算結果
for(int i=;i<N;i++){
if(a[i]+b[i] != c[i]){
printf("error:%d+%d=%d\n",a[i],b[i],c[i]);
}
}
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
return ;
}
在GPU上使用線程實作紋波效果
#include "cuda.h"
#include "../common/book.h"
#include "../common/cpu_anim.h"
#define DIM 1024
#define PI 3.1415926535897932
/* 核函數需要知道目前動畫時間以便于生成正确的幀。在CPUAnimBitmap的代碼中
* 将目前時間ticks傳遞給generate_frame()
* /
__global__ void kernel(unsigned char *ptr,int ticks){
int x = threadIdx.x + blockIdx.x*blockDim.x;
int y = threadIdx.y + blockIdx.y*blockDim.y;
int offset = x + y*blockDim.x*gridDim.x;
float fx = x - DIM/;
float fy = y - DIM/;
float d = sqrtf(fx*fx + fy*fy);
unsigned char gray = (unsigned char)(f + f*
cos(d/f - ticks/f)/
(d/f + f));
ptr[offset* + ] = gray;
ptr[offset* + ] = gray;
ptr[offset* + ] = gray;
ptr[offset* + ] = ;
}
struct DataBlock{
unsigned char *dev_bitmap;
CPUAnimBitmap *bitmap;
};
void generate_frame(DataBlock *d,int ticks){
dim3 blocks(DIM/,DIM/);
dim3 threads(,);
kernel<<<blocks,threads>>>(d->dev_bitmap,ticks);
cudaMemcpy(d->bitmap->get_ptr(),
d->dev_bitmap,
d->bitmap->image_size(),
cudaMemcpyDeviceToHost);
}
//釋放在GPU上配置設定的記憶體
void cleanup(DataBlock *d){
cudaFree(d->dev_bitmap);
}
int main(void)
{
DataBlock data;
CPUAnimBitmap bitmap(DIM,DIM,&data);
data.bitmap = &bitmap;
cudaMalloc((void**)&data.dev_bitmap,bitmap.image_size());
/*
* 這裡将一個指向generate_frame()的函數指針傳遞給anim_and_exit()。每當要生成一副新的動畫
* 時都要調用函數generate_frame();
*/
bitmap.anim_and_exit((void(*)(void*,int))generate_frame,//函數的指針調用
(void(*)(void*))cleanup);
return ;
}
共享記憶體和同步
到目前為止将線程塊分解為線程的目的隻是為了解決線程塊數量的硬體限制。這是一個很勉強的動機，因為這個完全可以由CUDA運作時在幕後完成。顯然還有其他原因需要将線程分解為多個線程。
CUDA C支援共享記憶體。這塊記憶體引出了對C語言的另一個擴充，這個擴充類似于
__device__
和
__global__
。在編寫代碼的時候，你可以将CUDA C的關鍵字
__share__
添加到變量的聲明中，使得這個變量存儲在共享記憶體中。那麼這樣做的目的是什麼了？
CUDA C編譯器對共享記憶體中的變量和普通變量将分别采取不同的處理方式。對于在GPU上自動啟動的多個線程塊，CUDA C編譯器都将建立該變量的一個副本。線程塊中的每個線程都将共享這塊記憶體，但是線程中無法看到也不能修改其他線程塊變量副本。這就實作了一種非常好的方式，使得一個線程塊中的多個線程能夠在計算上進行通信和協作。而且共享記憶體的緩存區駐留在實體GPU上，而不是駐留在GPU之外的系統記憶體中。是以，在通路共享記憶體時的延遲要遠遠低于通路普通緩沖區的延遲，使得共享記憶體像在每個線程塊的高速緩存或者中間結果暫存器上那樣高效。
線程之間進行通訊，那麼還需要一種機制來實作縣城之間的同步。例如：如果線程A将一個值寫入到共享記憶體，并且我們希望線程B對這個值進行一些操作，那麼隻有當線程A将寫入操作完成後，線程B才能開始執行它的操作。如果沒有同步，那麼将發生競态條件，在這種情況下，代碼執行結果的正确性将取決硬體的不确定性。
點積運算
矢量的點積運算計算公式如下：
(x1,x2,x3,x4)·(y1,y2,y3,y4)=x1y1+x2y2+x3y3+x4y4
這裡可以像矢量加法那樣，每個線程将兩個相應的元素相乘，然後移動到下兩個元素。由于最終的結果是所有乘積的總和，是以每個線程還要儲存它所計算的成積的加和。與矢量加法類似，線程每次對索引的增加值為線程的數量，這樣能確定不會遺漏任何元素。而且也不會将任何元素相乘兩次。具體代碼如下：
#include <stdio.h>
#define imin(a,b) (a<b?a:b)
const int N = *;
const int threadsPerBlock = ;
const int blocksPerGrid = imin(,(N+threadsPerBlock-)/threadsPerBlock);
__global__ void dot(float *a,float *b,float *c)
{
/*
* 代碼在這裡聲明了一個記憶體共享緩沖區，名字為cache，這個緩存區儲存每個
* 線程計算的加和值，這裡将數組大小聲明為threadsPerBlock，這樣線程塊中的
* 每個線程都能被配置設定足夠的記憶體，即線程塊的數量乘以threadsPerBlock。但對
* 于共享變量來說，由于編譯器将為每個線程塊生成一個共享變量的副本，是以
* 隻需要根據線程塊中線程的數量來配置設定記憶體。
*
*/
__shared__ float cache[threadsPerBlock];
/*
* 在配置設定共享記憶體後這裡通過線程塊索引和線程索引計算輸出數組中的一個全局偏移。
* 共享記憶體緩存中的偏移就等于線程索引。線程塊的索引與這個偏移無關，因為每個
* 線程塊都有該共享記憶體的私有副本
*/
int tid = threadIdx.x + blockIdx.x*blockDim.x
int cacheIndex = threadIdx.x;
float temp = ;
while(tid < N){
temp += a[tid] * b[tid];
tid += blockDim.x*gridDim.x
}
//設定catch中相應位置的值。這裡設定了共享記憶體緩沖區相應位置上的值，
//以便随後能将整個數組相加起來，并無需擔心某個特定位置商是否包含有效的資料
catche[cacheIndex] = temp;
/*
* 當算法執行到這一步時，我們需要将cache中所有的值相加起來。在執行這個運算時，
* 需要通過一個線程來讀取儲存在cache中的值，然而前面提到過，這可能是一種危險的
* 操作。我們需要某種方法來確定所有對共享數組cache[]的寫入操作在讀取的cache之前
* 完成。具體為使用__syncthreads();對線程同步。這個函數将確定線程塊中的每個線程都
* 執行完__syncthreads()前面的語句，才會執行後面的第一條語句，這時候，線程塊中的
* 其他線程肯定都執行完了__syncthreads()。
*/
//對線程塊同步
__syncthreads();
//通過歸約運算，計算每個線程塊内的和，這裡要求threadPerBlock必須是2的指數
int i = blockDim.x/;
while(i != ){
if(cacheIndex<i){
catche[catcheIndex] += catche[catcheIndex + i];
__syncthreads();
i /= ；
}
if(cacheIndex == )
c[blockIdx.x] = cache[];
}
}
int main(void)
{
float *a, *b, c, pratial_c;
float *dev_a,*dev_b,*dev_pratial_c;
//在CPU上為數組配置設定記憶體
a = (float*)malloc(N*sizeof(float));
b = (float*)malloc(N*sizeof(float));
partial_c = (float*)malloc(blocksPerGrid*sizeof(float));
//在GPU上配置設定記憶體
cudaMalloc((void**)&dev_a,N*sizeof(float));
cudaMalloc((void**)&dev_b,N*sizeof(float));
cudaMalloc((void**)&dev_pratial_c,blocksPerGrid*sizeof(float));
//填充主機記憶體
for(int i=;i< N;i++){
a[i] = i;
b[i] = i*;
}
//将數組‘a’和數組‘b’複制到GPU上
cudaMemcpy(dev_a,a,N*sizeof(float),cudaMemcpyHostToDevice);
cudaMemcpy(dev_b,b,N*sizeof(float),cudamemcpyHostToDevice);
dot<<<blocksPerGrid,threadsPerBlock>>>(dev_a,dev_b,dev_partial_c);
//将數組‘C’從GPU複制到CPU
cudaMemcpy(partial_c,dev_partial_c,blocksPerGrid*sizeof(float),cudaMemcpyDeviceToHost);
//在CPU上完成最後的求和運算
c = ;
for(int i = ;i<blocksPerGrid;i++){
c += partial_c[i];
}
#define sum_squares(x) (x*(x+1)*(2*x+1)/6)
printf("Does GPU Value %.6g = %.6g?\n",c,*sum_squares((float)(N-)));
//釋放GPU上的記憶體
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_partial_c);
//釋放CPU上的記憶體
free(a);
free(b);
free(partial_c);
}
（錯誤的）點積運算優化
在前面介紹了點積運算示例中的第二個
__syncthreads()
。現在我們将在這裡進一步分析這個函數的調用，并嘗試對這個函數做出改進。之是以需要第二個
__syncthreads()
，是因為在循環疊代中更新了共享記憶體變量
cache()
，并且在循環的下一次疊代開始之前，需要確定目前疊代中所有線程的更新操作都已經完成。
int i = blockDim.x/;
while(i != ){
if(cacheIndex < i){
cache[cacheIndex] += cache[cacheIndex + i];
}
__syncthreads();
i /= ;
}
在上面的程式中，我們發現隻有當cacheIndex小于i的時候才需要更新共享記憶體的緩存區cache[]。由于cacheIndex實際上就是等于threadIdx.x，是以這意味着隻有一部分線程會更新共享記憶體。由于調用
__syncthreads
的目的隻是為了確定這些更新在下一次疊代之前已經完成，是以如果将代碼修改為了等待那些需要寫入共享記憶體的線程，是不是就能獲得性能提升了。
int i = blockDim.x/;
while(i != ){
if(cacheIndex < i){
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
}
i /= ;
}
雖然這種優化代碼的初衷不錯，但是卻不起作用。實際上這種情況比優化之前的情況還要糟糕。再這樣對核函數修改之後，GPU将停止響應，是以不得不強行終止程式。在GPU中的線程塊中的每個線程依次通過代碼，每次一行。每個線程都要執行相同的指令，但是對不同的資料進行運算。然而，當每個線程執行的指令位于一個條件語句中，如
if()
語句中的時候，那麼将會出現什麼情況。顯然這時候不是每個線程都會執行這個指令。當某一些線程需要執行一條指令，而其他線程不需要執行時，這種情況就稱為線程的發散。在
__syncthreads()
情況下。應為CUDA架構将確定，除非線程塊中的每個線程都執行
__syncthreads()
，否則沒有任何線程能執行
__syncthreads()
之後的指令。遺憾的是。如果
__syncthreads()
處于發散分支中，那麼一些線程将永遠無法執行
__syncthreads()
。是以由于要確定每個線程執行完
__syncthreads()
後才能執行後面的語句，是以硬體将使這些線程保持等待。是以。如果在點積運算示例中将
__syncthreads()
放到if語句中，那麼任何cacheIndex大于或等于i的線程将永遠不能執行
__syncthreads()
。這将使處理器挂起，應為GPU在等待某個永遠不可能發生的事件。
基于共享記憶體的位圖
#include "cuda.h"
#include "../common/book.h"
#include "../common/cpu_bitmap.h"
#define DIM 1024
#define PI 3.1415926535897932f
__global__ void kernel( unsigned char *ptr ) {
// map from threadIdx/BlockIdx to pixel position
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
__shared__ float shared[][];
// now calculate the value at that position
const float period = ;
shared[threadIdx.x][threadIdx.y] =
* (sinf(x**PI/ period) + ) *
(sinf(y**PI/ period) + ) / ;
//由于最後是通過逆序複制的不加線程同步會導緻未被計算先給了一個随機的值
__syncthreads();
ptr[offset* + ] = ;
ptr[offset* + ] = shared[-threadIdx.x][-threadIdx.y];
ptr[offset* + ] = ;
ptr[offset* + ] = ;
}
// 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 grids(DIM/,DIM/);
dim3 threads(,);
kernel<<<grids,threads>>>( dev_bitmap );
HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
bitmap.image_size(),
cudaMemcpyDeviceToHost ) );
HANDLE_ERROR( cudaFree( dev_bitmap ) );
bitmap.display_and_exit();
}
本章小結
這裡了解了如何将線程塊進一步分解為更小的并行執行單元。這種并行執行單元也稱為線程。我們回顧前面矢量相加的示例，并介紹如何實作任意長度矢量的加法。我們還給出了一個規約運算的示例，以及如何通過 共享記憶體和同步來實作這個運算。這個示例說明了如何對GPU和CPU進行協作完成運算。最後，我們還給出了當忽略同步的時候給應用程式造成的問題。
到這裡你已經學習了CUDA C的大部分基礎概念，它與标準C在某些方面是相同的，但在大多數方面是不同的。此時，你可以思考你曾經遇到過的問題，以及那些問題可以通過CUDA C的并行實作。後面将講解CUDA C的更多功能以及一些進階的原子函數。
![CUDA學習筆記(5) 原子操作[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)