《CUDA By Example》【Chapter 05】執行緒協作 ？
阿新 • • 發佈：2019-02-17
概述
本章介紹程式碼在各個並行副本之間的通訊和協作。
1，瞭解不同執行緒之間的通訊機制；
2，瞭解並行執行執行緒的同步機制；
5.2 並行執行緒塊的分解
add<<<N,1>>>( dev_a, dev_b, dev_c);
//第一個引數為想要啟動的執行緒塊數量；
//第二個引數為CUDA執行時在每個執行緒塊中建立的執行緒數量；
//共建立了N個執行緒塊，N*1個並行執行緒；add<<N,2>>>( dev_a, dev_b, dev_c);
//共建立了N個執行緒塊，2N個並行執行緒；索引的使用（執行緒塊索引&執行緒索引）
add<<<N,1>>>( dev_a, dev_b, dev_c);
//此時使用了N個執行緒塊，每個執行緒塊1個執行緒；
//需要使用 執行緒索引: blockIdx.x
add<<<1,N>>>( dev_a, dev_b, dev_c)；
//此時使用了1個執行緒塊，1個執行緒塊有N個執行緒；
//需要使用 執行緒索引: threadIdx.x從並行執行緒塊到並行執行緒需要進行的兩處修改
1，blockIdx.x ==> threadIdx.x
2，add<<<N,1>>>() ==> add 使用執行緒實現GPU上的向量求和（單純使用執行緒索引）
#include "../common/book.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 在GPU上對更長的向量求和
因為硬體限制，執行緒塊的數量有上限，對於啟動核函式時每個執行緒塊中的執行緒數量，也有上限。具體地說，最大的執行緒數量不能超過裝置屬性結構中maxThreadsPerBlock域的值。對很多的GPU，這個限制值是每個執行緒塊最多512個執行緒。 對於並行執行緒大於512的向量相加，需結合線程和執行緒塊來計算。
從上例，改動兩個地方：核函式中的索引計算方法和核函式的呼叫方式
核函式中的索引計算方法
int tid = threadIdx.x + blockIdx.x * blockDim.x;
//blockDim, 對所有執行緒塊是個常數，為執行緒塊中每一維的執行緒數量。因為使用的一維執行緒塊，所以只使用blockDim.x。在chapter 4 中， gridDim有個類似的值，即線上程格中每一維的執行緒塊數量。gridDim是二維的，而blockDim實際上是三維的。CUDA執行時允許啟動一個二維執行緒格，並且執行緒格中的每個執行緒塊都是一個三維的執行緒陣列。不常用，但是支援。
核函式中的索引計算方法
對於N個執行緒的核函式呼叫如下。如每個執行緒塊有128個執行緒。
add<<<(N+127)/128, 128>>>( dev_a, dev_b, dev_c );在GPU上對任意長度的向量求和
概述
執行緒塊的數量也有硬體限制，如執行緒格每一維大小不能超過65535。
#include "../common/book.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, *b, *c;
int *dev_a, *dev_b, *dev_c;
// allocate the memory on the CPU
a = (int*)malloc( N * sizeof(int) );
b = (int*)malloc( N * sizeof(int) );
c = (int*)malloc( N * sizeof(int) );
// allocate the memory on the 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) ) );
// fill the arrays 'a' and 'b' on the CPU
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = 2 * i;
}
// copy the arrays 'a' and 'b' to the GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(int),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(int),
cudaMemcpyHostToDevice ) );
add<<<128,128>>>( dev_a, dev_b, dev_c );
// copy the array 'c' back from the GPU to the CPU
HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(int),
cudaMemcpyDeviceToHost ) );
// verify that the GPU did the work we requested
bool success = true;
for (int i=0; i<N; i++) {
if ((a[i] + b[i]) != c[i]) {
printf( "Error: %d + %d != %d\n", a[i], b[i], c[i] );
success = false;
}
}
if (success) printf( "We did it!\n" );
// free the memory we allocated on the GPU
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_c ) );
// free the memory we allocated on the CPU
free( a );
free( b );
free( c );
return 0;
}5.2.2 在GPU上使用執行緒實現濾波效果
#include "cuda.h"
#include "../common/book.h"
#include "../common/cpu_anim.h"
#define DIM 1024
#define PI 3.1415926535897932f
__global__ void kernel( unsigned char *ptr, int ticks ) {
// 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;
//????????? offset為什麼這麼計算？？？？？？
// now calculate the value at that position
float fx = x - DIM/2;
float fy = y - DIM/2;
float d = sqrtf( fx * fx + fy * fy );
unsigned char grey = (unsigned char)(128.0f + 127.0f *
cos(d/10.0f - ticks/7.0f) /
(d/10.0f + 1.0f));
ptr[offset*4 + 0] = grey;
ptr[offset*4 + 1] = grey;
ptr[offset*4 + 2] = grey;
ptr[offset*4 + 3] = 255;
}
struct DataBlock {
unsigned char *dev_bitmap;
CPUAnimBitmap *bitmap;
};
void generate_frame( DataBlock *d, int ticks ) {
dim3 blocks(DIM/16,DIM/16); //？？？？？
dim3 threads(16,16); //？？？？？
kernel<<<blocks,threads>>>( d->dev_bitmap, ticks );
HANDLE_ERROR( cudaMemcpy( d->bitmap->get_ptr(),
d->dev_bitmap,
d->bitmap->image_size(),
cudaMemcpyDeviceToHost ) );
}
// clean up memory allocated on the GPU
void cleanup( DataBlock *d ) {
HANDLE_ERROR( cudaFree( d->dev_bitmap ) );
}
int main( void ) {
DataBlock data;
CPUAnimBitmap bitmap( DIM, DIM, &data );
data.bitmap = &bitmap;
HANDLE_ERROR( cudaMalloc( (void**)&data.dev_bitmap,
bitmap.image_size() ) );
bitmap.anim_and_exit( (void (*)(void*,int))generate_frame,
(void (*)(void*))cleanup );
}5.3 共享記憶體和同步
執行緒塊分解為執行緒，一方面是為了解決執行緒塊數量的硬體限制。還有一個原因是便於利用共享記憶體實現執行緒同步。
對於GPU上啟動的每一個執行緒塊，CUDA C都將建立變數的一個副本。執行緒塊中的每個執行緒都共享這塊記憶體，但執行緒卻無法看到和修改其他執行緒塊的變數副本。這便可以實現執行緒塊內的執行緒通訊。 同時因為shared memory在GPU內部， 訪問速度可以很快。
使用執行緒塊中的共享記憶體，需要注意同步問題。
5.3.1 點積運算（dot product, 內積）
通常來說，歸約運算不能充分利用GPU的並行運算能力，一般會傳回CPU來運算。本例依舊使用GPU完成歸約計算。
#include "../common/book.h"
#define imin(a,b) (a<b?a:b)
const int N = 33 * 1024;
const int threadsPerBlock = 256;
const int blocksPerGrid =
imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );
__global__ void dot( float *a, float *b, float *c ) {
__shared__ float cache[threadsPerBlock];
//共享記憶體緩衝區；因為每個執行緒塊生成共享變數的一個副本，所以只需要根據執行緒塊中的執行緒的數量來分配記憶體
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
//共享的記憶體的索引和執行緒索引相同
float temp = 0;
while (tid < N) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
//可以解決任意長度的向量計算
}
// set the cache values
cache[cacheIndex] = temp;
// synchronize threads in this block
__syncthreads(); //確保執行緒塊中每個執行緒都執行完該函式前面的語句後，才會執行下一條語句
// for reductions, threadsPerBlock must be a power of 2
// because of the following code
//歸約運算
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
}
int main( void ) {
float *a, *b, c, *partial_c;
float *dev_a, *dev_b, *dev_partial_c;
// allocate memory on the cpu side
a = (float*)malloc( N*sizeof(float) );
b = (float*)malloc( N*sizeof(float) );
partial_c = (float*)malloc( blocksPerGrid*sizeof(float) );
// allocate the memory on the GPU
HANDLE_ERROR( cudaMalloc( (void**)&dev_a,
N*sizeof(float) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_b,
N*sizeof(float) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c,
blocksPerGrid*sizeof(float) ) );
// fill in the host memory with data
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = i*2;
}
// copy the arrays 'a' and 'b' to the GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float),
cudaMemcpyHostToDevice ) );
dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b,
dev_partial_c );
// copy the array 'c' back from the GPU to the CPU
HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
blocksPerGrid*sizeof(float),
cudaMemcpyDeviceToHost ) );
// finish up on the CPU side
c = 0;
for (int i=0; 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,
2 * sum_squares( (float)(N - 1) ) );
// free memory on the gpu side
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_partial_c ) );
// free memory on the cpu side
free( a );
free( b );
free( partial_c );
}5.3.2 （不正確的）點積運算優化
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}為了充分利用執行緒塊中的空閒執行緒， 對上面程式碼段進行下面的優化。 會發生執行緒發散（thread divergence），即當某些執行緒需要執行一條指令，而其他執行緒不需要執行時，這種情況就是執行緒發散。
正常情況下，發散的分支只會引起某些執行緒處於空閒狀態，而其他執行緒繼續執行分支中的程式碼。但是在__syncthreads()中，執行緒發散可能造成hang機。
所以，執行緒同步__syncthreads()需要保證不線上程發散中。
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
{
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
}
i /= 2;
}5.3.3 基於共享記憶體的點陣圖
#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[16][16];
// now calculate the value at that position
const float period = 128.0f;
shared[threadIdx.x][threadIdx.y] =
255 * (sinf(x*2.0f*PI/ period) + 1.0f) *
(sinf(y*2.0f*PI/ period) + 1.0f) / 4.0f;
// removing this syncthreads shows graphically what happens
// when it doesn't exist. this is an example of why we need it.
__syncthreads();
/如果沒有這個同步，輸出的影象會出現瑕疵
ptr[offset*4 + 0] = 0;
ptr[offset*4 + 1] = shared[15-threadIdx.x][15-threadIdx.y];
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 grids(DIM/16,DIM/16);
dim3 threads(16,16);
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();
}