CUDA編程入門
Hello World
首先一段程序寫個hello world
#include <stdio.h>
__global__ void hello(){
printf("Hello, threadIdx is:%d\n",threadIdx.x);
}
int main(){
hello<<<1,32>>>();
cudaDeviceReset();
}
編譯
nvcc -o main main.cu
運行
./main
結果
Kernels
kernel在cuda中指的是一個函數，當一個kernel被調用的時候，gpu會同時啟動很多個線程來執行這一個kernel，這樣就實現了並行化；每個線程執行這一kernel將通過線程號來對應輸入數據的下標，這樣保證每個thread執行的kernel一樣，但是處理的數據不一樣。
一個kernel在cuda中可以這么定義：
// kernel定義
__global__ void VecAdd(float* A,float* B, float C){
int i = threadIdx.x;
C[i] = A[i] + B[i]
}
// 調用
VecAddd<<<1,N>>>(A,B,C);
上面代碼中調用的時候通過指定<<<1,N>>>前者1表明使用1個block，每個block有N個threads，這樣在函數調用內部通過取內置變量threadIdx來得到thread的索引，以此對應數組下標，實現並行化。
__global__前綴表示在GPU上執行，其他對應的前綴有如下特點：
__global__:在GPU上執行，可以在CPU上被調用，也可以在GPU上被調用__device__:在GPU上執行，只能在GPU上被調用__host__:在CPU上執行，只能在CPU上被調用
實際上，可以先縷一縷GPU中內存的一些關系：
可以看到，GPU的整個內存被划分為了多個grid，每個grid有很多block，而每個block都共享同樣的內存，每個block的多條thread針對共享的內存並行執行不同的操作。
上面的例子中只是針對thredIdx是一維的情況，實際上，為了方便，我們可以直接將threadIdx變成多維的，以方便我們對矩陣或者更高維的數據的處理，這跟C語言里的二維數組等概念是一樣的，就是用x和y分別代表一個維度的下標。
__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]){
int x = threadIdx.x;
int y = threadIdx.y;
C[x][y] = A[x][y] + B[x][y];
}
int main(){
int num_blocks = 1;
dim3 threads_per_block(N,N);
MatAdd<<<num_blocks,threads_per_block>>>(A,B,C);
}
這個例子是相對很容易理解的了。需要知道的是，每個block的threads數目是有最大限制的，不能超過一定的值，書上說現在的GPU每個block最多是1024個thread。所以上面的例子其實並不能處理很大的矩陣，那么要處理很大的矩陣的時候該怎么辦呢？
可以用上block：
___global___ void MatAdd(float A[N][N], float B[N][N], float C[N][N]){
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
C[i][j] = A[i][j] + B[i][j];
}
int main(){
dim3 threads_per_block(16,16);
dim3 num_blocks(N/num_blocks.x,N/num_blocks.y);
MatAdd(A,B,C);
}
這個例子充分利用了block數量，這樣每個block是16*16也就是256個threads，假設N = 32，那么num_blocks = (2,2)也就是用了4個block。上面blockIdx.x取值是[0,1]，blockDim.x取值是16,threadIdx.x取值是[0,15].
假設每個block最多可以用1024個thread，那么其實可以處理$N_{max} = 32 * \sqrt{num_blocks} $階矩陣。
那這里就有個疑問了，不是說一個block一塊內存嗎，不同block之間不share內存，那這個代碼感覺所有block都share同一塊內存的樣子呀？其實這里跟share不share沒任何關系。。每個thread都在獨立計算一條語句，大家根本沒有說有那種你先我后的協作關系，所以談不上share，所以大家都是independent的，從global memory里讀入，然后寫到global memory。
Compute Capability
一般用x.y來表示一款GPU的計算能力，x相同的表示架構相同，y只是在這個架構上做了一些increasement。
現在GPU與計算能力的x之間有下圖這樣的對應：
| 7 | Volta |
|---|---|
| 6 | Pascal |
| 5 | Maxwell |
| 3 | Kepler |
| 2 | Fermi |
| 1 | Tesla |
這里注意不要與cuda版本號搞混了，計算能力描述的是gpu，而非軟件版本。
Memory Allocation
這節涉及之前省略的，之前的代碼沒有考慮A、B、C的內存申請，這一節討論的就是如何把host的內存(cpu data)轉到device里(gpu data)。
如果要申請一塊連續的GPU內存空間，直接可以通過調用cudaMalloc()來實現，調用cudaFree()釋放空間，那么把host上內存copy到gpu上是通過調用cudaMemcpy()來實現。
__global__ void VecAdd(float* A, float* B, float* C){
int i = blockIdx.x * blockDim.x + threadIdx.x;
// 根據下面的分析，i是可能大於等於N的，當N>threadsPerBlock且為threadsPerBlock的非整數倍的時候，這時候會引入多一個block來運算，其中多余的threads不應該參與運算，所以下面if需要判斷。
if (i < N){
C[i] = A[i] + B[i];
}
}
int main(){
int N = 32;
size_t size = N * sizeof(float);
float* h_A = (float*)malloc(size);
float* h_B = (float*)malloc(size);
float* h_C = (float*)malloc(size);
float* d_A;
float* d_B;
float* d_C;
cudaMalloc(&d_A,size);
cudaMalloc(&d_B,size);
cudaMalloc(&d_C,size);
cudaMemcpy(d_A, h_A,size,cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B,size,cudaMemcpyHostToDevice);
int threadsPerBlock = 256;
// 由於threadsPerBlock通常是定的，所以我們需要調整的一般是Blocks的數量
// 考慮到如果N < threadsPerBlock,那么N/threadsPerBlock就是0，但這是是少需要一個block
// 如果N = threadsPerBlock，這時候如果(N + threadsPerBlock)/threadsPerBlock就是2，但其實這時候只需要一個block
// 綜上，應該是下式才滿足條件
int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
VecAdd<<<blocksPerGrid,threadsPerBlock>>>(d_A,d_B,d_C,N);
cudaMemcpy(h_C,d_C,cudaMemcpyDeviceToHost);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
}
cudaMalloc ( void** devPtr, size_t size )
Parameters
devPtr- Pointer to allocated device memory
size- Requested allocation size in bytes
第一個參數是地址的指針，顯然cudaMalloc只能分配線性內存，線性內存訪問起來是相當塊的。在傳統C語言里，我們如果要取用malloc申請二維數組，可能先申請一些一維數組用來存儲首地址指針，每個首地址再去指向一個申請的一維數組，實現二維數組，在CUDA里可以通過cudaMallocPitch來申請二維數組，得到的地址仍然是線性的，因為訪問塊。
cudaMallocPitch ( void** devPtr, size_t* pitch, size_t width, size_t height )
Parameters
-
devPtr- Pointer to allocated pitched device memory
-
pitch- Pitch for allocation
-
width- Requested pitched allocation width (in bytes)
-
height- Requested pitched allocation height
pitch參數表示的是“寬度”，可以理解為一行有多少個元素乘以sizeof(float)，是返回值。
書上的一個例子為：
__global__ Mykernel(float* devPtr,size_t pitch, int width, int height){
for(int r = 0;r < height; ++r){
float* row = (float*)((char*)devPtr + r * pitch); // 轉char* 是因為原來的devPtr是float*的，加法運算會考慮stride。
for(int c=0;c<width;++c){
float val = row[c];
}
}
}
int width = 64,height = 64;
float* devPtr;
size_t pitch; // sizeof(float) * width
cudaMallocPitch(&devPtr,&pitch,width * sizeof(float),height); // height不用乘sizeof(float)是因為第一維的指針是void*
Mykernel<<<100,512>>>(devPtr,pitch,width,height);
Shared Memory
本節通過一個矩陣相乘的例子，來說明通過block內shared memory來減少從global memory讀取和寫入數據的次數，從而加速矩陣相乘的速度。
首先是沒有shared memory的代碼實現：
typedef struct{
int width;
int height;
float* elements;
} Matrix;
#define BLOCK_SIZE 16
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);
void MatMul(const Matrix A, const Matrix B, Matrix C){
Matrix d_A;
d_A.width = A.width;
d_A.height = A.height;
size_t size = d_A.width * d_A.height * sizeof(float);
cudaMalloc(&d_A.elements,size);
cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
Matrix d_B;
d_B.width = B.width;
d_B.height = B.height;
size_t size = d_B.width * d_B.height * sizeof(float);
cudaMalloc(&d_B.elements,size);
cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
Matrix d_C;
d_C.width = C.width;
d_C.height = C.height;
size_t size = d_C.width * d_C.height * sizeof(float);
cudaMalloc(&d_C.elements,size);
dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
// 列索引和行索引划分
dim3 dimGrid(B.width / dimBlock.x,A.height / dimBlock.y);
MatMulKernel<<<dimBlock,dimGrid>>>(d_A,d_B,d_C);
cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){
float Cvalue = 0;
// 這里求對應A矩陣的第row行以及B矩陣的第j列
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
for(int i=0;i<A.width;i++){
Cvalue += A.elements[row * A.width + i] * B.elements[i * B.width + col];
}
C.elements[row * C.width + col] = Cvalue;
}
代碼很好理解，每個thread負責通過一個循環來計算輸出矩陣C中的一個數值，需要注意的是，C矩陣中的每個數值計算，都需要一個thread從global memory中讀取數據，計算完成，然后寫回。下圖表達的就是這樣的一個計算實例：
所以總的讀取和寫入次數是非常頻繁的，為了減少訪問global memory的次數，可以對C矩陣分塊求取，C矩陣中的每個分塊矩陣由一個block負責計算，這block一次將數據讀入block內的共享內存，然后thread在共享內存上計算，最后寫入到global memory，這大大減少了訪問global memory的次數，加速了計算。
// 加的這個stride就是為了表示原矩陣width，這里的width和height可以是子矩陣的width和height
// elements是子矩陣的首地址
// M(row,col) = *(M.elements + row * M.stride + col)
typedef struct{
int width;
int height;
int stride;
float* elements;
} Matrix;
#define BLOCK_SIZE 16
__device__ void SetElement(Matrix A, int row, int col, float value){
// 雖然A是一個局部變量，但是A里面的elements指向的東西沒被深拷貝
A.elements[row * A.stride + col] = value; // 乘以A的stride是因為A跨行的stride和原矩陣一致
}
__device__ float GetElement(Matrix A,int row, int col){
return A.elements[row * A.stride + col];
}
__device__ Matrix getSubMatrix(Matrix A,int row, int col){
Matrix Asub;
Asub.width = BLOCK_SIZE;
Asub.height = BLOCK_SIZE;
Asub.stride = A.stride;
Asub.elements = A.elements + (BLOCK_SIZE * row) * stride + (BLOCK_SIZE*col) // 相當於乘以BLOCK_SIZE
return Asub;
}
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);
void MatMul(const Matrix A, const Matrix B, Matrix C){
Matrix d_A;
d_A.width = A.width;
d_A.height = A.height;
size_t size = d_A.width * d_A.height * sizeof(float);
cudaMalloc(&d_A.elements,size);
cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);
Matrix d_B;
d_B.width = B.width;
d_B.height = B.height;
size_t size = d_B.width * d_B.height * sizeof(float);
cudaMalloc(&d_B.elements,size);
cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);
Matrix d_C;
d_C.width = C.width;
d_C.height = C.height;
size_t size = d_C.width * d_C.height * sizeof(float);
cudaMalloc(&d_C.elements,size);
dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
// 列索引和行索引划分
dim3 dimGrid(B.width / dimBlock.x,A.height / dimBlock.y);
MatMulKernel<<<dimBlock,dimGrid>>>(d_A,d_B,d_C);
cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){
int block_row = blockIdx.x;
int block_rol = blockIdx.y;
Matrix Csub = GetSubMatrix(C,block_row,block_col);
float Cvalue = 0;
int row = threadIdx.x;
int col = threadIdx.y;
for(int i=0;i<(A.width / BLOCK_SIZE);++i){
// 這里不是直接一行算，A矩陣和B矩陣也是按block來，最后是Cvalue求和不影響結果
Matrix Asub = GetSubMatrix(A,block_row,i);
Matrix Bsub = GetSubMatrix(B,i,block_rol);
// 直接生命的對象 局部生命周期
__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
As[row][col] = GetElement(Asub,row,col);
Bs[row][col] = GetElement(Bsub,row,col);
// 保證整個子矩陣都被讀入
// 也就是一個block的所有thread都執行到這里
__syncthreads();
for(int e=0; e<BLOCK_SIZE;++e){
Cvalue += As[row][e] * Bs[e][col];
}
// 保證這個block算完，然后下一次訓練讀入下一輪的sub matrix
__syncthreads();
}
SetElement(Csub,row,col,Cvalue)
}
假設A矩陣是m行n列的，B矩陣是p行q列的，如果用前一種方法，A和B矩陣從global memory中被讀取的次數是np次，后一種方法是(n/block_size) * (p/block_size)次。(書上那種說法應該指的是同一行、列數據被讀取的次數，跟我這里不同)。
大致可用下圖表示：
Pinned Memory
Host端的內存分為頁內存(Pageable Memory)和鎖存(Pinned Memory)，前者在內存不足時會與虛擬內存交換，把數據放在虛擬內存中，后者內存是固定的。
而內存拷貝是一個非常耗時的操作，由於GPU知道內存的物理地址,因此就可以使用DMA技術來在GPU和CPU之間復制數據.當使用可分頁的內存進行復制時(使用malloc),CUDA驅動程序仍會通過dram把數據傳給GPU,這時復制操作會執行兩遍,第一遍從可分頁內存復制一塊到臨時的頁鎖定內存,第二遍是再從這個頁鎖定內存復制到GPU上.當從可分頁內存中執行復制時,復制速度將受限制於PCIE總線的傳輸速度和系統前段速度相對較低的一方.在某些系統中,這些總線在帶寬上有着巨大的差異,因此當在GPU和主機之間復制數據時,這種差異會使頁鎖定主機內存比標准可分頁的性能要高大約2倍.即使PCIE的速度於前端總線的速度相等,由於可分頁內訓需要更多一次的CPU參與復制操作,也會帶來額外的開銷.
所以在內存足夠的時候，可以使用pinned memory來加速數據讀取。