文章目錄
- 1. 為什麼需要使用GPU
- 2. GPU為什麼性能高
- 3. 如何運用GPU進行程式設計
-
- 3.1 NVIDIA GPU Architecture
- 3.2 Thread Hierarchy
- 3.3 Execution Model
- 3.4 kernel function
- 4. An example: Matrix Multiplication
1. 為什麼需要使用GPU
為什麼GPU（Graphics Processing Unit）程式設計越來越流行，主要是因為GPU相對于CPU的運算速度，記憶體帶寬均有較大的優勢，下面是摘自《CUDA C PROGRAMMING GUIDE》中的圖檔：
浮點數運算速度：
記憶體帶寬：
2. GPU為什麼性能高
這是因為GPU中硬體更多的用于data processing而不是data caching 或 flow control
NVIDIA GPU 更是采用了SIMT (Single-Instruction, Multiple-Thread)和Hardware Multithreading 技術來進行計算加速：
-
SIMT 相對于SIMD(Single Instruction, Multiple Data),前者主要采用線程并行的方式，後者主要采用資料并行的方式。
下面是一個采用SIMD進行運算的例子：
下面是一個SIMT的例子：void add(uint32_t *a, uint32_t *b, uint32_t *c, int n) { for(int i=0; i<n; i+=4) { //compute c[i], c[i+1], c[i+2], c[i+3] uint32x4_t a4 = vld1q_u32(a+i); uint32x4_t b4 = vld1q_u32(b+i); uint32x4_t c4 = vaddq_u32(a4,b4); vst1q_u32(c+i,c4); } }__global__ void add(float *a, float *b, float *c) { int i = blockIdx.x * blockDim.x + threadIdx.x; a[i]=b[i]+c[i]; //no loop! } - Hardware Multithreading技術主要是将程序的運作上下文一直儲存在硬體上，因而不存在運作上下文切換帶來開銷的問題（傳統的CPU多程序是将程序運作上下文儲存在記憶體中，程序切換時涉及到記憶體的讀取，因而開銷較大）
3. 如何運用GPU進行程式設計
既然GPU有這麼多的優勢，那麼如何使用GPU進行程式設計呢？由于GPU種類很多，不同的GPU都有不同的硬體實作以及相應的軟體接口。目前比較流行的是NVIDIA GPU, 這主要是因為其提供了一套易用的軟體接口CUDA， CUDA（Compute Unified Device Architecture）是NVIDIA公司基于其生産的圖形處理器GPU開發的一個并行計算平台和程式設計模型。
3.1 NVIDIA GPU Architecture
NVIDIA GPU的硬體架構一般如下，以GeForce8600 為例：
每個GPU中都有多個多流處理器Streaming Multiprocessors（簡稱SM，有時也直接叫做Multiprocessor), 每個Multiprocessors中有多個core,線程最終就是在這些core上運作的。
這些硬體資訊可以通過CUDA Runtime API 擷取，例如，我的Lenovo T440P上的GPU硬體資訊如下：
CUDA Device Query (Runtime API) version (CUDART static linking)
Detected 1 CUDA Capable device(s)
Device 0: "GeForce GT 730M"
CUDA Driver Version / Runtime Version 10.0 / 10.0
CUDA Capability Major/Minor version number: 3.5
Total amount of global memory: 984 MBytes (1031405568 bytes)
( 2) Multiprocessors, (192) CUDA Cores/MP: 384 CUDA Cores
GPU Max Clock rate: 758 MHz (0.76 GHz)
Memory Clock rate: 1001 Mhz
Memory Bus Width: 64-bit
L2 Cache Size: 524288 bytes
Maximum Texture Dimension Size (x,y,z) 1D=(65536), 2D=(65536, 65536), 3D=(4096, 4096, 4096)
Maximum Layered 1D Texture Size, (num) layers 1D=(16384), 2048 layers
Maximum Layered 2D Texture Size, (num) layers 2D=(16384, 16384), 2048 layers
Total amount of constant memory: 65536 bytes
Total amount of shared memory per block: 49152 bytes
Total number of registers available per block: 65536
Warp size: 32
Maximum number of threads per multiprocessor: 2048
Maximum number of threads per block: 1024
Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
Max dimension size of a grid size (x,y,z): (2147483647, 65535, 65535)
Maximum memory pitch: 2147483647 bytes
Texture alignment: 512 bytes
Concurrent copy and kernel execution: Yes with 1 copy engine(s)
Run time limit on kernels: Yes
Integrated GPU sharing Host Memory: No
Support host page-locked memory mapping: Yes
Alignment requirement for Surfaces: Yes
Device has ECC support: Disabled
Device supports Unified Addressing (UVA): Yes
Device supports Compute Preemption: No
Supports Cooperative Kernel Launch: No
Supports MultiDevice Co-op Kernel Launch: No
Device PCI Domain ID / Bus ID / location ID: 0 / 2 / 0
Compute Mode:
< Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >
deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 10.0, CUDA Runtime Version = 10.0, NumDevs = 1
Result = PASS
該GPU有2個Multiprocessor, 每個multiprocessor有192個core，總計384個core. 對于現在的Tesla型号的GPU,其core數為3584（56 * 64), 每個core都有其相對獨立的寄存器等，這是GPU高性能的基礎。
3.2 Thread Hierarchy
在NVIDIA GPU程式設計中，一個多線程的程式會采用分組的方式在GPU上運作，每個組稱為一個block，每個block中含有若幹個線程。每個thread block在一個Multiprocessor上運作；多個thread blocks可以在一個或多個Multiprocessor上運作。這樣做的好處是當增加GPU中Multiprocessor的個數時，程式性能可以随之提高。
Block在Grid中的排列形式可以是1D或2D（沒有3D的block)，每個block中有若幹線程，這些線程在block中的排列方式可以是1D/2D/3D，如下圖：
在GPU程式設計中，相應的概念均可以找到具體的實體實體:
- Grid 對應于GPU，一個GPU就是一個Grid，在多GPU的機器上，将會有多個Grid。
- Block對應（從屬于）MultiProcessors這個實體實體
- Thread對應于MultiProcessors下面的core這個實體實體，thread 運作在core上
具體的，當一個block運作在multiprocessor時，multiprocessor是以wrap為機關來排程block中的線程的，一個wrap一般是32個線程，這也就是我們為什麼說NVIDIA GPU采用SIMT的原因。wrap是來源于實際生活中的概念(織布中用的經，經紗)，下圖中所有的豎線即為一個wrap:
對應于上面硬體GeForce GT 730M，其線程相關參數如下：
- 每個Multiprocessor 最多可支援2048個線程；
- 每個thread block中最多可支援1024個線程；
- 每個thread block中維數方面x,y,z分别最多為1024,1024,64
- 每個grid中維數方面x,y,z分别最多為2147483647, 65535, 65535
Maximum number of threads per multiprocessor: 2048
Maximum number of threads per block: 1024
Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
Max dimension size of a grid size (x,y,z): (2147483647, 65535, 65535)
注意這裡grid size的z方向雖然最大可以是65535，但是在CUDA的實際程式設計接口中隻能是1.
3.3 Execution Model
采用CUDA程式設計時，程式的運作步驟一般如下：
1.準備GPU計算資料： 将資料從host記憶體拷貝到GPU記憶體
2.在GPU中運作程式
3.将計算結果從GPU記憶體拷貝到CPU記憶體
在GPU和CPU混合程式設計中，通常将GPU叫做device, 将CPU叫做host。如上步驟2中能夠在Host端被調用，在device端執行的函數叫kernel function。
3.4 kernel function
對于運作在device端的函數，一般以__global__和以__device__ 作為标志。以__device__作為标志的函數隻能在device上被調用；以__global__作為标志的函數可以在host端調用，也可以在device端調用，一般稱為kernel function, 調用kernel function時我們需要提供兩個參數：
- 以block為機關的，在grid内部block在x,y方向（不支援z方向）的維數B
- 以thread為機關的，在block内線程在x,y,z方向的維數T
kernel function調用的一般形式為：
myKernel<<< B, T >>>(arg1, … );
B,T在CUDA中采用如下類似的資料結構dim3:
struct dim3 {x; y; z;};
其提供了int到dim3的隐式類型轉換：
myKernel<<< 2, 3 >>>(arg1, … );
上面的參數等價于dim3 b(2,1,1) T(3,1,1)。CUDA為所有在device内運作的function提供了如下兩個内置變量gridDim和blockDim:
dim3 gridDim
dim3 blockDim
- 通過gridDim.x,gridDim.y,gridDim.z,擷取grid在x,y,z方向的維數，也就是block在grid内部x,y,z方向的個數，gridDim.z始終為1
- 通過blockDim.x,blockDim.y,blockDim.z,擷取block在x,y,z方向的維數，也就是線程在block内部x,y,z方向的個數
那麼程式中使用到的block數和單個block内部線程總數将分别是：
gridDim.x * gridDim.y*gridDim.z // girdDim.z = 1
blockDim.x * blockDim.y * blockDim.z
對于kernel function的調用，采用的是SIMT的方式，也就是說同一個function的函數指令将會運作在多個線程中，而線程又屬于某個block，我們怎麼擷取這些線程的索引(index)呢? CUDA 提供了兩個可以在kernel function内部使用的變量：
uint3 blockIdx
uint3 threadIdx
- 通過blockIdx.x, blockIdx.y擷取到目前block在grid内部x,y方向的索引
- 通過threadIdx.x, threadIdx.y, threadIdx.z擷取thread在block内部x,y,z方向的索引
對于2D Grid和2D block，線程在x,y方向的全局唯一ID就可以通過如下計算得到：
- x = blockIdx.x * blockDim.x + threadIdx.x;
- y = blockIdx.y * blockDim.y + threadIdx.y;
下面是一個2D Grid和2D block的示意圖，：
對于2D Grid和3D block的情形，有類似：
- x = blockIdx.x * blockDim.x + threadIdx.x;
- y = blockIdx.y * blockDim.y + threadIdx.y;
- z = blockIdx.z * blockDim.z + threadIdx.z;
注意前面提到過Grid的排列形式沒有3D的，隻有2D的，也就是說blockIdx.z = 0；
4. An example: Matrix Multiplication
下面通過矩陣相乘的例子來說明采用如何使用GPU進行程式設計，回憶一下，對于矩陣A,B,矩陣向乘的結果C中的元素是通過如下公式得到：
具體計算過程如下：
在C中，一般的實作如下：
void matrixMult (int a[N][N], int b[N][N], int c[N][N], int width)
{
for (int i = 0; i < width; i++) {
for (int j = 0; j < width; j++) {
int sum = 0;
for (int k = 0; k < width; k++) {
int m = a[i][k];
int n = b[k][j];
sum += m * n;
}
c[i][j] = sum;
}
}
}
其中，矩陣width是矩陣A的列數，顯然，上面算法的複雜度是O(N^3)。采用GPU程式設計隻需将上面的方法寫成kernel function的形式：
__global__ void matrixMult (int *a, int *b, int *c, int width) {
int k, sum = 0;
int col = threadIdx.x + blockDim.x * blockIdx.x;
int row = threadIdx.y + blockDim.y * blockIdx.y;
if(col < width && row < width) {
for (k = 0; k < width; k++) ｛
sum += a[row * width + k] * b[k * width + col];
｝
c[row * width + col] = sum;
}
}
對比一下C和GPU實作的線程數量和時間複雜度:
| 線程數量 | 時間複雜度 | |
|---|---|---|
| C | 1 | N^3 |
| GPU | N^2 | N |
較完整的GPU實作代碼如下：
#define N 16
#include <stdio.h>
__global__ void matrixMult (int *a, int *b, int *c, int width) {
int col = threadIdx.x + blockDim.x * blockIdx.x;
int row = threadIdx.y + blockDim.y * blockIdx.y;
if(col < width && row < width) {
for (k = 0; k < width; k++) {
sum += a[row * width + k] * b[k * width + col];
}
c[row * width + col] = sum;
}
int main() {
int a[N][N], b[N][N], c[N][N];
int *dev_a, *dev_b, *dev_c;
// initialize matrices a and b with appropriate values
int size = N * N * sizeof(int);
cudaMalloc((void **) &dev_a, size);
cudaMalloc((void **) &dev_b, size);
cudaMalloc((void **) &dev_c, size);
cudaMemcpy(dev_a, a, size, cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, size, cudaMemcpyHostToDevice);
dim3 dimGrid(1, 1);
dim3 dimBlock(N, N);
matrixMult<<<dimGrid, dimBlock>>>(dev_a, dev_b, dev_c, N);
cudaMemcpy(c, dev_c, size, cudaMemcpyDeviceToHost);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
}