本文來介紹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 = 0; int 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的很多架構特點）