CUDA教程一、認識cuda

什么是cuda

統一計算設備架構（Compute Unified Device Architecture, CUDA），是由NVIDIA推出的通用並行計算架構。解決的是用更加廉價的設備資源，實現更高效的並行計算。

點擊下面鏈接就可以下載cuda。我個人使用的是10.2版，截止到目前官方已經發布了11.0版。

有人就問了，std::thread它不香嗎，為什么要用cuda？

別忘了，cuda是英偉達（NVIDIA）推出的——這個英偉達是何方神聖？沒聽說過英偉達，也應該聽說過GPU了吧。沒錯，提出GPU概念的，正是英偉達。和中央處理器（Central Processing Unit, CPU）相對，圖形處理器（Graphics Processing Unit, GPU）是顯卡的核心芯片。而cuda正是暴露了英偉達開發的GPU的編程接口。

幾乎所有的編程語言，不使用特定框架，都只能實現CPU編程——std::thread也是將線程開在CPU中。而不同於每一位程序員都接觸過的CPU編程，GPU編程可以使用更多的流處理器、更多的線程數。

舉個例子，打開我們的任務管理器，我們可以查看自己電腦CPU的線程數：

 而GPU的核數和線程數更是豐富，只要合理使用，幾乎沒有用得完的情況：

如何合理使用，才最高效、利用率最高呢？

從矩陣相乘說起

舉一個簡單的例子，計算兩個矩陣AB相乘的結果：

 我們知道：

如果使用單線程編程，我們可以這么寫：

template <typename T>
void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) {  3     assert(A.Width() == B.Height());  4  C.Set_size(A.Height(), B.Width());  5     for(size_t i = 0; i < A.Height(); ++i) {  6         for(size_t j = 0; j < B.Width(); ++j) {  7             T ans = 0;  8             for(size_t k = 0; k < A.Width(); ++k) {  9                 ans += A[i][k] * B[k][j]; 10  } 11             C[i][j] = ans; 12  } 13  } 14 }

這是毋庸置疑的，相信即使是編程初學者也可以寫出。當然它在局部性上有待優化，不過優化空間非常非常小。

如果是多線程編程呢？

由於多個線程同時寫入同一個變量是很危險的，但同時讀取同一個變量不會出現問題，所以如果不想引入非常耗時的互斥鎖，我們就必須根據可寫入變量的數量來划分線程數——比如上式兩個矩陣的乘積是的方陣，因此我們可以分成四個線程來計算，每個線程只計算結果矩陣的一個元素。代碼如下：

template <typename T>
void thrd_MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */,  3  size_t row, size_t col) {  4 //每一個線程因row和col的不同而被區分
    T ans = 0;  6     for(size_t k = 0; k < A.Width(); ++k) {  7         ans += A[row][k] * B[k][col];  8  }  9     C[row][col] = ans; 10 } 11 
template <typename T>
void MatrixProduct(threadPool& tp, matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) { 14     assert(A.Width() == B.Height()); 15  C.Set_size(A.Height(), B.Width()); 16 //將每個計算單元的任務放進線程池中
    for(size_t i = 0; i < A.Height(); ++i) { 18         for(size_t j = 0; j < B.Width(); ++j) { 19  tp.addTask(thrd_MatrixProduct, C, A, B, i, j); 20  } 21  } 22 //線程池啟動並阻塞等待
 tp.start(); 24  tp.join(); 25 }

要注意的是，threadPool和matrix類都不是STL，需要自行實現。其中threadPool可以用STL封裝的線程std::thread或者C語言的pthread庫來實現。std::thread類的聲明中delete了拷貝構造函數和賦值運算符重載，這一點請在實現階段注意。

事實上，這個操作在數據量較小的時候看不出優勢；而數據量過大時，線程調度又是操作系統的難題。我們需要可以容納更多線程的廉價計算資源。

cuda正是給顯卡計算這一廉價而高效的並行計算方式提供了接口，同時也不需要線程池的維護。比如上述問題，用cuda實現，或許過程有點復雜，但核心還是相當容易的：

template <typename T>
__global__ void kern_MatrixProduct(T* C_ptr, const T* A_ptr, const T* B_ptr,  3  size_t A_Height, size_t A_Width, size_t B_Width) {  4 //獲取當前線程調用者的Id
    size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;  6     if(Idx >= A_Height * B_Width) return;  7     
//規定計算第i行j列的線程Id為(i * width + j)
    size_t row = Idx / B_Width; 10     size_t col = Idx % B_Width; 11 //核心計算
    T ans = 0; 13     for(size_t k = 0; k < A_Width; ++k) { 14         ans += A_ptr[row * A_Width + k] * B_ptr[k * B_Width + col]; 15  } 16     C_ptr[row * B_Width + col] = ans; 17 } 18 
template <typename T>
void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) { 21     assert(A.Width() == B.Height()); 22  C.Set_size(A.Height(), B.Width()); 23     
    T *A_ptr, *B_ptr, *C_ptr; 25 //在顯卡中申請內存（顯存）
    assert(cudaSuccess == cudaMalloc(&A_ptr, sizeof(T) * A.Height() * A.Width())); 27     assert(cudaSuccess == cudaMalloc(&B_ptr, sizeof(T) * B.Height() * B.Width())); 28     assert(cudaSuccess == cudaMalloc(&C_ptr, sizeof(T) * A.Height() * B.Width())); 29 //將數據從內存復制到顯存
    assert(cudaSuccess == cudaMemcpy(A_ptr, A.pos_ptr(), sizeof(T) * A.Height() * A.Width(), cudaMemcpyHostToDevice)); 31     assert(cudaSuccess == cudaMemcpy(B_ptr, B.pos_ptr(), sizeof(T) * B.Height() * B.Width(), cudaMemcpyHostToDevice)); 32 //調用核函數
    size_t thread_count = 1024; 34     size_t block_count = (A.Height() * B.Width() - 1) / thread_count + 1; 35     kern_MatrixProduct<<<block_count, thread_count>>> (C_ptr, A_ptr, B_ptr, A.Height(), A.Width(), B.Width()); 36  cudaDeviceSynchronize(); 37     assert(cudaSuccess == cudaGetLastError()); 38 //將數據從顯存復制到內存
    assert(cudaSuccess == cudaMemcpy(C.pos_ptr(), C_ptr, sizeof(T) * C.Height() * C.Width(), cudaMemcpyDeviceToHost)); 40 //釋放臨時內存
    assert(cudaSuccess == cudaFree(A_ptr)); 42     assert(cudaSuccess == cudaFree(B_ptr)); 43     assert(cudaSuccess == cudaFree(C_ptr)); 44 }

當然，如果你在實現matrix的時候就把其元素數組申請在顯存中，這個代碼就會變得更加容易：

template <typename T>
__global__ void kern_MatrixProduct(T* C_ptr, const T* A_ptr, const T* B_ptr,  3  size_t A_Height, size_t A_Width, size_t B_Width) {  4 //獲取當前線程調用者的Id
    size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;  6     if(Idx >= A_Height * B_Width) return;  7     
//規定計算第i行j列的線程Id為(i * width + j)
    size_t row = Idx / B_Width; 10     size_t col = Idx % B_Width; 11 //核心計算
    T ans = 0; 13     for(size_t k = 0; k < A_Width; ++k) { 14         ans += A_ptr[row * A_Width + k] * B_ptr[k * B_Width + col]; 15  } 16     C_ptr[row * B_Width + col] = ans; 17 } 18 
template <typename T>
void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) { 21     assert(A.Width() == B.Height()); 22  C.Set_size(A.Height(), B.Width()); 23 //調用核函數
    size_t thread_count = 1024; 25     size_t block_count = (A.Height() * B.Width() - 1) / thread_count + 1; 26     kern_MatrixProduct<<<block_count, thread_count>>> (C.pos_ptr(), A.pos_ptr(), B.pos_ptr(), A.Height(), A.Width(), B.Width()); 27  cudaDeviceSynchronize(); 28     assert(cudaSuccess == cudaGetLastError()); 29 }

最后一種寫法效率是最高的——第一種沒有使用並行，第二種的並行優勢不明顯（而且申請線程太多，消耗大量時間），第三種I/O的耗時過大。

Windows10系統下的實驗表明，當矩陣長寬小於10時，1、2、4三種寫法的函數都能在0.01秒內完成，第三種也能在0.02秒內計算完畢。而當矩陣長寬為512時，第一種在1.1s左右計算完畢，第二種10s內計算不出結果，第三種則在1.3s左右結束，而第四種只有0.03s。

如果說即使如此也看不出1、3、4三種方式的孰優孰劣的話，試想一下，如果是若干個矩陣連乘，效果會是如何？結果不言而喻：第三種最慢，因為頻繁調用I/O；第一種次之，因為沒有並行；第四種的優勢將會極其明顯，因為I/O少、並行效率高。

很多人好奇，具體是如何並行的呢？上文代碼中的三層尖括號（<<< >>>）是什么？莫急，讓我一一解釋。

首先我們來了解顯卡的結構。

顯卡的簡化結構

顯卡內部，有三級結構：網格（grid）、塊（block）、線程（thread）。

每個顯卡只有很少的網格，一個核函數目前只能運行在一個網格中，而一個網格里有多個塊，每個塊包含了若干線程。具體包含了多少，可以使用上文圖中的deviceQuery來查詢，通常是256、1024這樣的值。

 同樣是這個矩陣相乘的例子：

 上文的核函數在執行時，選取了一個網格中的一個塊進行並行計算：

 執行核函數時GPU一個執行塊的簡化示意圖，大多數線程都被直接return掉

關於顯卡的介紹到此為止，太過深入的知識我將不作介紹。

如何將計算部署在這些塊和線程之上呢？這就是核函數（kernel function）及其運行時參數（runtime parameter）的工作。

核函數及其運行時參數

相信大家也猜到了，核函數就是上文代碼里有__global__修飾的函數，其運行時參數就是放在三重尖括號<<< >>>之中的值。

每次調用核函數，都要指定其運行時參數，以便在顯卡中彈性地部署計算。而其函數參數和我們學過的C/C++寫法是一樣的，不過要寫在運行時參數的后面。

運行時參數有兩種寫法：

整數的寫法：第一個參數是部署的核數，第二個參數是每個核分配的線程數，即在調用時：

kernelFunc<<<b, t>>> (arg1, arg2);

　　注意不能超過每個核的最大線程數量，否則cudaGetLastError會返回一個

　　cudaErrorInvalidConfiguration錯誤。運行時通過blockIdx.x * blockDim.x + threadIdx.x來獲取線程ID，其區間范圍為。如果的值大於你的預期，請在不必要的線程中及時return，以免出現數組越界、訪問無效內存等情況。有人可能會問，我挨個指定每個核開辟多少線程不行嗎？不行。核函數的線程，往往都是要稍微開多一些，然后return掉不必要的。英偉達就是這么霸道，我不要你覺得，我要我覺得。

　　三維向量的寫法：第一個參數指定了網格中的維度，即如何分配塊；第二個參數指定了每個塊的維度，即如何分配線程；第三個參數可選，指定了核函數運行時每個block除了使用共享顯存外，還能再動態申請多少空間，默認是0，單位是字節——這和malloc函數參數的意義一樣；第四個參數也是可選的，指定了該核函數運行在哪一個流中，默認是0——這個流方便主顯異步並行，但相對麻煩，暫時不作介紹。

　　其中的類型是dim3，即三維無符號整數向量：x坐標表示該網格中沿x坐標的核數，即每行的核數；y坐標表示該網格中沿y坐標的核數，即每列的核數；z坐標只能為1，因為網格中每一個核都是平面擺放的。因此核函數部署的核數就是。

　　而的類型也是dim3：x坐標表示一個核中沿x坐標的線程數，即每行的線程數；y坐標表示一個核中沿y坐標的線程數，即每列的線程數；z坐標表示一個核中沿z坐標的線程數，即每個核的線程有多少個水平面。因此核函數在每個核開辟的線程數就是。線程ID的值的取值范圍則在區間中。

　　而這些值都是有最大值的，同樣也可以通過deviceQuery查詢，我的顯卡中五個參數的最大值分別為65536、65536、1024、1024、64。

　　而核函數運行在顯卡中，這就限制了它訪問內存中的變量——但我們可以將內存中的變量傳遞到顯存中。下一篇文章將會介紹主存與顯存的交互。