我今天分享的話題是：
CUDA memory & variables - different memory and variable types 的介紹 本文來介紹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
上一篇CUDA系列學習（一）An Introduction to GPU and CUDA中我們提到了memory由host memory和device memory組成，每部分尤其自己獨立的內存空間。Kernel跑在device memory上，所以runtime提供了分配，釋放，複製 device memory 和device <-->host 間transfer data的函數。
5.1 global arrays
global arrays:
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 <常用>
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)
這時每個thread有32個register
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內有以下資源：
這些決定了一個SMX上能同時運行多少個blocks（最多16個）。
參考：
1. CUDA C Programming Guide
2. different memory and variable types
3. CUDA 安裝與配置
4. CUDA調試工具——CUDA GDB
5. GPU工作方式
6. Fermi 架構白皮書（GPU繼承了Fermi的很多架構特點）
7. GTX460架構
以上就是這次的分享，謝謝大家！