一直很好奇GPU做矩陣運算是怎么并行加速的，今天看了一些粗淺的東西，并總結(jié)整理出來。
version：cuda 8

cuda C 中擴展的一些概念

主要包括函數(shù)聲明、變量聲明、內(nèi)存類型聲明、文理內(nèi)存、原子函數(shù)等，常用的有這么幾個：
參考（http://bbs.csdn.net/topics/390798229，原地址已經(jīng)失效）

• 主機
  將CPU及系統(tǒng)的內(nèi)存（內(nèi)存條）稱為主機。
• 設(shè)備
  將GPU及GPU本身的顯示內(nèi)存稱為設(shè)備。
• 線程(Thread)
  一般通過GPU的一個核進行處理。
• 線程塊(Block)
  1. 由多個線程組成（可以表示成一維，二維，三維，具體下面再細說）。
  2. 各block是并行執(zhí)行的，block間無法通信，也沒有執(zhí)行順序。
  3. 注意線程塊的數(shù)量限制為不超過65535（硬件限制）。

• 線程格(Grid)
  由多個線程塊組成（可以表示成一維，二維，三維，具體下面再細說）。

  
  
• 線程束
  在CUDA架構(gòu)中，線程束是指一個包含32個線程的集合，這個線程集合被“編織在一起”并且“步調(diào)一致”的形式執(zhí)行。在程序中的每一行，線程束中的每個線程都將在不同數(shù)據(jù)上執(zhí)行相同的命令。
• 核函數(shù)（Kernel）
  1. 在GPU上執(zhí)行的函數(shù)通常稱為核函數(shù)。
  2. 一般通過標識符global修飾，調(diào)用通過<<<參數(shù)1,參數(shù)2>>>，用于說明內(nèi)核函數(shù)中的線程數(shù)量，以及線程是如何組織的。
  3. 以線程格（Grid）的形式組織，每個線程格由若干個線程塊（block）組成，而每個線程塊又由若干個線程（thread）組成。
  4. 是以block為單位執(zhí)行的。
  5. 叧能在主機端代碼中調(diào)用。
  6. 調(diào)用時必須聲明內(nèi)核函數(shù)的執(zhí)行參數(shù)。
  7. 在編程時，必須先為kernel函數(shù)中用到的數(shù)組或變量分配好足夠的空間，再調(diào)用kernel函數(shù)，否則在GPU計算時會發(fā)生錯誤，例如越界或報錯，甚至導(dǎo)致藍屏和死機。
  //核函數(shù)聲明，前面的關(guān)鍵字global

__global__ void kernel( void ) {
}

函數(shù)修飾符
1. __global__，表明被修飾的函數(shù)在設(shè)備上執(zhí)行，但在主機上調(diào)用。
2. __device__，表明被修飾的函數(shù)在設(shè)備上執(zhí)行，但只能在其他__device__函數(shù)或者__global__函數(shù)中調(diào)用。

常用的GPU內(nèi)存函數(shù)

• cudaMalloc()
  1. 函數(shù)原型： cudaError_t cudaMalloc (void **devPtr, size_t size)。
  2. 函數(shù)用處：與C語言中的malloc函數(shù)一樣，只是此函數(shù)在GPU的內(nèi)存你分配內(nèi)存。
  3. 注意事項：
  3.1. 可以將cudaMalloc()分配的指針傳遞給在設(shè)備上執(zhí)行的函數(shù)；
  3.2. 可以在設(shè)備代碼中使用cudaMalloc()分配的指針進行設(shè)備內(nèi)存讀寫操作；
  3.3. 可以將cudaMalloc()分配的指針傳遞給在主機上執(zhí)行的函數(shù)；
  3.4. 不可以在主機代碼中使用cudaMalloc()分配的指針進行主機內(nèi)存讀寫操作（即不能進行解引用）。
• cudaMemcpy()
  1. 函數(shù)原型：cudaError_t cudaMemcpy (void *dst, const void *src, size_t count, cudaMemcpyKind kind)。
  2. 函數(shù)作用：與c語言中的memcpy函數(shù)一樣，只是此函數(shù)可以在主機內(nèi)存和GPU內(nèi)存之間互相拷貝數(shù)據(jù)。
  3. 函數(shù)參數(shù)：cudaMemcpyKind kind表示數(shù)據(jù)拷貝方向，如果kind賦值為cudaMemcpyDeviceToHost表示數(shù)據(jù)從設(shè)備內(nèi)存拷貝到主機內(nèi)存。
  4. 與C中的memcpy()一樣，以同步方式執(zhí)行，即當函數(shù)返回時，復(fù)制操作就已經(jīng)完成了，并且在輸出緩沖區(qū)中包含了復(fù)制進去的內(nèi)容。
  5. 相應(yīng)的有個異步方式執(zhí)行的函數(shù)cudaMemcpyAsync()，這個函數(shù)詳解請看下面的流一節(jié)有關(guān)內(nèi)容。
• cudaFree()
  1. 函數(shù)原型：cudaError_t cudaFree ( void* devPtr )。
  2. 函數(shù)作用：與c語言中的free()函數(shù)一樣，只是此函數(shù)釋放的是cudaMalloc()分配的內(nèi)存。
  下面實例用于解釋上面三個函數(shù)

GPU內(nèi)存分類

• 全局內(nèi)存
  通俗意義上的設(shè)備內(nèi)存。
• 共享內(nèi)存
  1. 位置：設(shè)備內(nèi)存。
  2. 形式：關(guān)鍵字__shared__添加到變量聲明中。如__shared__ float cache[10]。
  3. 目的：對于GPU上啟動的每個線程塊，CUDA C編譯器都將創(chuàng)建該共享變量的一個副本。線程塊中的每個線程都共享這塊內(nèi)存，但線程卻無法看到也不能修改其他線程塊的變量副本。這樣使得一個線程塊中的多個線程能夠在計算上通信和協(xié)作。
• 常量內(nèi)存
  1. 位置：設(shè)備內(nèi)存
  2. 形式：關(guān)鍵字__constant__添加到變量聲明中。如__constant__ float s[10];。
  3. 目的：為了提升性能。常量內(nèi)存采取了不同于標準全局內(nèi)存的處理方式。在某些情況下，用常量內(nèi)存替換全局內(nèi)存能有效地減少內(nèi)存帶寬。
  4. 特點：常量內(nèi)存用于保存在核函數(shù)執(zhí)行期間不會發(fā)生變化的數(shù)據(jù)。變量的訪問限制為只讀。NVIDIA硬件提供了64KB的常量內(nèi)存。不再需要cudaMalloc()或者cudaFree(),而是在編譯時，靜態(tài)地分配空間。
  5. 要求：當我們需要拷貝數(shù)據(jù)到常量內(nèi)存中應(yīng)該使用cudaMemcpyToSymbol()，而cudaMemcpy()會復(fù)制到全局內(nèi)存。
  6. 性能提升的原因：
  6.1. 對常量內(nèi)存的單次讀操作可以廣播到其他的“鄰近”線程。這將節(jié)約15次讀取操作。（為什么是15，因為“鄰近”指半個線程束，一個線程束包含32個線程的集合。）
  6.2. 常量內(nèi)存的數(shù)據(jù)將緩存起來，因此對相同地址的連續(xù)讀操作將不會產(chǎn)生額外的內(nèi)存通信量。

• 紋理內(nèi)存
  1. 位置：設(shè)備內(nèi)存
  2. 目的：能夠減少對內(nèi)存的請求并提供高效的內(nèi)存帶寬。是專門為那些在內(nèi)存訪問模式中存在大量空間局部性的圖形應(yīng)用程序設(shè)計，意味著一個線程讀取的位置可能與鄰近線程讀取的位置“非常接近”。如下圖：

  
  
  
  

  3. 紋理變量（引用）必須聲明為文件作用域內(nèi)的全局變量。
  4. 形式：分為一維紋理內(nèi)存 和 二維紋理內(nèi)存。
  4.1. 一維紋理內(nèi)存
  4.1.1. 用texture<類型>類型聲明，如texture<float> texIn。
  4.1.2. 通過cudaBindTexture()綁定到紋理內(nèi)存中。
  4.1.3. 通過tex1Dfetch()來讀取紋理內(nèi)存中的數(shù)據(jù)。
  4.1.4. 通過cudaUnbindTexture()取消綁定紋理內(nèi)存。
  4.2. 二維紋理內(nèi)存
  4.2.1. 用texture<類型,數(shù)字>類型聲明，如texture<float，2> texIn。
  4.2.2. 通過cudaBindTexture2D()綁定到紋理內(nèi)存中。
  4.2.3. 通過tex2D()來讀取紋理內(nèi)存中的數(shù)據(jù)。
  4.2.4. 通過cudaUnbindTexture()取消綁定紋理內(nèi)存。

• 固定內(nèi)存
  1. 位置：主機內(nèi)存。
  2. 概念：也稱為頁鎖定內(nèi)存或者不可分頁內(nèi)存，操作系統(tǒng)將不會對這塊內(nèi)存分頁并交換到磁盤上，從而確保了該內(nèi)存始終駐留在物理內(nèi)存中。因此操作系統(tǒng)能夠安全地使某個應(yīng)用程序訪問該內(nèi)存的物理地址，因為這塊內(nèi)存將不會破壞或者重新定位。
  3. 目的：提高訪問速度。由于GPU知道主機內(nèi)存的物理地址，因此可以通過“直接內(nèi)存訪問DMA（Direct Memory Access)技術(shù)來在GPU和主機之間復(fù)制數(shù)據(jù)。由于DMA在執(zhí)行復(fù)制時無需CPU介入。因此DMA復(fù)制過程中使用固定內(nèi)存是非常重要的。
  4. 缺點：使用固定內(nèi)存，將失去虛擬內(nèi)存的所有功能；系統(tǒng)將更快的耗盡內(nèi)存。
  5. 建議：對cudaMemcpy()函數(shù)調(diào)用中的源內(nèi)存或者目標內(nèi)存，才使用固定內(nèi)存，并且在不再需要使用它們時立即釋放。
  6. 形式：通過cudaHostAlloc()函數(shù)來分配；通過cudaFreeHost()釋放。
  7. 只能以異步方式對固定內(nèi)存進行復(fù)制操作。
• 原子性
  1. 概念：如果操作的執(zhí)行過程不能分解為更小的部分，我們將滿足這種條件限制的操作稱為原子操作。
  2. 形式：函數(shù)調(diào)用，如atomicAdd（addr,y)將生成一個原子的操作序列，這個操作序列包括讀取地址addr處的值，將y增加到這個值，以及將結(jié)果保存回地址addr。
  常用線程操作函數(shù)
  1. 同步方法__syncthreads()，這個函數(shù)的調(diào)用，將確保線程塊中的每個線程都執(zhí)行完__syscthreads()前面的語句后，才會執(zhí)行下一條語句。

cuda C 做矩陣乘法（Tiled 算法）

為什么看cuda C 做矩陣乘法運算呢？在深度神經(jīng)網(wǎng)絡(luò)中，全連接層、卷積層、池化層，幾乎我們可以想到的所有操作都離不開矩陣運算，卷積層最后其實也是轉(zhuǎn)化為矩陣的乘法操作進行優(yōu)化，在【conv2d 實現(xiàn) caffe&tensorflow】中有介紹原理。
參考視頻地址：https://www.youtube.com/watch?v=SqZaletdPCY

思想： 為了引入共享內(nèi)存的概念降低GPU帶寬使用，把要計算的兩個矩陣A B 先分解成BLOCK_SIZE=16大小的submatrix，每一個block結(jié)構(gòu)運算一個submatrix乘法，這樣在一個block中所有的線程是共享參數(shù)的，不用每次計算都從global memory中重新加載。

template <int BLOCK_SIZE> __global__ void
matrixMulCUDA(float *C, float *A, float *B, int wA, int wB)
{
 // Thread 所在 block 的 location 
 int bx = blockIdx.x;
 int by = blockIdx.y;

 // Thread 的location 
 int tx = threadIdx.x;
 int ty = threadIdx.y;

 // A矩陣16 * 16 子矩陣的起始下標
 int aBegin = wA * BLOCK_SIZE * by;

 // A矩陣16 * 16 子矩陣的終止下標（就是A矩陣一次運算一行，對應(yīng)著B 矩陣一次運算一列）
 int aEnd = aBegin + wA - 1;

 // A矩陣下標一次移動的步長， 子矩陣是16 * 16，一次處理一個子矩陣，那么步長顯然就是16了
 int aStep = BLOCK_SIZE;

 // B 矩陣子矩陣對應(yīng)的起始下標
 int bBegin = BLOCK_SIZE * bx;

 // B 矩陣子矩陣對應(yīng)的步長，一次移動16*widthB，同樣也是隔出16*16的子矩陣出來
 int bStep = BLOCK_SIZE * wB;

 // 累加，得到行 * 列的值
 float Csub = 0;

 // 循環(huán)次數(shù)等于widthA / 16，把長向量點積運算轉(zhuǎn)化為兩個短向量點積后的和
 for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep)
 {
 // 定義A的共享子矩陣變量,因為__shared__聲明，所以同一個block中的所有threads都可見，
 //每個thread填充一個元素，并計算一個行列乘積，減小帶寬使用
 __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];

 // 定義A的共享子矩陣變量
 __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

 // 每個block包含16 * 16 個線程，所以每個線程負責一個矩陣元素的拷貝（注意同步）
 As[ty][tx] = A[a + wA * ty + tx];
 Bs[ty][tx] = B[b + wB * ty + tx];

 // Synchronize to make sure the matrices are loaded
 __syncthreads();

 // 每個線程計算 子矩陣的行列乘積，大循環(huán)外邊還有累加，累加的是不同子矩陣點積和
 for (int k = 0; k < BLOCK_SIZE; ++k)
 {
 Csub += As[ty][k] * Bs[k][tx];
 }

 // 再次同步
 __syncthreads();
 }

 // 把結(jié)果賦值到C矩陣，計算結(jié)果對應(yīng)C下邊的過程
 int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
 C[c + wB * ty + tx] = Csub;
}

只看源代碼很難理解矩陣加速真正的原理，這是一個坑，還有是輸入矩陣的尺寸大小，只能是BLOCK_SIZE=16的整數(shù)倍，不然會出錯（實驗結(jié)果也表明確實出錯了，又是一個坑）。
為什么采用Tiled 算法呢？主要是不這么做GPU從global memory讀取數(shù)據(jù)的代價太大了。