1. 程式人生 > >CUDA系列學習(二)CUDA memory variables

CUDA系列學習(二)CUDA memory variables

               

本文來介紹CUDA的memory和變數存放,分為以下章節:

(一)、CPU Memory 結構

(二)、GPU Memory結構

(三)、CUDA Context

(四)、kernel設計

(五)、變數 & Memory

             5.1 global arrays

             5.2 global variables

             5.3 Constant variables

             5.4 Register

             5.5 Local Array

             5.6 Shared Memory

             5.7 Texture Memory

             5.8 總結        

(一)、CPU Memory 結構

CPU提速主要依靠區域性性原理,即時間區域性性和空間區域性性。我們先看一下CPU的記憶體結構:


Data Access

先複習一下資料在這幾級儲存中的傳輸。作為資料transfer的基本單位,cache line的典型大小為8*8(8個變數,每個8bytes)=64bytes. 當一個cache想要load資料到暫存器時,檢查cache中的line,如果hit了就get到資料,否則將整條line從主存中去出來,(通常通過LRU)替換cache中一條line。暫存器傳資料到cache也一樣的過程。

Importance of Locality


上圖中可見在CPU中memory<--->L3 Cache傳輸頻寬為20GB/s, 除以64bytes/line得到傳輸記錄速度約300M line/s,約為300M*8= 2.4G double/s. 一般地,浮點數操作需要兩個輸入+1個輸出,那麼loading 3個數(3 lines)的代價為 100Mflops。如果一個line中的全部8個variables都被用到,那麼每秒浮點操作可以達到800Mflops。而CPU工作站典型為10 Gflops。這就要靠時間區域性性來重用資料了。

(二)、GPU Memory結構


Data Access

  • Kepler GPU的cache line通常為128bytes(32個float or 16個double)。
  • 資料傳輸頻寬最高250GB/s
  • SMX的L2 cache統一1.5MB,L1 cache / shared memory有64KB
  • 沒有CPU中的全域性快取一致性,所以幾乎沒有兩塊block更新相同的全域性陣列元素。

Importance of Locality

GPU對浮點數的操作速度可達1Tflops。和上面CPU的計算類似,GPU中memory<--->L2Cache傳輸頻寬為250GB/s, 除以128bytes/line得到傳輸記錄速度約2G line/s,約為2G*16= 32G double/s. 一般地,浮點數操作需要兩個輸入+1個輸出,那麼loading 3個數(3 lines)的代價為 670Mflops。如果一個line中的全部16個variables都被用到,那麼每秒浮點操作可以達到11Gflops。

這樣的話每進行一次資料到device的傳輸需要45flops(45次浮點操作)才能達到500Gflops. 所以很多演算法基本上不是卡在計算瓶頸,而是傳輸頻寬

(三)、CUDA Context

一個CUDA Context類似於一個CPU程序。程式在Initialization的時候,runtime給每個device建立一個CUDA context,這個context在所有host threads中共享。driver API中的所有資源和action都封裝在一個CUDA context中,context被銷燬的時候系統自動清空這些資源,每個context擁有其自己的地址空間。所以,CUdeviceptr的value在不同context中會指向不同的記憶體空間。

一個host thread同一時刻只能用一個device context,每個host thread都有一個儲存當前contexts的stack。當一個context被cuCtxCreate()建立時,這個新的context被壓入棧(在棧頂),呼叫cuCtxPopCurrent() 可將這個context彈出來,然後這個context就會“漂”到其他host thread中再被壓入棧。

每個context都會維護一個count,表示有多少個threads在用。cuDtrCreate()令count = 1, cuCtxAttach()令count++,cuCtxDetach()令count--,cuCtxDestroy()令count = 0;一旦count=0,這個context就被銷燬。


(四)、kernel設計

我們在CUDA系列學習(一)中提到了GPU用的是SIMT cores,現在看一下它是如何進行執行緒管理的。每個SMX 多處理器在建立,管理,排程,執行的時候將threads每32個組成一組,稱為“wraps”。具體地,一個多處理器分配到多個blocks去執行的時候,它將blocks中的threads 分成wraps而且每個warp被一個warp scheduler來排程執行。一個warp一次執行一條相同指令,所以warp中所有threads同步執行是最有效的。那麼如果warp中的部分threads走上了資料相關的條件分支,warp就連續在各個branch上執行,暫停沒進入branch的threads。直到所有branch上的threads都執行完再合併了一起向下走。所以實現效能提升要注意儘量使warp內執行緒不要出現divergence。另外,注意這個branch divergence 之發生在warp內部;不同warp之間是獨立執行的。

看兩個kernel設計:

__global__ void kernel_1(float* x)int tid = threadIdx.x + blockDim.x * blockIdx.x; x[tid] = threadIdx.x;}__global__ void kernel_2(float* x)int tid = threadIdx.x + blockDim.x * blockIdx.x; x[1000*tid] = threadIdx.x;}

kernel_1中一個warp的32個thread訪問x的相鄰元素,即x[0]~x[31]在相同的cache line, 就是一個好的transfer;

kernel_2中訪問不連續記憶體,就要請求不同cache line,嚴重影響performance

(五)、變數 & Memory

 上一篇中我們提到了memory由host memory和device memory組成,每部分尤其自己獨立的記憶體空間。Kernel跑在device memory上,所以runtime提供了分配,釋放,複製 device memory 和device <-->host 間transfer data的函式。

5.1 global arrays

global arrays:

  • 儲存在/佔用device memory
  • 由host code(非kernel部分code)宣告
  • 一直存在,直到被host code釋放
  • 因為所有block執行順序不定,所以如果一個block修改了一個數組元素,其他block就不能再對該元素進行讀寫

5.2 global variables

宣告前加識別符號__device__,表示變數要放在device上了           e.g.  __device__ int reduction_lock=0;

__shared__(見4.6)和__constant__(見4.3)中至多有一個跟在__device__後面同時使用,標明用哪塊memory空間,如果這兩個都沒寫,則:

  • 變數可以被grid內的所有threads讀寫
  • 與application同生死
  • 也可以定義為array,但是必須指定size
  • 可以在host code中通過以下函式讀寫:

            1. cudaMemcpyToSymbol;     

            2. cudaMemcpyFromSymbol;

            3. cudaMemcpy + cudaGetSymbolAddress

Demo Code:

// float scalar__device__ float devData;float value = 3.14f;cudaMemcpyToSymbol(devData, &value, sizeof(float));//cudaMemcpyToSymbol(const char* symbol, const void* src, size_t count, size_t offset = 0, enum cudaMemcpyKind)// float array__device__ float* devPointer;float* ptr;cudaMalloc(&ptr, 256 * sizeof(float));cudaMemcpyToSymbol(devPointer, &ptr, sizeof(ptr));


5.3 Constant variables <常用>

  • 哪裡宣告隨便,宣告前加識別符號__constant__
  • 與application同生死
  • grid內所有thread可直接讀(不可update),在host code中通過以下函式初始化

            1. cudaMemcpyToSymbol;     

            2. cudaMemcpyFromSymbol;

            3. cudaMemcpy + cudaGetSymbolAddress

Demo Code:

__constant__ float constData[256];float data[256];cudaMemcpyToSymbol(constData, data, sizeof(data)); //cudaMemcpyToSymbol(const char* symbol, const void* src, size_t count, size_t offset = 0, enum cudaMemcpyKind)cudaMemcpyFromSymbol(data, constData, sizeof(data)); //cudaMemcpyFromSymbol(const char* dst, const void* src_symbol, size_t count, size_t offset = 0, enum cudaMemcpyKind)


5.4 Register

  • 預設一個kernel中的所有內部變數都存在register中
  • 64K 32-bit registers per SMX
  • up to 63 registers per thread (up to 255 for K20 / K40)

        這時有64K/63 = 1024個threads (256個threads for K20 / K40)

  • up to 2048 threads (at most 1024 per thread block)

        這時每個thread有32個register

  • not much difference between “fat” and “thin” threads
  • 如果程式需要更多的register呢?就“spill over”到L1 cache,這樣訪問速度就慢了,我們要儘量避免spill

5.5 Local Array

指kernel code中宣告的陣列。

  • 簡單情況下,編譯器會將小陣列float a[3]轉換成3個標量registers:a0,a1,a2作處理
  • 複雜的情況,會將array放到L1(16KB),只能放4096個32-bit的變數,如果有1024個執行緒,每個執行緒只能分配放4個變數。

5.6 Shared Memory

前面加識別符號__shared__  e.g.    __shared__   int  x_dim;  

  • 要佔用thread block的shared memory space.
  • 要比global memory快很多,所以只要有機會就把global memory整成shared memory
  • 與block同生死
  • thread block內所有threads共用(可讀可寫)
  • 啥時侯用呢?當所有threads訪問都是同一個值的時候,這樣就避免用register了

但是有問題就是,如果一個thread block有多個warp(上一篇blog中提到的概念,block中的thread每32個被分到一個warp,最後一個不足32個thread也沒關係,同樣形成一個warp),各warp執行指令順序是不定的,那麼久需要執行緒同步機制,用指令__syncthreads(); 插入一個“barrier”,所有wrap執行到這個barrier之前沒有thread/warp能夠越過去。

Kepler GPU給L1 Cache + shared memory總共64KB,可以分為16+48,32+32,48+16;這個split可以通過cudaFuncSetCacheConfig()或cudaDeviceSetCacheConfig()設定,預設給shared memroy 48KB。這個具體情況看程式了。

下面通過一個經典例子來看shared memory作用:矩陣乘法

目的:實現C=A*B,方法:c[i,j] = A[i,:] * B[:,j],

其中矩陣用row-major表示,即c[i,j] = *(c.elements + i*c.width + j)

1. 不用shared memory優化版:

設A為m*t的矩陣;B為t*n的矩陣;

每個執行緒讀取A的一行,B的一列,計算C的對應值;

所以這樣需要從global memory中讀n次A,m次B。

// Matrices are stored in row-major order:// M(row, col) = *(M.elements + row * M.width + col)typedef struct { int width; int height; float* elements;} Matrix;// Thread block size#define BLOCK_SIZE 16// Forward declaration of the matrix multiplication kernel__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);// Matrix multiplication - Host code// Matrix dimensions are assumed to be multiples of BLOCK_SIZEvoid MatMul(const Matrix A, const Matrix B, Matrix C)// Load A and B to device memory Matrix d_A; d_A.width = A.width; d_A.height = A.height; size_t size = A.width * 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 = B.width * B.height * sizeof(float); cudaMalloc(&d_B.elements, size); cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice); // Allocate C in device memory Matrix d_C; d_C.width = C.width; d_C.height = C.height; size = C.width * C.height * sizeof(float); cudaMalloc(&d_C.elements, size); // Invoke kernel dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE)dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y); MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C); // Read C from device memory cudaMemcpy(C.elements, Cd.elements, size, cudaMemcpyDeviceToHost); } // Free device memory cudaFree(d_A.elements); cudaFree(d_B.elements); cudaFree(d_C.elements);}// Matrix multiplication kernel called by MatMul()__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)// Each thread computes one element of C // by accumulating results into Cvalue float Cvalue = 0int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; for (int e = 0; e < A.width; ++e) Cvalue += A.elements[row * A.width + e]* B.elements[e * B.width + col]; C.elements[row * C.width + col] = Cvalue;}


2. 利用shared memory

每個thread block負責計算一個子矩陣Csub, 其中每個thread負責計算Csub中的一個元素。如下圖所示。為了將fit裝置資源,A,B都分割成很多block_size維的方形matrix,Csub將這些方形matrix的乘積求和而得。每次計算一個乘積時,先將兩個對應方形矩陣從global memory 載入 shared memory(一個thread負責載入A, B兩個sub matrix的元素),然後每個thread計算乘積的一個元素,再由每個thread將這些product加和,存入一個register,最後一次性寫入global memory。計算時注意同步,詳見程式碼。

設A為m*t的矩陣;B為t*n的矩陣;

這樣呢,A只從global memory讀了n/block_size次,B只讀了m/block_size次;


Kernel Code:

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)// Block row and column int blockRow = blockIdx.y; int blockCol = blockIdx.x; // Each thread block computes one sub-matrix Csub of C Matrix Csub = GetSubMatrix(C, blockRow, blockCol); // Each thread computes one element of Csub by accumulating results into Cvalue float Cvalue = 0// Thread row and column within Csub int row = threadIdx.y; int col = threadIdx.x; // Loop over all the sub-matrices of A and B that are // required to compute Csub // Multiply each pair of sub-matrices together // and accumulate the results for (int m = 0; m < (A.width / BLOCK_SIZE); ++m) {  // Get sub-matrix Asub of A  Matrix Asub = GetSubMatrix(A, blockRow, m);  // Get sub-matrix Bsub of B  Matrix Bsub = GetSubMatrix(B, m, blockCol);  // Shared memory used to store Asub and Bsub respectively  __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];  __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];  // Load Asub and Bsub from device memory to shared memory  // Each thread loads one element of each sub-matrix  As[row][col] = GetElement(Asub, row, col);  Bs[row][col] = GetElement(Bsub, row, col);  // Synchronize to make sure the sub-matrices are loaded  // before starting the computation  __syncthreads();  // Multiply Asub and Bsub together  for (int e = 0; e < BLOCK_SIZE; ++e)   Cvalue += As[row][e] * Bs[e][col];  // Synchronize to make sure that the preceding  // computation is done before loading two new  // sub-matrices of A and B in the next iteration  __syncthreads(); } // Write Csub to device memory // Each thread writes one element SetElement(Csub, row, col, Cvalue);}


Host Code:

// Invoke kerneldim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);


5.7 Texture memory

前面加識別符號const __restrict__, 之所以叫texture是因為之前用texture memory想服務於純graphics的應用。

不同於shared memory,對texture memory, 不同執行緒可以訪問到不同value。K20/K40中texture cache有48KB。

5.8 總結 

綜上,每個block內有以下資源:

  • threads
  • registers (registers per thread * number of threads)
  • shared memory

這些決定了一個SMX上能同時執行多少個blocks(最多16個)。

參考:

6. Fermi 架構白皮書(GPU繼承了Fermi的很多架構特點)