GPU程式設計(二): GPU架構瞭解一下!
目錄
- 前言
- GPU架構
- GPU處理單元
- 概念GPU
- GPU執行緒與SM
- GPU執行緒
- SM
- 加法
- 統一記憶體
- 乘法
- 最後
前言
在實際CUDA程式設計之前, 先來了解下GPU的結構. 和CPU相比顯得粗暴又強大(手動滑稽).
GPU架構
GPU處理單元
從這張GPU概念核心圖開始講起, 會發現和CPU核心是不同的, 少了三級快取, 分支預測等等. 但是增加了ALU的數量, 擴大了上下文儲存池(Pool of context storge).
可以看到, 上下文儲存池分成4份, 也就是說, 可以執行4條指令流, 比方說指令1阻塞, 立馬切換指令2, 指令2阻塞切換指令3, 這就起到了隱藏延遲的效果. 當然數量到底是多少是很講究的, 不是越多越好.
總的來看, 核心含8個ALU, 4組執行環境(Execution context), 每組有8個Ctx. 這樣, 一個這樣的核心可以併發(concurrent but interleaved)執行4條指令流(instruction streams), 32個併發程式片元(fragment).
概念GPU
複製16個上述的處理單元, 得到一個GPU. 實際肯定沒有這麼簡單的, 所以說是概念GPU.
這個GPU含16個處理單元, 128個ALU, 64組執行環境(Execution context), 512個併發程式片元(fragment).
祭出n多年前的卡皇GTX 480, 有480個CUDA核(也就是ALU), 記憶體頻寬177.4GB/s. 而GTX 980 Ti有2816個CUDA核, 記憶體頻寬336.5GB/s.
但是頻寬依舊是瓶頸, 雖然比CPU頻寬高了一個數量級, 但是可以看到, GTX 980 Ti的頻寬也就是多年前GTX 480的兩倍左右.
GPU執行緒與SM
由於目前還沒有完全依靠GPU執行得機器, 一般來說, 都是異構的, CPU+GPU. 這一點是要特別注意的, 也就是Host與Device.
GPU執行緒
在CUDA架構下, 顯示晶片執行時的最小單位是thread. 數個thread可以組成一個block. 一個block中的thread能存取同一塊共享的記憶體, 而且可以快速進行同步的動作. 不同block中的thread無法存取同一個共享的記憶體, 因此無法直接互通或進行同步. 因此, 不同block中的thread能合作的程度是比較低的. 上圖:
然後依據thread, block和grid, 有著不同的儲存. 核心就是thread. 可以結合下圖進行理解:
- 每個處理器上有一組本地32位暫存器(Registers);
- 並行資料快取或共享儲存器(Shared Memory), 由所有標量處理器核心共享, 共享儲存器空間就位於此處;
- 只讀固定快取(Constant Cache), 由所有標量處理器核心共享, 可加速從固定儲存器空間進行的讀取操作(這是裝置儲存器的一個只讀區域);
- 一個只讀紋理快取(Texture Cache), 由所有標量處理器核心共享, 加速從紋理儲存器空間進行的讀取操作(這是裝置儲存器的一個只讀區域), 每個多處理器都會通過實現不同定址模型和資料過濾的紋理單元訪問紋理快取.
SM
看到上圖, GPU硬體的一個核心元件是SM, SM是英文名是Streaming Multiprocessor, 翻譯過來就是流式多處理器. SM的核心元件包括CUDA核心, (其實就是ALU, 如上圖綠色小塊就是一個CUDA核心)共享記憶體, 暫存器等, SM可以併發地執行數百個執行緒, 併發能力就取決於SM所擁有的資源數. 當一個kernel被執行時, 它的gird中的執行緒塊被分配到SM上, 一個執行緒塊只能在一個SM上被排程. SM一般可以排程多個執行緒塊, 這要看SM本身的能力. 那麼有可能一個kernel的各個執行緒塊被分配多個SM, 所以grid只是邏輯層, 而SM才是執行的物理層.
下圖是我GT 750M的顯示卡資訊:
SM採用的是SIMT(Single-Instruction, Multiple-Thread, 單指令多執行緒)架構, 基本的執行單元是執行緒束(wraps), 執行緒束包含32個執行緒, 這些執行緒同時執行相同的指令, 但是每個執行緒都包含自己的指令地址計數器和暫存器狀態,也有自己獨立的執行路徑.
加法
試著用CUDA程式設計做一個矩陣加法:
#include <stdio.h>
__global__ void add(float * x, float *y, float * z, int n){
int index = threadIdx.x + blockIdx.x * blockDim.x;
int stride = blockDim.x * gridDim.x;
for (int i = index; i < n; i += stride){
z[i] = x[i] + y[i];
}
}
int main(){
int N = 1 << 20;
int nBytes = N * sizeof (float);
float *x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
for (int i = 0; i < N; i++){
x[i] = 10.0;
y[i] = 20.0;
}
float *d_x, *d_y, *d_z;
cudaMalloc((void**)&d_x, nBytes);
cudaMalloc((void**)&d_y, nBytes);
cudaMalloc((void**)&d_z, nBytes);
cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
dim3 blockSize(256);
// 4096
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
add << < gridSize, blockSize >> >(d_x, d_y, d_z, N);
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyDeviceToHost);
float maxError = 0.0;
for (int i = 0; i < N; i++){
maxError = fmax(maxError, (float)(fabs(z[i] - 30.0)));
}
printf ("max default: %.4f\n", maxError);
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
free(x);
free(y);
free(z);
return 0;
}
由於我是用mac ssh訪問linux主機的, 所以看中文都是亂碼的, 就沒打註釋. 簡單說下.
- 邏輯部分:
申請1M的float, 放入10.0. 申請1M的float, 放入20.0, 然後加起來. 但是我們不存在直接看結果的. 就迴圈計算誤差值, 輸出最大的那個誤差值. 最後看到是0就代表全部計算正確了.- CUDA部分:
cudaMalloc((void**)&d_x, nBytes);是很搶眼的, 意思也很簡單, 在GPU中申請空間, 而不是CPU.
用cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);將CPU中的資料放入到GPU, 注意第二個是源資料, 第三個是方向.
dim3 blockSize(256);猜猜也知道了, 就是申請256個block.dim3 gridSize()同理.
最後cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyDeviceToHost);從GPU中把結果拷貝回CPU, 注意第三個引數和之前的不同.
記得釋放申請的空間.
統一記憶體
時不我待. 在高版本的CUDA中, 升級了申請空間的操作. CUDA 6.x引入統一記憶體(Unified Memory). 具體內容建議查閱我給出的連結, 說的非常細緻. 簡單來說, 就是申請一次就好, 不用先CPU後GPU, 再拷貝來拷貝去, 太傻.
#include <stdio.h>
__global__ void add(float * x, float *y, float * z, int n){
// 同之前, 略
}
int main()
{
int N = 1 << 20;
int nBytes = N * sizeof(float);
float *x, *y, *z;
cudaMallocManaged((void**)&x, nBytes);
cudaMallocManaged((void**)&y, nBytes);
cudaMallocManaged((void**)&z, nBytes);
for (int i = 0; i < N; ++i)
{
x[i] = 10.0;
y[i] = 20.0;
}
dim3 blockSize(256);
// 4096
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
add << < gridSize, blockSize >> >(x, y, z, N);
cudaDeviceSynchronize();
float maxError = 0.0;
for (int i = 0; i < N; i++){
maxError = fmax(maxError, (float)(fabs(z[i] - 30.0)));
}
printf ("max default: %.4f\n", maxError);
cudaFree(x);
cudaFree(y);
cudaFree(z);
return 0;
}
之前看起來有多蠢, 現在就有多簡潔(手動滑稽). 注意到
cudaMallocManaged((void**)&x, nBytes);, 這樣的申請就無需再拷貝操作了.
cudaDeviceSynchronize();同步device, 保證結果能正確訪問. 具體原理請參看之前連結.
乘法
學過線性代數的都知道矩陣乘法是那種不難單頭疼的工作, 原來可能用MATLAB偷懶, 現在可以考慮CUDA了(手動滑稽).
先定義一個Matrix:
struct Matrix
{
int width;
int height;
float *elements;
};
定義矩陣操作函式和乘法函式:
__device__ float getElement(Matrix *A, int row, int col)
{
return A->elements[row * A->width + col];
}
__device__ void setElement(Matrix *A, int row, int col, float value)
{
A->elements[row * A->width + col] = value;
}
__global__ void matMulKernel(Matrix *A, Matrix *B, Matrix *C)
{
float Cvalue = 0.0;
int row = threadIdx.y + blockIdx.y * blockDim.y;
int col = threadIdx.x + blockIdx.x * blockDim.x;
for (int i = 0; i < A->width; ++i)
{
Cvalue += getElement(A, row, i) * getElement(B, i, col);
}
setElement(C, row, col, Cvalue);
}
主函式並不難寫, 主要是理解排程原理:
int main()
{
int width = 1 << 10;
int height = 1 << 10;
Matrix *A, *B, *C;
cudaMallocManaged((void**)&A, sizeof(Matrix));
cudaMallocManaged((void**)&B, sizeof(Matrix));
cudaMallocManaged((void**)&C, sizeof(Matrix));
int nBytes = width * height * sizeof(float);
cudaMallocManaged((void**)&A->elements, nBytes);
cudaMallocManaged((void**)&B->elements, nBytes);
cudaMallocManaged((void**)&C->elements, nBytes);
A->height = height;
A->width = width;
B->height = height;
B->width = width;
C->height = height;
C->width = width;
for (int i = 0; i < width * height; ++i)
{
A->elements[i] = 1.0;
B->elements[i] = 2.0;
}
dim3 blockSize(32, 32);
dim3 gridSize((width + blockSize.x - 1) / blockSize.x,
(height + blockSize.y - 1) / blockSize.y);
matMulKernel << < gridSize, blockSize >> >(A, B, C);
cudaDeviceSynchronize();
float maxError = 0.0;
for (int i = 0; i < width * height; ++i)
maxError = fmax(maxError, fabs(C->elements[i] - 2 * width));
printf ("max fault: %f\n", maxError);
return 0;
}
和加法差不多, 就是多些矩陣操作. 不多說了~
然後矩陣計算程式碼參考這篇文章.
最後
CUDA近期就先更到這裡, 下次一更新應該是核心必須懂部分的. 喜歡記得點贊哦, 有意見或者建議評論區見~