CUDA 高效能平行計算入門

【原文：https://blog.csdn.net/cyhbrilliant/article/details/79434090】

(UPDATE IN 2018.3.8) 
1.更新pitch索引操作的描述

• 概述

  • 什麼是CUDA？ 
    CUDA（Compute Unified Device Architecture）是 NVIDIA公司開發的一種計算架構，可以利用NVIDIA系列顯示卡對一些複雜的計算進行並行加速。

  • 為什麼要用CUDA加速？ 
    在科學計算領域所要用到的計算往往不是我們熟知的普通矩陣，而是千維甚至萬維的矩陣，而普通的CPU序列計算往往不能滿足與科學計算所要求的效能。最好的例子就是深度學習，可以說深度學習起源於感知機但正真發展在於計算能力的提高，也就是顯示卡計算的興起。深度學習的計算都是基於矩陣的計算，而普通一個識別手寫數字的神經網路都有上千個隱含節點，故CUDA效能優化是一門重要的技術。 
    

  • 為什麼顯示卡可以加速？ 
    • 首先，顯示卡可以加速最大的原因是其含有上千個CUDA核心，而CPU的核心往往都在各位數，這就決定了顯示卡可以高度並行的對一些高維度矩陣進行計算。CPU的I/O需要數百上千個週期，序列的計算大部分時間都消耗在I/O上，而顯示卡則不然，NVIDIA顯示卡的架構是類似於SIMD架構的SIMT架構，T指的是Threads，也就是單指令多執行緒架構，上一個執行緒在進行運算操作時下一個執行緒就開始I/O操作，類似於指令流水線的形式，當Threads數量足夠多時，就可以遮蔽I/O所帶來的大量開銷，所以本身架構決定了CPU與顯示卡計算方式的區別。
    • CPU就好比是一個身強體壯的大漢一下可以挑100公斤糧食，顯示卡好比是隻能挑起10公斤糧食的小矮人，但有20個。現在有2000公斤糧食，I/O相當於需要1分鐘時間把糧食裝進桶才能挑走（每次只能一個人裝桶），運算時間相當於挑一次到目的地需要20分鐘。CPU每趟需要21分鐘，共20趟，總共需要420分鐘。顯示卡第一趟也需要21分鐘，且只能挑10公斤糧食，但第一趟在剛裝完桶時，第二個小人就開始裝桶了，在第一個小人送完時最後一個小人也裝完了桶，接著第一個小人又可以裝桶了，所以從共需要2000/10=200個輪迴，但每個輪迴只需要1分鐘，加上最後的小人需要20分鐘送完，總計200×1+20=220分鐘，比CPU快近一倍的時間。當運算需求越大時，GPU的優勢越明顯。 
      

• GPU硬體介紹

  說了這麼多接下去就開始正式入題，想要做CUDA程式設計就必須先了解自己電腦顯示卡的硬體，才能針對自己的顯示卡做必要的優化。當然我指的是NVIDIA系列顯示卡，Intel和AMD就不要過來搗亂了…… 
  
  真正進入顯示卡計算時代的當從Tesla架構開始算起，經歷了Fermi，Kapler，到現在的Maxwell，Pascal，一代比一代計算能力強。這裡就有人要問了，計算能力強的顯示卡一定快嗎？ 
  錯了，NVIDIA官方所說的計算能力指的是單個CUDA核心的架構運算能力，幾千個Fermi核心當然可以隨隨便便秒了幾百個的Maxwell核心（同頻率下…）。而且這裡計算能力並沒有把頻寬給考慮進去，往往限制CUDA運算速度的不是峰值效能，而是視訊記憶體頻寬。 
  想要知道自己的硬體明細當然可以用軟體的方式直接獲取，LINUX系統可以下載CUDA-Z，WINDOWS系統可以下載GPU-Z，MAC土豪請自行尋找軟體。 
  我們以CUDA-Z為例 
  

• 在Core選項卡中我們可以看到上圖所示 
  (ps:請無視840m垃圾的效能) 
  Compute Capability指的是計算能力：1.x 2.x 3.x 5.x….等等，計算能力越高不光是單個核心的計算能力強，並且可以在CUDA程式設計的核函式中支援更多的操作。 
  Clock Rate指的是單個核心的執行頻率：即我們每個CUDA核心每秒所能操作的次數。相當於CPU的主頻。 
  Multiprocessors指的是官方所說的SM數，可以當成許多CUDA核心的集合。 
  CUDA Cores指的是共有多少個計算核心或者說是計算單元：在Kapler架構中，一個SM擁有192個計算單元，而在Maxwell架構中，每個SM只有128個計算單元。所以可以看到840m共有3個SM，有3*128=384個Cores。 
  Warp Size指的是GPU在排程運算時的最小執行緒數單位！！！非常重要，涉及到優化問題。 
  Threads Per Block指的是每個Block可排程的最大執行緒數，具體的後面會解釋。

• 然後切換到Memory選項卡 
   
  (ps:請再次無視840m垃圾的效能) 
  Toltal Global即視訊記憶體大小。 
  Bus Width指的是視訊記憶體位寬，單通道視訊記憶體32bit，64bit就是雙通道。 
  Clock Rate就是視訊記憶體頻率，900Mhz是低壓的GDDR3視訊記憶體，等效頻率就是通道數×視訊記憶體頻率，840m等效1800Mhz 
  Shared Per Block是共享記憶體，請牢記！！！！優化重點。在Maxwell架構中一個SM擁有16KB共享記憶體，840m有3個SM就是48KB。 
  Total Constant常量記憶體

• 最後切換到Performence選項卡 
   
  Host && Device資料交換速率在編寫CUDA程式時儘量少從主存中Copy矩陣到視訊記憶體，這個操作走的PCI通道，非常的費時間。即使是最新的PCI-E3.0頻寬也遠遠不及視訊記憶體頻寬，所以在編寫程式時儘早對資料交換進行。 
  Single-precision Float即單浮點運算效能，639GFlops是每秒可以執行639G次操作，在Maxwell架構中，每個核心具有兩個乘加器，即同時可以進行兩次乘法和加法操作。

• 下面具體介紹硬體架構

  1. 記憶體結構 
     （盜用中關村的美圖，懶得去角標，如若侵權請速與本人聯絡） 
      
     在編寫CUDA程式時常用的儲存結構：

     • __Host__修飾符為PC主存，即CPU可排程的記憶體。
     • __Global__修飾符為Device記憶體，也就是視訊記憶體，顯示卡可以排程的儲存空間。
     • __Shared__修飾符為Shared記憶體，共享記憶體，速度與L2 Cache相同級別，且延遲很低，讀取週期很短，但要注意Bank Conflict問題（後面會講到）。
     • __Device__修飾符通常用在核函式外，儲存在Device記憶體上，作為全域性變數。
     • __Constant__常量記憶體。 
       資料拷貝時，Shared先從 Host 走 PCI-E 通道拷貝到 Device 上，GPU讀取指令後先從L1、L2 Cache中尋找資料地址，若沒有再從 Device 上尋找，當然 Device 的讀取週期漫長。一般的做法是拷貝到 Shared Memory 上，GPU直接從 Shared Memory 取地址呼叫資料。

  2. 排程結構 
     在呼叫kernal函式時總體為一個Grid，Grid中含有Block，一個SM在執行時自動分配呼叫一些Block，每個Block中有大量的Thread。 
     GPU在運算時以一個Warp為單位，即32個Threads為單位，後面我們可以進行驗證排程過程。 
     Block可以是一維的二維的三維的，Thread也是如此，一般我們選用二維作為排程結構，圖中給出來索引的方式。 
     
     要注意的是：千萬不要讓Thread的x和y乘起來大於1024，因為一個block最多隻能排程1024個Threads。

• 程式效能分析

  想要寫出一個結合效能和準確率的程式時我們通常需要對這個程式進行效能分析，確保程式可以充分的利用硬體，而不是被記憶體的存取等原因阻塞計算。 
  一般CUDA程式分析看兩個： 
  1.GFlpos 
  2.存取頻寬 
  每個顯示卡都具有上述兩個瓶頸，可以說是不可逾越的鴻溝，當然顯示卡越好瓶口越大… 
  我們首先分析一下我們自己顯示卡的瓶頸在哪裡，當然無視最後我算出的結果，因為真的很低，這個瓶頸效能也可以從各種軟體當中獲取到。

  • GFlops瓶頸，單浮點數運算峰值效能，剛剛的CUDA-Z已經為我們測算出來了平均效能，我們自己手算一遍峰值效能： 
    840m在預設頻率下核心速率為1029MHZ，可以Boosting到1124MHZ（當然不可能一直超頻，所以只算預設頻率下的計算峰值效能），每顆核心可以同時執行2次乘加，也就是一顆CUDA Core每秒可以執行1029M × 2 = 2058M = 2.058G次乘加，一共有384個CUDA Cores，所以最後每秒可以執行的乘加數為2.058G × 384 = 790.272 GFLOPS。
  • 記憶體瓶頸，即頻寬瓶頸，為什麼會有這個瓶頸，因為GPU在運算時需要從Device Memory上取運算元，而取數操作非常慢，所以若程式執行緒可以保持一直並行取數的狀態，並將這些數用於計算，則效率會相當高，在GPU-Z中會顯示記憶體頻寬大小，這與視訊記憶體位寬和視訊記憶體頻率有關，我們也來手動計算一下： 
    840m的預設視訊記憶體顆粒頻率為900MHZ，等效工作頻率為1800MHZ，一次可以取64Bit大小的資料，我們在計算時一般都利用單浮點數進行計算，單浮點數有4個Byte，每個Byte有8個Bit，也就是一個Float型數為32Bit，也就是840m每次操作可以取2個Float型的數，共8Byte大小，故每秒可以取1800 × 8 = 14400MHZ = 14.4GHZ。Linux系統等效頻率可以在NVIDIA X SERVER裡面檢視。當然裡面你也能看到PCI-E及其緩慢的傳輸速率。 
    

• 在Ubuntu上配置CUDA

  建議在Linux上進行程式設計，因為你的工程可以直接MAKEFILE或者直接命令列進行NVCC編譯而不需要像Windows上那樣配置在VS201x上。 
  系統建議用16.04版本，因為這和之前的14.04一樣時LTS版本，有著較長的生命週期並且比較穩定，但為什麼我不推薦14.04版本，因為現在許多新版本的框架已經不支援老版本Ubuntu系統了。

  • 第一步要裝顯示卡驅動，Ubuntu相對方便很多。 
    開啟設定 
    
    找到最底下的軟體和更新 
    開啟附加驅動選項卡 
     
    選擇最新的NVIDIA binary driver，如果是筆記本帶有Inter整合顯示卡的使用者千萬千萬在選擇NVIDIA驅動的同時下面會有Intel的驅動，選擇成除預設外的另外一個，如果你忘了選擇Intel驅動就開始安裝NVIDIA驅動然後重啟的話…嗯…你可能需要重新安裝你的Ubuntu系統，別問我怎麼知道的，慘痛的親身經歷～ 
    安裝好後會重啟電腦，這時Dash中會出現我剛剛所說的NVIDIA X SERVER，當你不用獨立顯示卡時或者外出為了省電時可以在這個軟體裡面進行熱切換顯示卡，具體的操作在PRIME PROFILES選項卡中。 
    

  • 接著要配置CUDA 
    這裡我們要安裝CUDA8.0,為什麼不安裝9.0？ 
    因為現在許多深度學習框架對CUDA9.0的支援並不好，所以最穩定最保險就是安裝8.0版本，當然如果要進行深度學習開發還需要安裝Cudnn等附加包，這裡就不再贅述了。 
    首先要取NVIDIA官網下載CUDA8.0，選擇LINUX-x86_64-Ubuntu-16.04-runfile[local]，也就是本地的run檔案。 
    
    下載好以後用命令列執行run檔案。

    //執行run檔案
    sudo sh cuda_8.xx.xx_linux.run
    //編輯環境變數
    sudo gedit /etc/profile
    //在末尾加入以下環境變數
    export PATH=/usr/local/cuda-8.0/bin:$PATH
    export LD_LIBRARY_PATH=/usr/local/cuda-8.0/lib64$LD_LIBRARY_PATH

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7

  完成以上步驟後重啟電腦 
  檢測自己的CUDA是否配置成功只需ctrl+alt+T開啟命令列輸入

  nvcc

  • 1

  如果跳出來 
   
  或者並沒有報錯 
  則說明CUDA已經安裝完畢 
  至此所有提前需要了解的知識已經講完了，當然我預設讀者是掌握了C，C++並有一些程式設計功底的人，畢竟並行程式設計要考慮的東西遠遠比在CPU上寫程式來的多的多。

• 額外的程式環境

  在這裡我推薦直接用Sublime Text進行編寫，畢竟看著要比Vim舒服寫。而且在編寫較長程式時右上角的縮圖滾動是我認為比別的文字編輯器做的都好的地方。 
  
  這裡有官網下載傳送門：Download Sublime Text

• 小試牛刀

  終於要開始寫第一個CUDA程式了，清空腦袋，開始吧！ 
  建立一個CUDA檔案，CUDA檔案的字尾名是.cu，我們取名叫basic.cu 
  用sublime text開啟。 
  由於編譯器不能認識.cu字尾名，這樣語法就不能高亮顯示，我們在軟體的右下角選擇到C++模式。 
  一個CUDA程式包含兩個部分，一個是序列部分，一個時CUDA並行部分，序列部分在CPU上執行，CUDA部分則會拷貝到顯示卡上執行。 
  首先我們要引入標頭檔案，這裡我們使用Cpp，C，CUDA混編模式，所以需要引入三個標頭檔案：

#include <stdio.h>
#include <iostream>
#include <cuda_runtime.h>
using namespace std;

• 1
• 2
• 3
• 4

其中cuda_runtime.h是cuda函式所在的宣告檔案 
接下來我們要完成的任務是：將一個矩陣輸入到Global記憶體中，利用GPU全部加1後返回到Host記憶體進行輸出。 
第一步是需要在CPU中建立一個矩陣，我們一般使用一維動態陣列開闢，用二維的方式進行索引。 
先利用Malloc函式在CPU上開闢一塊空間，並全部賦值為1。

int size = 5;
float *a = (float*)malloc(sizeof(float)*size*size);
for (int i = 0; i < size; i++) {
    for (int j = 0; j < size; j++) {
        a[i*size+j] = 1.0f;
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7

然後需要在GPU上同樣開闢一個相同大小的空間以存放矩陣，這裡使用cudaMalloc函式。

float *a_cuda;
cudaMalloc((void**)&a_cuda,sizeof(float)*size*size);

• 1
• 2

接著，我們將矩陣從CPU上copy到GPU上。

cudaMemcpy(a_cuda,a,sizeof(float)*size*size,cudaMemcpyHostToDevice);

• 1

這時的a_cuda指向的是GPU上Device Memory上的一塊地址。 
那GPU要如何才能執行這一塊記憶體中的資料呢？ 
對，沒錯，就是使用核函式，也叫作Kernel函式。 
核函式的使用語法如下：

Function<<<griddim,blockdim,extern shared memory,GPU stream>>>(param...);

• 1

中間的引數可以控制核函式執行所佔用的資源。 
griddim表示呼叫的block塊數 
blockdim表示呼叫的thread數 
後面兩個引數分別表示動態定義共享記憶體大小和可使用的SM處理器數。 
那說到這裡，如何定義kernel呢？ 
kernal函式用__global__修飾符來修飾 
下面我們就來定義一個矩陣每個元素都加 1 的kernel函式 
在定義核函式之前先要考慮好需要呼叫多少block和thread，這裡時5×5的矩陣，我們可以使用1個block和25個thread排列為5×5thread陣列。 
核函式定義如下：

__global__ void addone(float *a) {  
    int tix = threadIdx.x;
    int tiy = threadIdx.y;
    int bdx = blockDim.x;
    int bdy = blockDim.y;
    a[tix*bdy+tiy] += 1;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7

其中threadIdx和blockDim是CUDA的內建變數，指示這當前正在排程的執行緒號和執行緒的數量。 
每個執行緒都會呼叫一次，故只需要將a矩陣對應位置的數值+1即可。 
接著在主函式中呼叫核函式。

dim3 grid(1, 1, 1), block(5, 5, 1);
addone<<<grid,block>>>(a_cuda);

• 1
• 2

dim3是一個CUDA內建的變數，是一個三維指示變數，分別對應這x，y，z三維，利用這個變數我們可以控制程式所呼叫的block和thread總量。 
由於沒有用到動態的shared memory也不去控制呼叫的SM核心數，所以後面兩個引數保持預設值。 
最後將執行完成的a_cuda拷貝到Host記憶體中。

cudaMemcpy(a,a_cuda,sizeof(float)*size*size,cudaMemcpyDeviceToHost);

• 1

為了矩陣輸出的方便，自己寫一個矩陣列印函式。

void MatrixPrint(float *mat, int rows, int cols) {
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            cout << setw(2) << mat[i*cols+j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9

這些也需要對應的格式化輸出標頭檔案

#include <iomanip>

• 1

到這裡就完成了整個CUDA程式的編寫。 
看一下整體程式碼：

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <cuda_runtime.h>
using namespace std;

void MatrixPrint(float *mat, int rows, int cols) {
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            cout << setw(2) << mat[i*cols+j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}

__global__ void addone(float *a) {
    int tix = threadIdx.x;
    int tiy = threadIdx.y;
    int bdx = blockDim.x;
    int bdy = blockDim.y;
    a[tix*bdy+tiy] += 1;
}
int main()
{
    int size = 5;
    float *a = (float*)malloc(sizeof(float)*size*size);
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            a[i*size+j] = 1.0f;
        }
    }
    MatrixPrint(a,size,size);
    float *a_cuda;
    cudaMalloc((void**)&a_cuda,sizeof(float)*size*size);
    cudaMemcpy(a_cuda,a,sizeof(float)*size*size,cudaMemcpyHostToDevice);

    dim3 grid(1, 1, 1), block(5, 5, 1);
    addone<<<grid,block>>>(a_cuda);
    cudaMemcpy(a,a_cuda,sizeof(float)*size*size,cudaMemcpyDeviceToHost);
    MatrixPrint(a,size,size);
    return 0;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43

那要如何執行這個程式呢？？？ 
將命令列cd到.cu檔案所在目錄，利用nvcc編譯器進行編譯，當然你要知道你的顯示卡計算能力時多少。 
整體編譯語句風格與GCC類似：

nvcc -gencode=arch=compute_50,code=\"sm_50,compute_50\" -o basic basic.cu

• 1

-gencode表示為計算能力xx而生成程式 
如果跳出的只有warning而沒有error的話說明程式通過編譯，可以執行。

./basic

• 1

直接執行即可 
 
發現在CUDA執行下已經將矩陣元素從1變為2。

• 進階之路

  基本的語法相信讀者已經初步掌握了，接下去就要對計算效能提出更高的要求。 
  之後的所有優化我們都圍繞著矩陣相乘來做，最後會對比CUBALS效能，CUBLAS是NVIDIA用其不開源的編譯器編譯出來的線性代數庫，彙編編寫的底層，執行速度飛一般的快……但是，沒有關係，經過我們的優化可以在一定範圍內超越CUBLAS的執行的速度。

• 矩陣相乘（第“3的0次方”講）

  大矩陣相乘是一個經典的可並行化的問題，也是一個在科學計算領域用的最多的一個演算法。 
  首先我們要捋清楚思路，矩陣是如何相乘的。 
   
  公式是這樣的，然而看起來還是不夠明瞭 
   
  這樣看起來就舒服多了 
  假如你還是沒有看懂或者忘了的話，線性代數知識可能需要好好回顧了。 
  那麼，我們需要先定義CPU程式碼，看看CPU是如何完成矩陣相乘的，當然這對於大多數程式設計師來說不過幾行程式碼的問題，非常簡單，在這裡我們直接將其封裝成函式：

void MatrixMul_host(float *a, int a_rows, int a_cols, float *b, int b_rows, int b_cols, float *c) {
    for (int i = 0; i < a_rows; i++) {
        for (int j = 0; j < b_cols; j++) {
            float t = 0;
            for (int k = 0; k < b_rows; k++) {
                t += a[i*a_cols+k]*b[k*b_cols+j];
            }
            c[i*b_cols+j] = t;
        }
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11

相信大家都可以看懂的～～～ 
當然想要執行這個函式還需要申請空間然後對矩陣賦值…and so on 
我們也封裝一些這樣的函式，簡便起見我們將所有的矩陣都賦值為-1和1兩個值。 
先是矩陣生成函式

void MatrixRandBin(float *mat, int rows, int cols) {
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if ((float)rand()/RAND_MAX > 0.5) {
                mat[i*cols+j] = 1.0f;
            }else {
                mat[i*cols+j] = -1.0f;
            }
        }
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11

然後是矩陣的比較函式，用於比較GPU最後跑出來的結果和CPU跑出來的是否一致，也就是驗證正確性，這相當重要，不能優化著優化著把結果搞錯了。

float MatrixCompare(float *a,float *b,int rows,int cols){
    float err=0;
    for (int i=0;i<rows;i++){
        for (int j=0;j<cols;j++){
            err+=abs(a[i*cols+j]-b[i*cols+j]);  
        }
    }
    return err;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9

到這裡為止，我們的輔助函式就定義完成了，接下去開始實現GPU版的矩陣相乘。 
在寫所有程式之前都應該想清楚程式如何執行，執行緒如何排程。 
+ 讓人一想就能想到的演算法是，既然有block，既然有thread，那麼每個block負責1行的計算，每個thread負責幾列的計算應該是比較好的。

thread1	thread2	thread3	……	thread1	thread2	……
block2
block3
block4
……
block100
block101
……

這時如何寫這個kernel函式呢？ 
Device Memory有一個特點，連續讀取的速度遠遠高於隨機讀取的速度，那什麼叫連續讀取，這裡就涉及到執行緒的排程問題。單個執行緒連續讀取一塊記憶體算是連續讀取嗎？錯了，在GPU執行時，是一個執行緒讀取完一組資料後直接呼叫下一個執行緒讀取，而非一個執行緒連續讀取，所以按照執行緒號排列的索引才是連續的索引。具體操作看kernel函式：

__global__ void MatrixMul_device(float *a, int a_rows, int a_cols, float *b, int b_rows, int b_cols, float *c) {
    int tix = threadIdx.x;
    int tiy = threadIdx.y;
    int bix = blockIdx.x;
    int biy = blockIdx.y;
    int bdx = blockDim.x;
    int bdy = blockDim.y;
    int gdx = gridDim.x;
    int gdy = gridDim.y;

    for (int i = tix; i < b_cols; i += bdx) {
        float sum = 0;
        for (int k = 0; k < a_cols; k++) {
            sum += a[bix*a_rows+k]*b[k*b_cols+i];
        }
        c[bix*a_cols+i] = sum;
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18

這個演算法看上去非常的簡潔，然後我們同樣的先要在CPU上開闢空間，用1000×1000的矩陣相乘來做實驗比較好，因為太大了CPU的驗證函式真的跑不動== 
然後也是同樣的拷貝到Device Memory上

int Matrixsize=1000;
float *a_host;
float *a_device;
float *b_host;
float *b_device;
float *result_host;
float *result_device;
float *result_cpu;
a_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
b_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
result_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
result_cpu = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
srand(0);
MatrixRandBin(a_host,Matrixsize,Matrixsize);
MatrixRandBin(b_host,Matrixsize,Matrixsize);
cudaMalloc((void**)&a_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMalloc((void**)&b_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMalloc((void**)&result_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMemcpy(a_device,a_host,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyHostToDevice);
cudaMemcpy(b_device,b_host,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyHostToDevice);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20

result_device用來儲存矩陣相乘後的GPU結果，所以也需要在GPU上開闢一段空間。 
接著呼叫kernel函式

cudaEvent_t start_device, stop_device;
float time_device;
cudaEventCreate(&start_device);
cudaEventCreate(&stop_device);
cudaEventRecord( start_device, 0 );
dim3 gridsize(1000,1,1);
dim3 blocksize(256,1,1);
MatrixMul_device<<<gridsize,blocksize>>>(a_device,Matrixsize,Matrixsize,b_device,Matrixsize,Matrixsize,result_device);
cudaEventRecord( stop_device, 0 );
cudaEventSynchronize( stop_device );
cudaEventElapsedTime( &time_device, start_device, stop_device );
cudaEventDestroy( start_device );
cudaEventDestroy( stop_device );
cout<<"gputime="<<time_device<<"ms"<<endl;

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14

上下兩個程式碼段包圍的都是用來測試GPU計算所用的時間的。 
因為有1000行，所以需要呼叫1000個block，thread的呼叫是有技巧的，讀者可以試試不同的值帶來的影響。第一點，設定thread最好在32的倍數，因為GPU是以warp作為排程單位，設定33這種，實際還是要呼叫2個warp，實則浪費了31個執行緒的計算能力。第二點，thread並不是開的越多越好，thread少，則程式並行度不夠，運算時沒有其他的warp進行I/O操作。thread多了，每個SM中暫存器數量有限，thread越多，所能夠並行的block就越少，最後還是會導致程式執行緩慢，達不到頻寬瓶頸。 
將執行的結果copy到Host中

cudaMemcpy(result_host, result_device,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyDeviceToHost);
cudaFree(a_device);
cudaFree(b_device);
cudaFree(result_device);
clock_t start_host = clock();
MatrixMul_host(a_host,Matrixsize,Matrixsize,b_host,Matrixsize,Matrixsize,result_cpu);
cout<<"cputime="<<(double)(clock() - start_host)/1000<<"ms"<<endl;
float err=MatrixCompare(result_cpu,result_host,Matrixsize,Matrixsize);
cout<<"err in gpu and cpu = "<<err<<endl;

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9

最後執行一下整個程式： 
 
可以看到，比CPU快了25倍！！！ 
讓我們來計算一下他的記憶體使用頻寬： 
一共1000行×1000列×中間1000次乘加×每次乘加需要取2個float型資料×每個float4個Byte=8GB資料 
8GB資料/用了0.2秒=40GB/S頻寬。。。。。 
？？？？？？？？？？？？？？？？？？？？？？？？？？？？ 
840m峰值頻寬只有14.4GB/S，怎麼會算出來40GB/S呢，難道哪裡有問題？？ 
其實這麼算並沒有錯，但是在計算時並沒有考慮到一個問題，就是GPU自動會將一些重複使用度高的資料快取到L1、L2 cache中，而快取到cache中所需要的頻寬幾乎可以忽略，所以在這個過程中所需要在Device Memory讀取的資料遠遠小於8GB。 
如何讓其GPU可以更好的利用CACHE呢？最簡單的做法就是不要頻繁的對快取的資料進行交換，並且儘可能一次讓其多讀取幾次，怎麼來實現？只需要減少thread數量就可以了，但是也是在保證並行執行的基礎上，由於呼叫了1000個block，似乎不必太關心並行度，這時我們將thread數減少到32再執行： 
 
OMGGGGGGGG!!!! 
比上一個版本快了將近3倍～～～～～～～ 
當然這只是在這樣簡單的程式上，GPU對程式的自主優化，但是我們有沒有想過自己來對這些進行優化呢？L1 Cache畢竟是有限的，還有和L1 Cache速度相仿的shared Memory 沒有用呢。當然，這些我們提到後面再說，接下來一講要說明Device Memory的對齊讀取。

矩陣相乘（第“In(e的2次方)”講）

在上一個版本的程式中，對GPU的Cache優化和Device Memory的連續存取有了初步的瞭解，但是GPU對Device Memory還有一個特點，GPU對它是按照記憶體對齊讀取的，在Compute Capability 5.0中，記憶體按照128個float進行對齊，比如說讀取第（0，513）元素時候可以從0的地址開始讀取512個地址，只需要少量定址就可以了，但是在讀取（1,513）時候還能這麼讀取嗎？ 
1000×1000的矩陣，第二行是1001個元素，已經不是和128對齊了，故GPU不能用PITCH方式進行讀取，而還是從0地址開始讀取，讀取1000×1+512個元素才能定址到，若是讀取最後一個則需要將整個矩陣元素全部定址才能讀取，所以記憶體對齊非常重要。若是每一行都與128的倍數對齊的話，CPU就可以按照PITCH讀取，也就是讀（1,513）時只需要從第二行開頭開始讀取，讀取512個元素即可，若是讀取最後一個元素，只需要用1000次索引到最後一行的開頭，然後再索引全部最後一行即可，時間複雜讀從o(m*n)變成了o(m+n)。 
怎麼對齊呢？ 
可以手動對齊，也就是（（1000-1）/128+1）× 128 = 1024。也就是需要開闢1000×1024的矩陣才能對齊記憶體，難道每次開闢時候還需要這麼麻煩的計算？ 
當然NVIDIA已經為我們考慮好了，記憶體對齊後的一行成為一個Pitch，利用cudaMallocPitch函式，CUDA會為我們自動分配Pitch不需要自己手動計算。在copy視訊記憶體時需要用cudaMemcpy2D函式，這個函式會自動匹配Pitch來進行copy。此時記憶體拷貝變成如下：

size_t pitch_a_device, pitch_b_device, pitch_result_device;
cudaMallocPitch((void**)&a_device , &pitch_a_device , sizeof(float) *Matrixsize , Matrixsize);
cudaMallocPitch((void**)&b_device , &pitch_b_device , sizeof(float) *Matrixsize , Matrixsize);
cudaMallocPitch((void**)&result_device , &pitch_result_device , sizeof(float) *Matrixsize , Matrixsize);
cudaMemcpy2D(a_device,pitch_a_device,a_host,sizeof(float) *Matrixsize ,sizeof(float) *Matrixsize, Matrixsize,cudaMemcpyHostToDevice);
cudaMemcpy2D(b_device,pitch_b_device,b_host,sizeof(float) *Matrixsize ,sizeof(float) *Matrixsize, Matrixsize,cudaMemcpyHostToDevice);

• 1
• 2
• 3
• 4
• 5
• 6

這些size_t就是以Byte為單位的Pitch長度。如果要轉換為float型多大空間則需要/sizeof(float)來得到。 
拷貝到Host的函式是

cudaMemcpy2D(result_host,sizeof(float) *Matrixsize, result_device,pitch_result_device,sizeof(float) *Matrixsize , Matrixsize,cudaMemcpyDeviceToHost);

• 1

kernal函式的索引需要做出改變

__global__ void MatrixMul_device(float *a,int a_pitch,int a_rows,int a_cols,float *b,int b_pitch,int b_rows,int b_cols,float *c,int c_pitch){
    int tix=threadIdx.x;
    int tiy=threadIdx.y;
    int bix=blockIdx.x;
    int biy=blockIdx.y;
    int bdx=blockDim.x;
    int bdy=blockDim.y;
    int gdx=gridDim.x;
    int gdy=gridDim.y;
    int i,k;
    for(i=tix;i<b_cols;i+=bdx){
        float sum=0;
        for(k=0;k<a_cols;k++){
            sum+=a[bix*a_pitch+k]*b[k*b_pitch+i];
        }
        c[bix*a_pitch+i]=sum;
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18

用picth來作為索引的列寬。 
呼叫核函式也需要做出改變。

MatrixMul_device<<<gridsize,blocksize>>>(a_device,pitch_a_device/sizeof(float),Matrixsize,Matrixsize,b_device,pitch_b_device/sizeof(float),Matrixsize,Matrixsize,result_device,pitch_result_device/sizeof(float));

• 1

同樣的用nvcc編譯一下.cu檔案。 
看一下結果。 
 
又減少了1/5的時間，由於Cache快取了整個操作，實際在Device Memory上進行操作非常少，所以看起來提升不明顯，其實整個操作可以帶來大約3倍頻寬的提升。

• #### 矩陣相乘（第“(int)(100/33)”講） 
  這一講我們將呼叫類似於L1 Cache速度的shared memory來手動快取一些操作，來更多的挖掘整個顯示卡的計算能力。 
  在用shared memory時我們的目標是將重複計算快取到shared memory中，因為本身快取操作也是一種讀取操作，如果沒有大量的重複呼叫這些地址則帶來的提升收效甚微。 
  在這種方式操作的矩陣相乘中，A每一個行的都需要與B的每一列進行對應相乘，每次都需要讀取A所操作的一行，在這之中就有大量的重複讀取，若能將A所乘的一行讀取到shared memory中，就可以減少大量的頻寬佔用。 
  再次改變核函式內部：

__shared__ float a_shared[1000];
for(int j=tix;j<a_cols;j+=bdx){
    a_shared[j]=a[bix*a_pitch+j];
}
__syncthreads();
int i,k;
for(i=tix;i<b_cols;i+=bdx){
    float sum=0;
    for(k=0;k<a_cols;k++){
        sum+=a_shared[k]*b[k*b_pitch+i];
    }
    c[bix*a_pitch+i]=sum;
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13

首先申請1000空間的shared memory，用__shared__修飾。 
然後利用執行緒將所需要操作的空間copy進去，此時需要同步執行緒時間才能進入下一步操作，因為執行緒在這一步操作與下一步操作無關，下一步的操作讀取必須實在shared memory讀取完成後的，所以需要用__syncthreads()同步block內執行緒，當這個block內所有執行緒都完成copy操作後再進行乘法操作即可。 
這一版本的程式需要多少時間花費呢？ 
 
又減少了不少的時間，這時L1快取的作用減少了許多，因為我們已經把大部分重複操作都讀進共享記憶體中，可以看出我們自己寫的記憶體優化操作其實要比GPU自動進行的好，當然也是由於L1、L2的原因，並不能直接計算頻寬，實際上，這個操作可以帶來大約2-5倍的頻寬節省，節省非常的大。此時已經比cpu快了100倍。當然還可以更快速。 
這三講把矩陣相乘優化了，但是在優化記憶體時也僅僅把A部分的矩陣快取至shared memory，B矩陣部分並沒有過多的優化，可以藉助分塊矩陣演算法進行優化，這個會在以後講到，下一講將針對一個比較有用的矩陣進行特殊的優化。

• 矩陣相乘（第“2×(最小自然數+2)”講）

  除了利用分塊矩陣演算法將B矩陣也快取進shared memory似乎就沒有其他可以對記憶體頻寬帶來收益的操作了。所以程式的解決必須開拓思維，分析整個程式的核心部分。矩陣乘法的核心部分就是A的行和B的列對應相乘部分，這部分如何可以減少時間呢？ 
  在深度學習中，幾乎都用的單浮點型矩陣相乘，而如果矩陣只有1與-1值呢？我們稱為二進位制矩陣，用其進行的操作必然會帶來精度的降低，但是如果用單浮點矩陣在需要實時運算的且裝置能力較低的嵌入式中使用大型矩陣相乘操作，無疑會帶來

  1. 視訊記憶體空間
  2. 計算耗時

  這些問題，用二進位制矩陣來代替現有的矩陣框架，怎麼相乘？1和-1相乘為-1不就相當於1與0進行XNOR操作嗎？這時矩陣只需要佔用1位進行儲存，視訊記憶體空間減少32倍，也就是利用int型可以存下32個數字，並且只需要1步同或操作就可以完成32個數字的乘法操作，相當迅捷，這就為我們接下來的優化提供了思路。 
  第一步還是先捋清整個演算法的思路： 
  1.申請記憶體空間 
  2.copy到視訊記憶體 
  3.A按行，B按列，每32個數字利用二進位制表示儲存成int型 
  4.A與B相乘，把乘法操作轉化為同或操作 
  5.統計異或結果的和時利用popcount演算法 
  6.得到結果，返回Host 
  7.比較結果 
  開始一步步實現： 
  我們先利用之前不用shared memory版本的程式框架來實現。 
  先申請記憶體空間並copy到視訊記憶體：

int Matrixsize=1000;
float *a_host;
float *b_host;
float *result_host;
a_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
b_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
result_host = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
srand(0);
MatrixRandBin(a_host,Matrixsize,Matrixsize);
MatrixRandBin(b_host,Matrixsize,Matrixsize);
// cout<<MatrixCopysize<<endl;

float *a_device;
float *b_device;
float *result_device;
cudaMalloc((void**)&a_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMalloc((void**)&b_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMalloc((void**)&result_device,sizeof(float) *Matrixsize * Matrixsize);
cudaMemcpy(a_device,a_host,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyHostToDevice);
cudaMemcpy(b_device,b_host,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyHostToDevice);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20

建立一個函式封裝整個操作：

void MatrixMul_device(float *a,float *b,int a_rows,int a_cols,int b_cols,float *result){}

• 1

第一個要實現的就是將普通矩陣A和B按32位轉換成int型。 
先是開始的準備階段，進行pitch對齊，由於需要32位轉換一次，所以A的每一行（B中每一列）也必須是32的倍數，空餘的進行0填充，如何進行0填充？ 
比較好的思路就是先申請一塊32倍數對齊的空間，然後全部置零，再將原矩陣copy進去。

int BINSIZE=32;//size of bin2int, 32 means 0000 0000 0000 0000 0000 0000 0000 0000
int MaxBS=(a_cols-1)/BINSIZE+1;
int a_cols_Copysize=MaxBS*BINSIZE;
float *a_device;//a_rows * a_cols_Copysize
float *b_device;//a_cols_Copysize * b_cols
size_t pitch_a_device, pitch_b_device;
cudaMallocPitch((void**)&a_device , &pitch_a_device , sizeof(float) *a_cols_Copysize , a_rows);
cudaMallocPitch((void**)&b_device , &pitch_b_device , sizeof(float) *b_cols , a_cols_Copysize);

cudaMemset(a_device, 0, pitch_a_device * a_rows);
cudaMemset(b_device, 0, pitch_b_device * a_cols_Copysize);
cudaMemcpy2D(a_device,pitch_a_device,a,sizeof(float) *a_cols ,sizeof(float) *a_cols, a_rows,cudaMemcpyDeviceToDevice);
cudaMemcpy2D(b_device,pitch_b_device,b,sizeof(float) *b_cols ,sizeof(float) *b_cols, a_cols,cudaMemcpyDeviceToDevice);

int *a_device_bin;
int *b_device_bin;
size_t pitch_a_device_bin, pitch_b_device_bin;
cudaMallocPitch((void**)&a_device_bin , &pitch_a_device_bin , sizeof(int) *MaxBS , a_rows);
cudaMallocPitch((void**)&b_device_bin , &pitch_b_device_bin , sizeof(int) *b_cols , MaxBS);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19

然後就是轉換操作，由於A，B矩陣的轉換不同，故需要寫兩個函式。

//horizontal
__global__ void AMatrix2Bin(float *a,int *a_bin,int a_rows,int pitch_a,int pitch_a_bin,int MaxBS,int BINSIZE){
    int tix=threadIdx.x;
    // int tiy=threadIdx.y;
    int bix=blockIdx.x;
    // int biy=blockIdx.y;
    int bdx=blockDim.x;
    // int bdy=blockDim.y;
    int gdx=gridDim.x;
    // int gdy=gridDim.y;


    int maxThreads=MaxBS*a_rows;
    for(int id = bix*bdx+tix; id < maxThreads; id+=gdx*bdx) {
        int rid=id/MaxBS;
        int cid=id%MaxBS;

        int Integer=0;
        int base=1;
        for (int i=0;i<BINSIZE;i++){
            if (a[rid*pitch_a+(cid+1)*BINSIZE-1-i]==1.f){
                Integer+=base;
            }
            base=base<<1;
        }

        a_bin[rid*pitch_a_bin+cid]=Integer;
    }

}
//vetical
__global__ void BMatrix2Bin(float *b,int *b_bin,int b_cols,int pitch_b,int pitch_b_bin,int MaxBS,int BINSIZE){
    int tix=threadIdx.x;
    // int tiy=threadIdx.y;
    int bix=blockIdx.x;
    // int biy=blockIdx.y;
    int bdx=blockDim.x;
    // int bdy=blockDim.y;
    int gdx=gridDim.x;
    // int gdy=gridDim.y;

    int maxThreads=MaxBS*b_cols;
    for(int id = bix*bdx+tix; id < maxThreads; id+=gdx*bdx) {
        int cid=id/MaxBS;
        int rid=id%MaxBS;

        int Integer=0;
        int base=1;
        for (int i=0;i<BINSIZE;i++){
            if (b[((rid+1)*BINSIZE-1-i)*pitch_b+cid]==1.f){
                Integer+=base;
            }
            base=base<<1;
        }

        b_bin[rid*pitch_b_bin+cid]=Integer;
    }

}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56
• 57
• 58
• 59

其中有一個小技巧，就是利用>>左移符號，左移一位相當乘2,當然這樣操作會比直接利用浮點型計算要快。 
這個函式並不需要多少時間執行，時間複雜度並不高，且記憶體讀取次數也非常少，故不將其計入總執行時間。 
呼叫兩個函式。

AMatrix2Bin<<<GS_BIN,BS_BIN>>>(a_device , a_device_bin , a_rows ,
    pitch_a_device/sizeof(float) , pitch_a_device_bin/sizeof(int) , MaxBS , BINSIZE);
BMatrix2Bin<<<GS_BIN,BS_BIN>>>(b_device , b_device_bin , b_cols ,
    pitch_b_device/sizeof(float) , pitch_b_device_bin/sizeof(int) , MaxBS , BINSIZE);

• 1
• 2
• 3
• 4

接著就是寫主kernel函式，當然需要先準備popcount函式，popcount函式版本非常多，本人採用利用儲存空間來換取效能的popcount函式，先將256個計數位放到Device memory中，這裡需要用到cudaMemcpyToSymbol（）函式，這個函式可以用來copy全域性記憶體和常量記憶體等。然後在Device memory上建立一個kernel可呼叫的popcount函式。 
申明全域性變數和函式：

__device__ unsigned char __popcount_tab_device[256];//__constant__ is slower than __device__
__device__ int popcount (int x) {
  return __popcount_tab_device[(x >>  0) & 0xff]  
  + __popcount_tab_device[(x >>  8) & 0xff]  
  + __popcount_tab_device[(x >> 16) & 0xff] 
  + __popcount_tab_device[(x >> 24) & 0xff];
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7

拷貝計數位：

const unsigned char __popcount_tab[] = {
  0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8,
};
cudaMemcpyToSymbol(__popcount_tab_device, __popcount_tab, sizeof(__popcount_tab));

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11

呼叫主kernel的準備工作，建立結果容器：

float *result_device;//a_rows * b_cols
size_t pitch_result_device;
cudaMallocPitch((void**)&result_device , &pitch_result_device , sizeof(float) *b_cols , a_rows);

• 1
• 2
• 3

最主要的工作就是建立核函式：

__global__ void MatrixMulXnor(int *a,int *b,int a_rows,int a_cols,
    int b_cols,float *result,int pitch_a,int pitch_b,
    int pitch_result,int BINSIZE,int RealMidSize){

    int tix=threadIdx.x;
    // int tiy=threadIdx.y;
    int bix=blockIdx.x;
    // int biy=blockIdx.y;
    int bdx=blockDim.x;
    // int bdy=blockDim.y;
    int gdx=gridDim.x;
    // int gdy=gridDim.y;

    int rest=(BINSIZE*a_cols-RealMidSize);

    for(int j=tix;j<b_cols;j+=bdx){
        // printf("i=%d ; j=%d\n",i,j);
        int sum=0;
        for(int k=0;k<a_cols;k++){
            int bin=(a[bix*pitch_a+k]^b[k*pitch_b+j]);
            int negnum=popcount(bin);
            int posnum=BINSIZE-negnum;
            //calculate ignores the rest of BINSIZE if the Matsize cant devided by BINSIZE ,it can cause err
            //(10/00)'(01/00) should be 0000 but it is 0011,so 1+1 is trash in the result.and it can cause a_rows*b_cols times. 
            sum+=(posnum-negnum);
        }
        result[bix*pitch_result+j]=sum-rest;
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29

為什麼需要最後減去rest，rest就是那些填充0的地方，由於只有XOR異或運算子，所以運算出的結果是相反的，但本身是0的地方也會相反成1,這時popcount函式就會多計算末尾填充的那些位置，需要減去他們才能得到正確結果。 
呼叫一下：

cudaEvent_t start_device, stop_device;
float time_device;
cudaEventCreate(&start_device);
cudaEventCreate(&stop_device);
cudaEventRecord( start_device, 0 );

dim3 BS_MM(32,1,1);
dim3 GS_MM(1000,1,1);
MatrixMulXnor<<<GS_MM,BS_MM>>>(a_device_bin , b_device_bin , a_rows , MaxBS , b_cols ,
 result_device , pitch_a_device_bin/sizeof(int) , pitch_b_device_bin/sizeof(int) , 
 pitch_result_device/sizeof(float) , BINSIZE , a_cols);

cudaEventRecord( stop_device, 0 );
cudaEventSynchronize( stop_device );
cudaEventElapsedTime( &time_device, start_device, stop_device );
cudaEventDestroy( start_device );
cudaEventDestroy( stop_device );
cout<<"gputime="<<time_device<<"ms"<<endl;

cudaMemcpy2D(result,sizeof(float) *b_cols, result_device,pitch_result_device,sizeof(float) *b_cols , a_rows ,cudaMemcpyDeviceToDevice);

cudaFree(a_device);
cudaFree(b_device);
cudaFree(a_device_bin);
cudaFree(b_device_bin);
cudaFree(result_device);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26

到這裡函式就算封裝完成了，只需要外部呼叫一下即可：

//run in gpu warp in C code
MatrixMul_device(a_device,b_device,Matrixsize,Matrixsize,Matrixsize,result_device);

cudaMemcpy(result_host, result_device,sizeof(float) *Matrixsize * Matrixsize,cudaMemcpyDeviceToHost);
cudaFree(a_device);
cudaFree(b_device);
cudaFree(result_device);
// MatrixPrint(result_host,Matrixsize,Matrixsize);

//run in cpu
float *result_cpu;
result_cpu = (float*) malloc(sizeof(float) * Matrixsize * Matrixsize);
clock_t start_host = clock();
MatrixMul_host(a_host,Matrixsize,Matrixsize,b_host,Matrixsize,Matrixsize,result_cpu);
cout<<"cputime="<<(double)(clock() - start_host)/1000<<"ms"<<endl;
// MatrixPrint(result_cpu,Matrixsize,Matrixsize);


//compare value of gpu and cpu
float err=MatrixCompare(result_cpu,result_host,Matrixsize,Matrixsize);
cout<<"err in gpu and cpu = "<<err<<endl;

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22

執行一下試試，此時你會非常吃驚的發現，沒有寫快取且沒有做任何優化的竟然比上面做了大量優化工作的演算法還要快近一倍。 
 
由於中間對應相乘的計算減少了，GPU自動優化的空間大量減少，此時我們可以試著計算一下這個程式的頻寬使用，當然真實的肯定還是要比這個低的多： 
1000行×1000列×1000/32次核心計算×每次2個int型記憶體讀取×int為4位元組/0.03S運算=8GB/S的頻寬。。。真實的可能連一半都不到，畢竟在14.4GB峰值上連40GB都能搞出來的優化，真實的肯定大打折扣。一定要記住之前的40GB/S只是算出來的，真實的其實在GPU cache以後不需要讀取那麼多次記憶體，所以40只是一個虛的數字。

矩陣相乘（第“5×(1000的0次方)”講）

當然這一講肯定是試著將普通矩陣相乘最好的優化方式套到二進位制矩陣相乘，也就是加上shared memory來快取A的一行而已，但是現在一行才30個int型作用，可想而知，速度提升肯定不會很大。在這裡直接上改動的核函式和結果，具體不做講解。

__global__ void MatrixMulXnor(int *a,int *b,int a_rows,int a_cols,
    int b_cols,float *result,int pitch_a,int pitch_b,
    int pitch_result,int BINSIZE,int RealMidSize){
    int tix=threadIdx.x;
    // int tiy=threadIdx.y;
    int bix=blockIdx.x;
    // int biy=blockIdx.y;
    int bdx=blockDim.x;
    // int bdy=blockDim.y;
    int gdx=gridDim.x;
    // int gdy=gridDim.y;
    extern __shared__ int a_shared[];
    for(int j=tix;j<a_cols;j+=bdx){
        // printf("i=%d ; j=%d\n",i,j);
        a_shared[j]=a[bix*pitch_a+j];
    }
    __syncthreads();
    int rest=(BINSIZE*a_cols-RealMidSize);

    for(int j=tix;j<b_cols;j+=bdx){
        // printf("i=%d ; j=%d\n",i,j);
        int sum=0;
        for(int k=0;k<a_cols;k++){
            int bin=(a_shared[k]^b[k*pitch_b+j]);
            int negnum=popcount(bin);
            int posnum=BINSIZE-negnum;
            //calculate ignores the rest of BINSIZE if the Matsize can't devided by BINSIZE ,it can cause err
            //(10/00)'(01/00) should be 0000 but it is 0011,so 1+1 is trash in the result.and it mislead a_rows*b_cols times. 
            sum+=(posnum-negnum);
        }
        result[bix*pitch_result+j]=sum-rest;
    }
}

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33

結果如下，並沒有多少改進，所以必須尋求它法。 

矩陣相乘 倒數第二講

之前曾經提到過分塊矩陣乘法，即將矩陣分成可以相乘的幾個塊，可以按照乘法規則運算： 
 
 
這樣矩陣乘法核心部分就可以轉化成將一部分A矩陣和一部分B矩陣做乘法然後相加，從而A、B都可以存於shared memory，且不像上一版本只能將A的一行大約32個存於shared memory導致頻寬提升不大。 
分析演算法的改動思路： 
1.思考需要多大的矩陣作為分塊矩陣： 
由於執行緒最多開1024個故最多用32×32的分塊矩陣，一次操作需要快取32×32×2(A與B)×4(int型大小) = 8KB，一個SM只能呼叫16KBshared memory，直接導致程式並行度下降，故選擇開16×16作為分塊矩陣大小。 
2.思考迴圈結構對於GFlops的影響： 
迴圈結構也是需要耗費計算能力的，一般程式的核心語句計算能力如果佔整個kernel程式越多，則最後越接近峰值計算能力，而迴圈則導致每次計算完乘法就需要計算迴圈的條件，導致大量計算能力被分出去，如果是固定的結構且迴圈數量並不大，直接展開迴圈進行編寫程式效率會大幅度提高，而16×16的矩陣相乘正式如此。 
3.幾乎每次都會讀取popcount的計數位，而其存放在Device memory中，將其存放在shared memory也會帶來巨量的頻寬提升。 
開始寫封裝程式： 
為了避免分塊矩陣在邊緣計算溢位（用if語句判斷將導致分支結構，而GPU並沒有分支預測能力，所以分支結構要儘量減少），直接將矩陣擴充為16的倍數，方便操作。

int RectBlockSize = 16;
dim3 RectBlockNum_a_bin((a_rows-1)/RectBlockSize+1, (MaxBlocks-1)/RectBlockSize+1, 1);//with block multiply
dim3 RectBlockNum_b_bin((MaxBlocks-1)/RectBlockSize+1, (b_cols-1)/RectBlockSize+1, 1);
int *a_bin;
int *b_bin;
size_t Pitch_a_bin, Pitch_b_bin;
cudaMallocPitch((void**)&a_bin , &Pitch_a_bin , sizeof(int)*RectBlockSize*RectBlockNum_a_bin.y, RectBlockSize*RectBlockNum_a_bin.x);
cudaMallocPitch((void**)&b_bin , &Pitch_b_bin , sizeof(int)*RectBlockSize*RectBlockNum_b_bin.y, RectBlockSize*RectBlockNum_b_bin.x);
cudaMemset(a_bin, 0, Pitch_a_bin*RectBlockSize*RectBlockNum_a_bin.x);
cudaMemset(b_bin, 0, Pitch_b_bin*RectBlockSize*RectBlockNum_b_bin.x);
dim3 BS_BIN(512,1,1);
dim3 GS_BIN(6,1,1);
AMatrix2Bin<<< GS_BIN, BS_BIN >>>(a_copy, a_bin, 
    Pitch_a_copy/sizeof(float), Pitch_a_bin/sizeof(int), a_rows, MaxBlocks, BINSIZE);
BMatrix2Bin<<< GS_BIN, BS_BIN >>>(b_copy, b_bin, 
    Pitch_b_copy/sizeof(float), Pitch_b_bin/sizeof(int), b_cols, MaxBlocks, BINSIZE);
cudaFree(a_copy);
cudaFree(b_copy);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18

申請矩陣計算結果空間並將popcount計數位copy到Device memory，因為想copy到shared memory就必須先將資料存於視訊記憶體。

float *result_bin;//a_rows * b_cols
size_t Pitch_result_bin;
cudaMallocPitch((void**)&result_bin , &Pitch_result_bin , sizeof(float)*RectBlockSize*RectBlockNum_b_bin.y, RectBlockSize*RectBlockNum_a_bin.x);
const unsigned char __popcount_tab[] = {
  0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8,
};
unsigned char *__popcount_tab_copy;
cudaMalloc((void**)&__popcount_tab_copy, sizeof(__popcount_tab));
cudaMemcpy(__popcount_tab_copy, __popcount_tab, sizeof(__popcount_tab), cudaMemcpyHostToDevice);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16

最重要的建立kernel函式：

_global__ void MatrixMulXnor(int *a, int *b, float *result, unsigned char *__popcount_tab,
    int pitch_a, int pitch_b, int pitch_result,
    int midBlocks, int BINSIZE, int RealMidSize)

• 1
• 2
• 3

最後呼叫並將結果copy出來：

cudaEvent_t start_device, stop_device;
float time_device;
cudaEventCreate(&start_device);
cudaEventCreate(&stop_device);
cudaEventRecord(start_device, 0);
dim3 BS_MM(RectBlockSize, RectBlockSize, 1);
dim3 GS_MM(RectBlockNum_a_bin.x, RectBlockNum_b_bin.y, 1);
MatrixMulXnor<<< GS_MM, BS_MM >>>(a_bin, b_bin, result_bin, __popcount_tab_copy,
    Pitch_a_bin/sizeof(int), Pitch_b_bin/sizeof(int), Pitch_result_bin/sizeof(float),
    RectBlockNum_a_bin.y, BINSIZE, a_cols);
cudaEventRecord( stop_device, 0 );
cudaEventSynchronize( stop_device );
cudaEventElapsedTime( &time_device, start_device, stop_device );
cudaEventDestroy( start_device );
cudaEventDestroy( stop_device );
cout<<"gputime="<<time_device<<"ms"<<endl;
cudaMemcpy2D(result,sizeof(float) *b_cols, result_bin,Pitch_result_bin,sizeof(float) *b_cols , a_rows ,cudaMemcpyDeviceToDevice);
cudaFree(a_bin);
cudaFree(b_bin);
cudaFree(result_bin);

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20

執行一下，真的會嚇一跳！！！ 
 
比上一個版本快了將近8倍，比CPU快了1000倍！！！這只是在840m的顯示卡上進行的運算啊。 
等會等會，難道你以為這樣就結束了？？？？？？？ 
錯了，這個程式有一個問題還是沒有考慮到！！！！！！ 
就是shared memory的bank conflict問題。

矩陣相乘 最後一講

什麼是bank conflict？ 
在讀取shared memory時，cpu是按照bank進行讀取，bank一般分為32個。當一個warp讀取到相同的bank的不同元素時就會造成bank conflict。 
現在有一個shared mem A[2][32]，當索引為[tix]時，由於一個warp排程按照32個執行緒，故會利用兩個warp排程完第一行以後再排程第二行，每一列為一個bank的話，一個warp剛好佔據每個bank，不會bank conflict。但當索引為[2*tix]時，warp中的第16個執行緒與第1個執行緒都指向A[x][0]，讀取了同一個bank中的不同元素，造成兩個執行緒必須序列執行，其他執行緒也一樣，這就減少了一半的並行度。當然每次都讀取同一個bank的相同元素並不會導致bank conflict，因為GPU具有廣播機制，會將地址廣播給其他執行緒。 
在上一版本程式中，在索引時是tix × 列寬 + tiy，出現tix的倍數則說明導致了bank conflict。解決辦法有： 
1.將share memory多申明一位，這就讓本身衝突的bank讀取進行了錯位處理。 
2.x與y索引倒置，也就是讓warp對tiy進行而非對tix進行。 
顯然第二種方法最簡便，只需在定義tix和tiy的時候進行調換，即tix為原本的tiy，tiy為原本的tix。 
這也就是為什麼在使用cublas時候需要對矩陣進行轉置再寫入函式，因為cublas是按列主序進行存取，這就避免了bank conflict現象，這也就是為什麼CUDA更加支援FORTRAN語言的原因，FORTRAN語言就是按列主序進行矩陣儲存的。 
看一下結果： 
 
提升三分之一的速度，非常的有效，因為只需要改兩個變數而已。 
最後計算一下頻寬： 
由於全部利用shared memory快取，故頻寬接近與真實頻寬。 
（16×16×4×2）×ceil（1000/16）×ceil（1000/16）×ceil（1000/32/16）/0.003S=5.59GB/S 
看到這裡可能有人就會說 太低了 離14.4GB/S還遠，但其實這個頻寬已經不取決於讀取效能了，我們可以用8×8的分塊矩陣試試。 
附上8×8的運算結果圖： 
 
與上面16×16的耗時幾乎相同，但是我們計算一下頻寬： 
（8×8×4×2）×ceil（1000/8）×ceil（1000/8）×ceil（1000/32/8）/0.003S=11.1GB/S 
這幾乎已經接近峰值頻寬了，因為在其他地方還有許多記憶體交換，比如說popconnt計數位等，當然也可以再小用4×4： 
 
這時耗時已經上升了，說明頻寬已經達到峰值了，計算一下： 
（4×4×4×2）×ceil（1000/4）×ceil（1000/4）×ceil（1000/32/4）/0.0098S=13.6GB/S 
可以看到，加上其他的運算操作，已經達到14.4GB/S的峰值頻寬了。

到這裡，矩陣優化就講完了，最後是對比CUBLAS的效能。在中型、小型矩陣方面（5000×5000以下），大幅度領先CUBLAS，但不得不佩服CUBLAS在超大型矩陣的優化能力，幾乎時間複雜度無變化。過些時間會有比較圖與更加精確的資料對比。

感謝大家的閱讀，若有錯誤，歡迎批評職責，若有問題，歡迎交流。