前言
CUDA並行程序設計系列是本人在學習CUDA時整理的資料，內容大都來源於對《CUDA並行程序設計：GPU編程指南》、《GPU高性能編程CUDA實戰》和CUDA Toolkit Documentation的整理。通過本系列整體介紹CUDA並行程序設計。內容包括GPU簡介、CUDA簡介、環境搭建、線程模型、內存、原子操作、同步、流和多GPU架構等。
本系列目錄：
- 【CUDA並行程序設計系列（1）】GPU技術簡介
- 【CUDA並行程序設計系列（2）】CUDA簡介及CUDA初步編程
- 【CUDA並行程序設計系列（3）】CUDA線程模型
- 【CUDA並行程序設計系列（4）】CUDA內存
- 【CUDA並行程序設計系列（5）】CUDA原子操作與同步
- 【CUDA並行程序設計系列（6）】CUDA流與多GPU
- 關於CUDA的一些學習資料
本章將介紹CUDA的內存結構，通過實例展示寄存器和共享內存的使用。
CUDA內存結構
GPU的內存結構和CPU類似，但也存在一些區別，GPU的內存中可讀寫的有：寄存器(registers)、Local memory、共享內存（shared memory）和全局內存（global memory），只讀的有：常量內存（constant memory）和紋理內存（texture memory）。
CUDA Toolkit Document給出的的內存結構如下圖所示：
每個線程都有獨立的寄存器和Local memory，同一個block的所有線程共享一個共享內存，全局內存、常量內存和紋理內存是所有線程都可訪問的。全局內存、常量內存和紋理內存對程序的優化有特殊作用。
寄存器
與CPU不同，GPU的每個SM（流多處理器）有成千上萬個寄存器，在GPU技術簡介中已經提到，SM類似於CPU的核，每個SM擁有多個SP（流處理器），所有的工作都是在SP上處理的，GPU的每個SM可能有8~192個SP，這就意味着，SM可同時運行這些數目的線程。
寄存器是每個線程私有的，並且GPU沒有使用寄存器重命名機制，而是致力於為每一個線程都分配真實的寄存器，CUDA上下文切換機制非常高效，幾乎是零開銷。當然，這些細節對程序員是完全透明的。
和CPU一樣，訪問寄存器的速度是非常快的，所以應盡量優先使用寄存器。無論是CPU還是GPU，通過寄存器的優化方式都會使程序的執行速度得到很大提高。
舉一個例子：
for (int i = 0; i < size; ++i)
{
sum += array[i];
}
sum如果存於內存中，則需要做size次讀/寫內存的操作，而如果把sum設置為局部變量，把最終結果寫回內存，編譯器會將其放入寄存器中，這樣只需1次內存寫操作，將大大節約運行時間。
Local memory
Local memory和寄存器類似，也是線程私有的，訪問速度比寄存器稍微慢一點。事實上，是由編譯器在寄存器全部使用完的時候自動分配的。在優化程序的時候可以考慮減少block的線程數量以使每個線程有更多的寄存器可使用，這樣可減少Local memory的使用，從而加快運行速度。
共享內存
共享內存允許同一個block中的線程讀寫這一段內存，但線程無法看到也無法修改其它block的共享內存。共享內存緩沖區駐留在物理GPU上，所以訪問速度也是很快的。事實上，共享內存的速度幾乎在所有的GPU中都一致（而全局內存在低端顯卡的速度只有高端顯卡的1/10），因此，在任何顯卡中，除了使用寄存器，還要更有效地使用共享內存。
共享內存的存在就可使運行線程塊中的多個線程之間相互通信。共享內存的一個應用場景是線程塊中多個線程需要共同操作某一數據。考慮一個矢量點積運算的例子：
（x1, x2, x3, x4 ) * (y1, y2, y3, y4) = x1y1 + x2y2 + x3y3 + x4y4
和矢量加法一樣，矢量點積也可以在GPU上並行計算，每個線程將兩個相應的元素相乘，然后移到下兩個元素，線程每次增加的索引為總線程的數量，下面是實現這一步的代碼：
const int N = 33 * 1024;
const int threadsPerBlock = 256;
__global__ void dot( float *a, float *b, float *c )
{
__shared__ float cache[threadsPerBlock];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0;
while (tid < N)
{
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
cache[cacheIndex] = temp;
}
CUDA C使用__shared__修飾符申明共享內存的變量。在每個線程中分別計算相應元素的乘積之和，並保存在共享內存變量cache對應的索引中，可以看出，如果只有一個block，那么所有線程結束后，對cache求和就是最終結果。當然，實際會有很多個block，所以需要對所有block中的cache求和，由於共享內存在block之間是不能訪問的，所以需要在各個block中分部求和，並把部分和保存在數組中，最后在CPU上求和。block中分部求和代碼如下：
__global__ void dot( float *a, float *b, float *c ) {
__shared__ float cache[threadsPerBlock];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0;
while (tid < N) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
cache[cacheIndex] = temp;
//同步
__syncthreads();
//分部求和
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
}
__syncthreads()是線程同步函數，調用這個函數確保在線程塊中所有的線程都執行完__syncthreads()之前的代碼，在執行后面的代碼，當然，這會損失一定性能。
當執行__syncthreads()之后的代碼，我們就能確定cache已經計算好了，下面只需要對cache求和就可以了，最簡單的就是用一個for循環計算。但是，這相當只有一個線程在起作用，線程塊其它線程都在做無用功，
使用規約運行是一個更好地選擇，即每個線程將cache中的兩個值相加起來，然后結果保存會cache中，規約的思想如下圖所示。
按這種方法，每次將會使數據減少一半，只需執行log2(threadsPerBlock)個步驟后，就能得到cache中所有值的總和。
最后使用如下代碼將結果保存在c中：
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
這是因為只有一個值需要寫入，用一個線程來操作就行了，如果不加if，那么每個線程都將執行一次寫內存操作，浪費大量的運行時間。
最后，只需要在CPU上把c中的值求和就得到了最終結果。下面給出完整代碼：
#include <stdio.h>
const int N = 33 * 1024;
const int threadsPerBlock = 256;
const int blocksPerGrid = (N + threadsPerBlock -1) / threadsPerBlock;
__global__ void dot( float *a, float *b, float *c )
{
__shared__ float cache[threadsPerBlock];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0;
while (tid < N)
{
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
cache[cacheIndex] = temp;
//同步
__syncthreads();
//規約求和
int i = blockDim.x/2;
while (i != 0)
{
if (cacheIndex < i)
{
cache[cacheIndex] += cache[cacheIndex + i];
}
__syncthreads();
i /= 2;
}
if (cacheIndex == 0)
{
c[blockIdx.x] = cache[0];
}
}
int main(int argc, char const *argv[])
{
float *a, *b, *partial_c;
float *dev_a, *dev_b, *dev_partial_c;
a = (float*)malloc( N*sizeof(float) );
b = (float*)malloc( N*sizeof(float) );
partial_c = (float*)malloc( blocksPerGrid*sizeof(float));
cudaMalloc(&dev_a, N*sizeof(float));
cudaMalloc(&dev_b, N*sizeof(float));
cudaMalloc(&dev_partial_c, blocksPerGrid*sizeof(float));
for (int i=0; i < N; ++i)
{
a[i] = i;
b[i] = i * 2;
}
cudaMemcpy(dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice);
dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b, dev_partial_c );
cudaMemcpy(partial_c, dev_partial_c, blocksPerGrid*sizeof(float), cudaMemcpyDeviceToHost);
int c = 0;
for (int i=0; i < blocksPerGrid; ++i)
{
c += partial_c[i];
}
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_partial_c);
free(a);
free(b);
free(partial_c);
return 0;
}
常量內存
常量內存，通過它的名字就可以猜到它是只讀內存。常量內存其實只是全局內存的一種虛擬地址形式，並沒有特殊保留的常量內存塊。內存的大小為64KB。常量內存可以在編譯時申明為常量內存，使用修飾符__constant__申明，也可以在運行時通過主機端定義為只讀內存。常量只是從GPU內存的角度而言的，CPU在運行時可以通過調用cudaCopyToSymbol來改變常量內存中的內容。
全局內存
GPU的全局內存之所以是全局內存，主要是因為GPU與CPU都可以對它進行寫操作，任何設備都可以通過PCI-E總線對其進行訪問。在多GPU系統同，GPU之間可以不通過CPU直接將數據從一塊GPU卡傳輸到另一塊GPU卡上。在調用核函數之前，使用cudaMemcpy函數就是把CPU上的數據傳輸到GPU的全局內存上。
紋理內存
和常量內存一樣，紋理內存也是一種只讀內存，在特定的訪問模式中，紋理內存能夠提升程序的性能並減少內存流量。紋理內存最初是為圖形處理程序而設計，不過同樣也可以用於通用計算。由於紋理內存的使用非常特殊，有時使用紋理內存是費力不討好的事情。因此，對於紋理內存，只有在應用程序真正需要的時候才對其進行了解。主要應該掌握全局內存、共享內存和寄存器的使用。
參考文獻
- 庫克. CUDA並行程序設計. 機械工業出版社, 2014.
- 桑德斯. GPU高性能編程CUDA實戰. 機械工業出版社, 2011.
- CUDA C Programming Guide
- CUDA Toolkit Documentation
- R. Couturier, Ed., Designing Scientific Applications on GPUs, CRC Press, 2013.