CUDA編程指南閱讀筆記

轉自http://blog.csdn.net/csgxy123/article/category/1500471

隨着多核CPU和衆核GPU的到來，並行編程已經獲得了業界愈來愈多的重視，CPU-GPU異構程序可以極大提升現有計算機系統的運算性能，對於科學計算等運算密集型程序有着很是重要的意義。這一系列文章是根據《CUDA C語言編程指南》來整理的，該指南是NVIDIA公司提供的CUDA學習資料，介紹了CUDA編程最基本最核心的概念，是學習CUDA必不可少的閱讀材料。

初學CUDA，筆記錯誤之處在所不免，還請發現問題的諸位讀者不吝賜教。

1. 什麼是CUDA？

      CUDA全稱是Compute Unified Device Architecture，中文名稱即統一計算設備架構，它是NVIDIA公司提出了一種通用的並行計算平臺和編程模型。使用CUDA，咱們能夠開發出同時在CPU和GPU上運行的通用計算程序，更加高效地利用現有硬件進行計算。爲了簡化並行計算學習，CUDA爲程序員提供了一個類C語言的開發環境以及一些其它的如FORTRAN、DirectCOmpute、OpenACC的高級語言/編程接口來開發CUDA程序。

2. CUDA編程模型如何擴展？

     咱們知道，不一樣的GPU擁有不一樣的核心數目，在覈心較多的系統上CUDA程序運行的時間較短，而在覈心較少的系統上CUDA程序的執行時間較多。那麼，CUDA是如何作到的呢？

      並行編程的中心思想是分而治之：將大問題劃分爲一些小問題，再把這些小問題交給相應的處理單元並行地進行處理。在CUDA中，這一思想便體如今它的具備兩個層次的問題劃分模型。一個問題能夠首先被粗粒度地劃分爲若干較小的子問題，CUDA使用被稱爲塊（Block）的單元來處理它們，每一個塊都由一些CUDA線程組成，線程是CUDA中最小的處理單元，將這些較小的子問題進一步劃分爲若干更小的細粒度的問題，咱們即可以使用線程來解決這些問題了。對於一個普通的NVIDIA GPU，其CUDA線程數目一般能達到數千個甚至更多，所以，這樣的問題劃分模型即可以成倍地提高計算機的運算性能。

      GPU是由多個流水多處理器構成的，流水處理器以塊（Block）爲基本調度單元，所以，對於流水處理器較多的GPU，它一次能夠處理的塊（Block）更多，從而運算速度更快，時間更短。而反之對於流水處理器較少的GPU，其運算速度便會較慢。這一原理能夠經過下圖形象地看出來：

3. CUDA基本概念（上）

      本節將介紹CUDA的一些基本的編程概念，該節用到的例子來自於CUDA Sample中的VectorAdd項目。

3.1 內核（Kernels）

      CUDA C是C語言的一個擴展，它容許程序員定義一種被稱爲內核函數（Kernel Functions）的C函數，內核函數運行在GPU上，一旦啓動，CUDA中的每個線程都將會同時並行地執行內核函數中的代碼。

      內核函數使用關鍵字__global__來聲明，運行該函數的CUDA線程數則經過<<<...>>>執行配置語法來設置。（參見章節"C語言擴展"），每個執行內核函數的線程都由一個惟一的線程ID，這一ID能夠經過在內核函數中訪問threadIdx變量來獲得。

      下面經過一些示例代碼來展現剛剛提到的這些概念該如何應用在編程中：

 // Kernel definition
 __global__ void VecAdd(float* A, float* B, float* C) {
     int i = threadIdx.x;
     C[i] = A[i] + B[i];
 }

 int main() {
     ...
     // Kernel invocation with N threads
     VecAdd<<<1, N>>>(A, B, C);
    ...
 }

在上面的代碼中，N個線程將會並行地同時執行加法運算。

3.2 線程層次（Thread Hierarchy）

      CUDA的每個線程都有其線程ID，線程的ID信息由變量threadIdx給出。threadIdx是CUDA C語言的內建變量，一般它用一個三維數組來表示。使用三維數組的方便之處在於能夠很方便地表示一維、二維和三維線程索引，進而方便地表示一維、二維和三維線程塊（thread block）。這樣，不管是數組、矩陣仍是體積的計算，均可以很容易地使用CUDA進行運算。

      線程的索引與線程ID之間存在着直接的換算關係，對於一個索引爲(x, y, z)的線程來講：

      一、若是線程塊（block）是一維的，則線程ID = x

      二、若是線程塊是二維的，假設塊尺寸爲（Dx，Dy），那麼線程ID = x + y * Dx

      三、若是線程塊是三維的，設其尺寸爲（Dx，Dy，Dz），那麼線程ID = x + y * Dx + z * Dx * Dy

      下面的例子展現了兩個NxN矩陣相加的CUDA實現：

 // Kernel definition
 __global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {
     int i = threadIdx.x;
     int j = threadIdx.y;
     C[i][j] = A[i][j] + B[i][j];
 }

 int main() {
     ...
     // Kernel invocation with one block of N * N * 1 threads
     int numBlocks = 1;
     dim3 threadsPerBlock(N, N);
     MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
     ...
 }

     每一個線程塊（block）中的線程數量是有限制的，由於依據前面所說，同一線程塊（block）中的全部線程都會被分配到同一個處理器核上運行，共享有限的存儲資源，所以對於當前的GPU，線程塊所能包含的最大線程數目爲1024。

      上面的例子中numBlocks表明線程塊的數量，這裏的值爲1。在通常的CUDA程序中，這個值一般大於1，也就是說將會有多個線程塊被分配到多個處理器核中同時進行處理，這樣就大大提升了程序的並行性。

      在CUDA中，線程塊包含在線程格（grid）當中，線程格能夠是一維、二維或者三維的，線程格的尺寸通常根據待處理數據的規模或者處理器的數量來指定。線程格中所包含的線程塊數目一般遠遠大於GPU處理器核心的數目。下圖展現了線程格（grid）、線程塊（block）以及線程（thread）之間的關係：

      內核函數的調用能夠簡化爲kernel<<<A，B>>>(parameters)，在尖括號中，A表明線程格（grid）的尺寸，它能夠是三維的，用類型dim3表示，也能夠是一維的，用int類型表示。B表明線程塊（block）的尺寸，它與A相似，也可分別用dim3或int類型表示。

      在內核函數內部，CUDA爲咱們內建了一些變量用於訪問線程格、線程塊的尺寸和索引等信息，它們是：

      1. gridDim：表明線程格（grid）的尺寸，gridDim.x爲x軸尺寸，gridDim.y、gridDim.z相似。拿上圖來講，它的gridDim.x = 3，gridDim.y = 2，gridDim.z = 1。

      2. blockIdx：表明線程塊（block）在線程格（grid）中的索引值，拿上圖來講，Block（1,1）的索引值爲：blockIdx.x = 1，blockIdx.y = 1。

      3. blockDim：表明線程塊（block）的尺寸，blockDIm.x爲x軸尺寸，其它依此類推。拿上圖來講，注意到Block（1,1）包含了4 * 3個線程，所以blockDim.x = 4, blockDim.y = 3。

      4. threadIdx：線程索引，前面章節已經詳細探討過了，這裏再也不贅述。

      明白了這些變量的含義，那麼下面的矩陣加法程序便不難理解了：

 // Kernel definition
 __global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {
     int i = blockIdx.x * blockDim.x + threadIdx.x;
     int j = blockIdx.y * blockDim.y + threadIdx.y;
     if (i < N && j < N)
         C[i][j] = A[i][j] + B[i][j];
 }

 int main() {
     ...
     // Kernel invocation
     dim3 threadsPerBlock(16, 16);
     dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
     MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
     ...
 }

      在上面的程序中，線程塊（block）的尺寸是16x16，這是CUDA編程中一個很是廣泛的選擇。線程格（grid）包含了足夠多的線程塊（block）來進行計算。

      線程塊（block）是獨立執行的，在執行的過程當中線程塊之間互不干擾，所以它們的執行順序是隨機的。

      同一線程塊中的線程能夠經過訪問共享內存（shared memory）或者經過同步函數__syncthreads()來協調合做。這些概念將在之後的章節中詳細解釋。

3. CUDA基本概念（下）

3.3 內存層次（Memory Hierarchy）

      在GPU上CUDA線程能夠訪問到的存儲資源有不少，每一個CUDA線程擁有獨立的本地內存（local Memory）；每個線程塊（block）都有其獨立的共享內存（shared memory），共享內存對於線程塊中的每一個線程都是可見的，它與線程塊具備相同的生存時間；同時，還有一片稱爲全局內存（global memory）的區域對全部的CUDA線程都是可訪問的。

      除了上述三種存儲資源之外，CUDA還提供了兩種只讀內存空間：常量內存（constant memory）和紋理內存（texture memory），同全局內存相似，全部的CUDA線程均可以訪問它們。對於一些特殊格式的數據，紋理內存提供多種尋址模式以及數據過濾方法來操做內存。這兩類存儲資源主要用於一些特殊的內存使用場合。

      一個程序啓動內核函數之後，全局內存、常量內存以及紋理內存將會一直存在直到該程序結束。下面是CUDA的內存層次圖：

3.4 異構編程（Heterogeneous Programming）

      CUDA的異構編程模型假定CUDA線程都運行在一個可被看作CPU協處理器的芯片上，這就使得CUDA內核函數能夠和CPU端C程序的運行並行運行，從而加快程序的運行效率。爲了達到這個效果，CUDA程序須要管理兩大塊由DRAM構成的內存區域：CPU端能夠訪問到的主機內存（host memory）以及GPU端供CUDA內核訪問到的設備內存（device memory），設備內存主要由全局內存、常量內存以及紋理內存構成。如今，CUDA程序的運行機制便很明瞭了：CPU端代碼生成原始數據，經過CUDA運行時函數庫將這些原始數據傳輸到GPU上，在CPU端啓動CUDA內核函數進行運算，而後將運算結果從設備端傳輸到主機端，計算任務便完成了。

4. CUDA C語言編程接口

    異構程序設計跟傳統的串行程序設計差異是很大的，學習起來也是很是不容易的。NVIDIA很是夠意思，爲了簡化CUDA的學習曲線，它採用了絕大多數程序員都熟悉的C語言做爲其根基，CUDA C是NVIDIA爲程序員提供的一類編程接口，它其實是一個C語言的擴展，在C的基礎上增長了一些新的語法和變量，而且提供了功能豐富的庫函數，方便程序員使用GPU進行異構計算。

        除了前面章節提到的CUDA最基本、最核心的概念之外，CUDA C呈現給程序員的接口主要由兩大類API構成，它們分別是CUDA Runtime API和CUDA Driver API，Runtime API其實是對於Driver API的封裝，其目的天然是方便程序員的代碼編寫工做。Driver API爲用戶提供了更細一層的控制手段，經過它能夠控制諸如CUDA Contexts（一種相似主機進程的概念）以及CUDA Modules（相似主機動態加載庫的概念）等更加底層的CUDA模塊。

4.1 NVCC編譯器

        任何一種程序設計語言都須要相應的編譯器將其編譯爲二進制代碼，進而在目標機器上獲得執行。對於異構計算而言，這一過程與傳統程序設計語言是有一些區別的。爲何？由於CUDA它本質上不是一種語言，而是一種異構計算的編程模型，使用CUDA C寫出的代碼須要在兩種體系結構徹底不一樣的設備上執行：一、CPU；二、GPU。所以，CUDA C的編譯器所作的工做就有點略多了。一方面，它須要將源代碼中運行在GPU端的代碼編譯獲得能在CUDA設備上運行的二進制程序。另外一方面，它也須要將源代碼中運行在CPU端的程序編譯獲得能在主機CPU上運行的二進制程序。最後，它須要把這兩部分有機地結合起來，使得兩部分代碼可以協調運行。

        CUDA C爲咱們提供了這樣的編譯器，它即是NVCC。嚴格意義上來說，NVCC並不能稱做編譯器，NVIDIA稱其爲編譯器驅動（Compiler Driver），本節咱們暫且使用編譯器來描述NVCC。使用nvcc命令行工具咱們能夠簡化CUDA程序的編譯過程，NVCC編譯器的工做過程主要能夠劃分爲兩個階段：離線編譯（Offline Compilation）和即時編譯（Just-in-Time Compilation）。

離線編譯（Offline Compilation）

    下面這幅圖簡單說明了離線編譯的過程：

        在CUDA源代碼中，既包含在GPU設備上執行的代碼，也包括在主機CPU上執行的代碼。所以，NVCC的第一步工做即是將兩者分離開來，這一過程結束以後：

        1. 運行於設備端的代碼將被NVCC工具編譯爲PTX代碼（GPU的彙編代碼）或者cubin對象（二進制GPU代碼）；

        2. 運行於主機端的代碼將被NVCC工具改寫，將其中的內核啓動語法（如<<<...>>>）改寫爲一系列的CUDA Runtime函數，並利用外部編譯工具（gcc for linux，或者vc compiler for windows）來編譯這部分代碼，以獲得運行於CPU上的可執行程序。

        完事以後，NVCC將自動把輸出的兩個二進制文件連接起來，獲得異構程序的二進制代碼。

即時編譯（Just-in-time Compile）

        任何在運行時被CUDA程序加載的PTX代碼都會被顯卡的驅動程序進一步編譯成設備相關的二進制可執行代碼。這一過程被稱做即時編譯（just-in-time compilation）。即時編譯增長了程序的裝載時間，可是也使得編譯好的程序能夠重新的顯卡驅動中得到性能提高。同時到目前爲止，這一方法是保證編譯好的程序在還未問世的GPU上運行的惟一解決方案。

        在即時編譯的過程當中，顯卡驅動將會自動緩存PTX代碼的編譯結果，以免屢次調用同一程序帶來的重複編譯開銷。NVIDIA把這部分緩存稱做計算緩存（compute cache），當顯卡驅動升級時，這部分緩存將會自動清空，以使得程序可以自動得到新驅動爲即時編譯過程帶來的性能提高。

        有一些環境變量能夠用來控制即時編譯過程：

        1. 設置CUDA_CACHE_DISABLE爲1將會關閉緩存功能

        2. CUDA_CACHE_MAXSIZE變量用於指定計算緩存的字節大小，默認狀況下它的值是32MB，它最大能夠被設置爲4GB。任何大於緩存最大值得二進制代碼將不會被緩存。在須要的狀況下，一些舊的二進制代碼可能被丟棄以騰出空間緩存新的二進制代碼。

        3. CUDA_CACHE_PATH變量用於指定計算緩存的存儲目錄地址，它的缺省值以下：

        4. 設置CUDA_FORCE_PTX_JIT爲1會強制顯卡驅動忽略應用程序中的二進制代碼而且即時編譯程序中的嵌入PTX代碼。若是一個內核函數沒有嵌入的PTX代碼，那麼它將會裝載失敗。該變量能夠用來確認程序中存在嵌入的PTX代碼。同時，使用即時編譯（just-in-time Compilation）技術也可確保程序的向前兼容性。

4.2 兼容性

      一、二進制兼容性

      二進制代碼是設備相關的，使用NVCC編譯器編譯時，若指定-code選項，則會編譯產生目標設備的二進制cubin對象。例如，編譯時使用-code=sm_13會產生適用於計算能力1.3的二進制代碼。二進制代碼在CUDA計算設備上具備小版本的向前兼容性，可是在大版本上不具有兼容性。也就是說，對於計算能力X.y的硬件，使用-code=sm_Xy編譯後，程序可以運行於計算能力X.z（其中z>=y）的硬件上，但不能運行在計算能力M.n（M!=X）的硬件上。

      二、PTX代碼兼容性

      不一樣計算能力的設備所支持的PTX指令條數是不一樣的，一些PTX指令只在擁有較高計算能力的設備上被支持。例如，全局內存（global Memory）的原子操做指令只能用於計算能力不小於1.1的設備；雙精度浮點運算指令只能用於計算能力不小於1.3的設備。在將C語言編譯爲PTX代碼時，NVCC使用-arch編譯選項指定PTX代碼目標設備的計算能力。所以，要想使用雙精度運算，編譯時必須使用選項-arch=sm_13（或使用更高的計算能力），不然NVCC會自動將雙精度操做降級爲單精度操做。

      爲某一特定設備產生的PTX代碼，在運行時老是可以被具備更高計算能力的設備JIT編譯爲可執行的二進制代碼。

      三、應用程序兼容性

      執行CUDA程序有兩種方式，一種是直接加載編譯好的CUDA二進制代碼運行，另外一種是首先加載程序中的PTX代碼，再執行JIT編譯獲得二進制的設備可執行文件，而後運行。特別須要注意的是，爲了讓程序運行具備更高計算能力的將來設備上，必須讓程序加載PTX代碼。

      事實上，在一個CUDA C程序中能夠嵌入不止一個版本的PTX/二進制代碼。那麼，具體執行時哪個版本的PTX或者二進制代碼會獲得執行呢？答案是：最兼容的那個版本。例如編譯一個名爲x.cu的CUDA源代碼：

將會產生兼容計算能力1.1硬件的二進制代碼（第一排的-gencode選項）以及兼容計算能力1.1設備的PTX和二進制代碼，這些代碼都將會嵌入到編譯後的目標文件中。

      主機端將會產生一些額外的代碼，在程序運行時，這些代碼會自動決定裝載哪個版本的代碼來執行。對於上面的例子：

• 計算能力1.0的設備運行該程序將會裝載1.0版本的二進制代碼
• 計算能力1.一、1.2或者1.3的設備運行該程序將會裝載1.1版本的二進制代碼
• 計算能力2.0或者更高的設備運行該程序將會裝載1.1版本的PTX代碼進而對其進行JIT編譯獲得相應設備的二進制代碼

      同時，x.cu還能夠在程序中使用一些特殊的宏來改變不一樣設備的代碼執行路徑。例如，對於計算能力1.1的設備而言，宏__CUDA_ARCH__等於110，在程序中能夠對該宏的值進行判斷，而後分支執行程序。

      NVCC用戶手冊列出了不少-arch，-code和-gencode等編譯選項的簡化書寫形式。例如，-arch=sm_13就是-arch=compute_13 -code=compute13, sm_13的簡化形式。更多詳盡的內容請參閱該手冊。

      四、C/C++兼容性

      NVCC編譯器前端使用C++語法啊規則來處理CUDA源文件。在主機端，CUDA支持完整的C++語法；而在設備端，只有部分C++語法是被支持的。這方面更爲詳盡的討論請參見《CUDA C程序設計指南》的C/C++語言支持章節。

      五、64位兼容性

      64位版本的nvcc編譯器將設備代碼編譯爲64位模式，即指針是64位的。運行64位設備代碼的先決條件是主機端代碼必須也使用64位模式進行編譯。一樣，32位版本的nvcc將設備代碼編譯爲32位模式，這些代碼也必須與相應的32位主機端代碼相配合方能運行。

      32位nvcc編譯器可使用-m64編譯選項將設備代碼編譯爲64位模式。同時64位nvcc編譯器也可以使用-m32編譯選項將設備代碼編譯爲32位模式。

4.3 CUDA C Runtime

        CUDA C Runtime使用cudart動態連接庫實現（cudart.dll或者cudart.so），運行時中全部的入口函數都以cuda爲前綴。

4.3.1 初始化

        CUDA C Runtime函數庫沒有明確的初始化函數，在程序第一次調用Runtime庫函數時它會自動初始化。所以，在記錄Runtime函數調用時間和理解程序中第一個Runtime調用返回的錯誤代碼時，須要將初始化考慮在內。

        在初始化期間，Runtime將會爲系統中每個設備建立一個CUDA上下文（相似CPU中進程的數據結構），這個上下文是設備的基本上下文，它被程序中全部的主機線程所共享。建立過程在後臺運行，而且，Runtime將隱藏基本上下文使之對Runtime API這一層次的程序員不可見。

        當一個主機線程調用cudaDeviceReset()函數時，它將會銷燬線程當前控制設備的基本上下文。也就是說，當線程下一次調用runtime函數時將會重啓初始化，一個新的CUDA基本上下文將被建立出來。

4.3.2 設備內存

        正如前面異構計算章節所講，CUDA編程模型假定系統是由主機和設備構成的，它們分別具備本身獨立的內存空間。Runtime負責設備內存的分配，回收，拷貝以及在主機和設備間傳輸數據的工做。

        設備內存能夠有兩種分配方式：線性內存或者CUDA數組

        CUDA數組是一塊不透明的內存空間，它主要被優化用於紋理存取。

        線性內存空間與平時咱們訪問的內存相似，對於計算能力1.x的設備來講，它存在於一個32位的地址空間。對於更高計算能力的設備而言，它存在於一個40位的地址空間中。所以，單獨分配的實體可使用指針來相互應用。

        咱們一般使用cudaMalloc()函數分配線性內存空間，使用cudaFree()函數釋放線性內存空間，使用cudaMemcpy()函數在主機和設備之間傳輸數據。下面是CUDA Vector Add代碼示例的一些片斷：

 // Device code
 __global__ void VecAdd(float *A, float *B, float *C, int N) {
     int i = blockDim.x * blockIdx.x + threadIdx.x;
     if (i < N)
         C[i] = A[i] + B[i];
 }

 // Host code
 int main() {
     int N = ...;
     size_t size = N * sizeof(float);

     // Allocate input vectors h_A and h_B in host memory
     float *h_A = (float*)malloc(size);
     float *h_B = (float*)malloc(size);

     // Initialize input vectors
     ...

     // Allocate vectors in device memory
     float *d_A, *d_B, *d_C;
     cudaMalloc(&d_A, size);
     cudaMalloc(&d_B, size);
     cudaMalloc(&d_C, size);

     // Copy vectors from host memory to device memory
     cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
     cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

     // Invoke kernel
     int threadsPerBlock = 256;
     int blocksPerGrid = (N +threadsPerBlock - 1) / threadsPerBlock;
     VecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

     // Copy result from device memory to host Memory
     cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

     // Free device memory
     cudaFree(d_A);
     cudaFree(d_B);
     cudaFree(d_C);

     // Free host memory
     ...
 }

片斷展現了設備內存的分配，傳輸以及回收過程。

        除了上面展現的方法，咱們還可使用cudaMallocPitch()和cudaMalloc3D()函數來分配線性內存。這些函數可以確保分配的內存知足設備內存訪問的對齊要求，對於行地址的訪問以及多維數組間的數據傳輸提供高性能保證，所以很是適合對於二維和三維數組內存空間的分配。下面的代碼片斷展現了分配和使用尺寸爲width x height的二維數組的技術：

 // Host code
 int width = 64, height = 64;
 float *devPtr;
 size_t pitch;
 cudaMallocPitch(&devPtr, &pitch, width * sizeof(float), height);
 MyKernel<<<100, 512>>>(devPtr, pitch, width, height);

 // Device code
 __global__ void MyKernel(float* devPtr, size_t pitch, int width, int height) {
     for (int r = 0; r < height; ++r) {
         float* row = (float*)((char*)devPtr + r * pitch);
         for (int c = 0; c < width; ++c) {
             float element = row[c];
         }
     }
 }

        下面的代碼片斷展現了一個尺寸爲width x height x depth的三維數組的分配和使用方法：

 // Host code
 int width = 64, height = 64, depth = 64;
 cudaExtent extent = make_cudaExtent(width * sizeof(float), height, depth);
 cudaPitchedPtr devPitchedPtr;
 cudaMalloc3D(&devPitchedPtr, extent);
 MyKernel<<<100, 512>>>(devPitchedPtr, width, height, depth);

 // Device code
 __global__ void MyKernel(cudaPitchedPtr devPitchedPtr, int width, int height, int depth) {
     char* devPtr = devPitchedPtr.ptr;
     size_t pitch = devPitchedPtr.pitch;
     size_t slicePitch = pitch * height;
     for (int z = 0; z < depth; ++z) {
         char* slice = devPtr + z * slicePitch;
         for (int y = 0; y < height; ++y) {
             float* row = (float*)(slice + y * pitch);
             for (int x = 0; x < width; ++x)
                 float element = row[x];
         }
     }
 }

更多詳細的內容請查閱參考手冊。

        下面的代碼示例展現了多種使用Runtime API訪問全局變量的技術：

 __constant__ float constData[256];
 float data[256];
 cudaMemcpyToSymbol(constData, data, sizeof(data));
 cudaMemcpyFromSymbol(data, constData, sizeof(data));

 __device__ float devData;
 float value = 3.14f;
 cudaMemcpyToSymbol(devData, &value, sizeof(float));

 __device__ float* devPointer;
 float* ptr;
 cudaMalloc(&ptr, 256 * sizeof(float));
 cudaMemcpyToSymbol(devPointer, &ptr, sizeof(ptr));

使用cudaGetSymbolAddress()函數能夠得到被聲明存儲在全局內存中的變量地址。爲了得到分配內存的大小，可使用cudaGetSymbolSize()函數。

4.3.3 共享內存（Shared Memory）

        共享內存是CUDA設備中很是重要的一個存儲區域，有效地使用共享內存能夠充分利用CUDA設備的潛能，極大提高程序性能。那麼，共享內存有哪些特色呢？

        一、共享內存（shared Memory）是集成在GPU處理器芯片上的（on-chip），所以相比於存在於顯存顆粒中的全局內存（global Memory）和本地內存（local Memory），它具備更高的傳輸帶寬，通常狀況下，共享內存的帶寬大約是全局內存帶寬的7-10倍。

        二、共享內存的容量很小。根據NVIDIA官方文檔的說法，在計算能力1.x的設備中，每個流多處理器（Streaming Multiprocessor）上的共享內存容量爲16KB。對於計算能力2.x、3.0及3.5的設備該參數爲48KB。所以共享內存是稀有資源。

        三、共享內存在物理上被劃分爲不少塊，每一塊被稱爲一個存儲體（bank）。在同一時刻，CUDA設備能夠同時訪問多個存儲體。所以，若是一次針對共享內存的訪存操做須要讀取n個地址，而這n個地址剛好分佈在n個不一樣的存儲體（bank）中，那麼只須要一個存取週期就能夠完成n個地址的訪存任務了。對於計算能力1.x的設備，共享內存被平均劃分爲16個存儲體。而對於計算能力2.x、3.0及3.5的設備此參數爲32。在共享內存中，相鄰兩塊32bit的數據分別屬於相鄰的兩個存儲體。存儲體每兩個時鐘週期能夠傳輸32位數據。

        四、共享內存既能夠靜態分配，也能夠動態分配。

        從共享內存的這些特色中咱們能夠看出，它實際上至關於一個程序員能夠操控的緩存（cache），下面，咱們使用矩陣乘法的例子來講明如何有效使用共享內存。

        首先，咱們使用最直觀的方法來完成矩陣乘法C = A x B：讀取A的每一行和B的每一列，順次完成計算任務。矩陣乘法的示意圖以下所示：

下面是矩陣乘法的CUDA C主要實現代碼：

 // 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_SIZE
 void 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, d_c.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.xl
     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;
 }

能夠看出，爲了計算矩陣C的任何一個元素，程序都須要從全局內存（global memory）中得到矩陣A的一行和矩陣B的一列。所以，完成這一計算矩陣A被讀取了B.width次，矩陣B被讀取了A.height次。
        如今咱們來使用共享內存（shared memory）實現矩陣乘法。假設矩陣C能夠被劃分爲若干個較小的子方陣Csub，咱們使用一個線程塊（thread block）來負責某一子方陣的計算，線程塊中的每個線程（thread）正好負責子方陣Csub中一個元素的計算。這樣劃分後，任何一個結果子方陣Csub'（尺寸爲block_size * block_size）都是與該方陣具備相同行索引的尺寸爲A.width * block_size的A的子矩陣Asub和與該方陣具備相同列索引的尺寸爲block_size * B.height的B的子矩陣Bsub相乘所獲得。

        爲了匹配設備的計算資源，兩個子矩陣Asub和Bsub被劃分爲儘量多的分離的維度爲block_size的子方陣，Csub的值即是這些子矩陣相乘後相加所獲得的結果。子矩陣乘法的執行順序都是首先將它們從全局內存（global memory）拷貝到共享內存（shared memory）（線程塊中的每個線程正好負責方陣一個元素的拷貝），而後由線程本身完成相應元素的計算任務，利用寄存器存儲局部結果，最後將寄存器的內容與新獲得的計算結果依此累加起來獲得最終運算結果並將其傳輸到全局內存（global memory）中。

        經過使用這種分治的計算策略，共享內存獲得了很好的利用，採用這種方案計算完成時全局內存中矩陣A被訪問的次數爲B.width / block_size，矩陣B被訪問的次數爲A.height / block_size，很明顯，這爲咱們節省了很是多的全局內存帶寬。優化後的矩陣計算示意圖以下所示：

        爲了提高計算效率，咱們爲類型Matrix增長了一個成員變量stride。__device__函數用來得到和設置子矩陣的元素。下面是優化後的代碼：

 // Matrices are stored in row-major order;
 // M(row, col) = *(M.elements + row * M.stride + col)
 typedef struct {
     int width;
     int height;
     int stride;
     float* elements;
 } Matrix;

 // Get a matrix element
 __device__ float GetElement(const Matrix A, int row, int col) {
     return A.elements[row * A.stride + col];
 }

 // Set a matrix element
 __device__ void SetElement(Matrix A, int row, int col, float value) {
     A.elements[row * A.stride + col] = value;
 }

 // Get the BLOCK_SIZExBLOCK_SIZE sub-matrix Asub of A that is
 // located col sub-matrices to the right and row sub-matrices down
 // from the upper-left corner of A
 __device__ Matrix GetSubMatrix(Matrix A, int row, int col) {
     Matrix Asub;
     Asub.width = BLOCK_SIZE;
     Asub.height = BLOCK_SIZE;
     Asub.stride = A.stride;
     Asub.elements = &A.elements[A.stride * BLOCK_SIZE * row + BLOCK_SIZE * col];
     return Asub;
 }

 // 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_SIZE
 void MatMul(const Matrix A, const Matrix B, Matrix C) {
     // Load A and B to device memory
     Matrix d_A;
     d_A.width = d_A.stride = 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 = d_B.stride = 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 = d_C.stride = 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, d_C.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) {
     // 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;

     // Look 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);
 }

異步並行執行

主機和設備間並行執行

        爲了支持主機和設備的並行執行，CUDA提供了一些異步函數。異步是指設備在完成功能執行以前就將控制權交還給主機線程，以便主機線程繼續執行。這些函數有：

        一、內涵啓動（Kernel Launches）；

        二、同一設備內存中兩個地址塊之間的數據傳輸；

        三、從主機到設備的小於或等於64KB的一次數據塊傳輸；

        四、使用Async前綴函數進行的數據傳輸；

        五、內存置位函數調用（Memory set functions calls）。

程序員能夠經過設置環境變量CUDA_LAUNCH_BLOCKING來開啓或者關閉內核啓動（Kernel Launch）的異步功能。可是這一方法僅限於調試，在任何產品代碼中不該當關閉異步內核啓動。

        內核啓動在下面這些狀況下則是同步的：

        一、應用程序經過調試器或者內存檢查器運行在計算能力爲1.x的設備上。

        二、硬件計數器信息正被性能分析器收集。

將內核啓動與數據傳輸重疊起來

        對於一些計算能力等於或高於1.1的設備，它們能夠將內核啓動任務和鎖頁內存到設備內存的數據傳輸任務並行執行。應用程序能夠檢查設備屬性中的asyncEngineCount項來肯定設備是否支持這一功能。當該項值大於0時表明設備支持這一層次的並行。對於計算能力1.x的設備，該功能不支持經過cudaMallocPitch()函數分配的CUDA數組或2D數組。

並行內核執行

        一些計算能力2.x或更高的設備能夠同時並行執行多個內核函數。應用程序能夠檢查設備屬性中的concurrentKernels項來肯定設備是否支持這一功能，值爲1表明支持。運算能力3.5的設備在同一時刻可以並行執行的最大內核函數數量爲32，運算能力小於3.5的硬件則最多支持同時啓動16個內核函數的執行。同時須要注意的是，在一個CUDA上下文中的內核函數不能與另外一個CUDA上下文中的內核函數同時執行。使用不少紋理內存或者大量本地內存的內核函數也極可能沒法與其它內核函數並行執行。

並行數據傳輸

        一些計算能力爲2.x或更高的設備能夠將鎖頁內存到設備內存的數據傳輸和設備內存到鎖頁內存的數據傳輸並行執行。應用程序可檢查設備屬性中的asyncEngineCount項來肯定這一功能的支持程度，等於2時表示支持。

流（Streams）

        應用程序經過流來管理並行。一個流是一個順次執行的命令序列。不一樣的流之間並行執行，沒有固定的執行順序。

一、流的建立與銷燬

        定義一個流的過程一般包括：建立一個流對象，而後指定它爲內核啓動或者主機設備間數據傳輸的流參數。下面的一段代碼建立了兩個流而且在鎖頁內存中分配了一塊float類型的數組hostPtr：

 cudaStream_t stream[2];
 for (int i = 0; i < 2; ++i)
     cudaStreamCreate(&stream[i]);
 float *hostPtr;
 cudaMallocHost(&hostPtr, 2 * size);

下面的代碼定義了每個流的行爲：從主機端拷貝數據到設備端，內核啓動，從設備端拷貝數據到主機端：

 for (int i = 0; i < 2; ++i) {
     cudaMemcpyAsync(inputDevPtr + i * size, hostPtr + i * size, size, cudaMemcpyHostToDevice, stream[i]);
     MyKernel<<<100, 512, 0, stream[i]>>>(outputDevPtr + i * size, inputDevPtr + i * size, size);
     cudaMemcpyAsync(hostPtr + i * size, outputDevPtr + i * size, size, cudaMemcpyDeviceToHost, stream[i]);
 }

這部分代碼中有一點須要注意：爲了並行化數據拷貝和內核執行，主機端內存必須分配爲鎖頁（page-locked）內存。

要銷燬一個流須要調用函數cudaStreamDestroy()

 for (int i = 0; i < 2; ++i)
     cudaStreamDestroy(stream[i]);

cudaStreamDestroy()函數等待以前流中的指令序列運行完成，而後銷燬指定流，將控制權返還給主機端。

二、默認流（Default stream）

        在內核啓動或者數據拷貝過程當中若是不指定流，或者設置流參數爲0，則相應的指令將會運行在默認流上，它們也所以而順次執行。

三、明同步（Explicit Synchronization）

        在CUDA中有不少種方式能夠用來同步流的執行：

        cudaDeviceSynchronize()函數使得主機端線程阻塞直到全部流中的指令執行完成。

       cudaStreamSynchronize()函數將一個流對象做爲輸入參數，用以等待指定流中的全部指令執行完成。

       cudaStreamWaitEvent()函數將一個流對象和一個事件做爲輸入參數，它將延遲該函數調用後在指定流中全部新加入的命令的執行直到指定的事件完成爲止。流參數能夠爲0，在該情形下全部流中的任何新加入的指令都必須等待指定事件的發生，而後才能夠執行。

       cudaStreamQuery()函數爲應用程序提供了一個檢測指定流中以前指令是否執行完成的方法。

        爲了不同步帶來的性能降低，全部上述同步函數最好用於計時目的或者分離錯誤的內核執行或數據拷貝。

四、暗同步（Implicit Synchronization）

        若是任何一個流中正在執行如下操做，那麼其它流是不能與其並行運行的：

        a. 分配鎖頁內存空間

        b. 設備內存分配

        c. 設備內存置位

        d. 同一設備兩個不一樣地址間正在進行數據拷貝

        e. 默認流中有指令正在執行

        f. L1/shared內存配置的轉換

        對於支持並行內核執行而且計算能力3.0或如下的設備來講，任何一個須要檢查依賴性以肯定流內核啓動是否完成的操做：

        a. 只有當前CUDA上下文中全部流中全部以前的內核啓動以後纔可以啓動執行。

        b. 將會阻塞全部當前CUDA上下文中的任意流中新加入的內核調用直到內核檢查完成。

        須要進行依賴性檢查的操做包括執行檢查的內核啓動所在流中的其它指令以及任何在該流上對cudaStreamQuery()函數的調用。所以，應用程序能夠遵守如下指導原則來提高潛在並行性：

        （1）全部非依賴操做應當比依賴性操做提早進行

        （2）任何類型的同步越遲越好

五、重疊行爲（Overlapping Behavior）

        兩個流間重疊行爲的數量取決於如下幾個因素：

        （1）每一個流中命令發出的次序

        （2）設備是否支持內核啓動與數據傳輸並行

        （3）設備是否支持多內核並行啓動

        （4）設備是否支持多數據傳輸並行

        例如，在不支持並行數據傳輸的設備上，「流的建立與銷燬」章節中代碼樣例中的操做就不能並行，由於在stream[0]中發出設備端到主機端的數據拷貝後，stream[1]又發出主機端到設備端的數據拷貝命令，這兩個命令式不能重疊執行的。假設設備支持數據傳輸與內核啓動並行，那麼以下代碼：

 for (int i = 0; i < 2; ++i)
     cudaMemcpyAsync(inputDevPtr + i * size, hostPtr + i * size, size, cudaMemcpyHostToDevice, stream[i]);
 for (int i = 0; i < 2; ++i)
     MyKernel<<<100, 512, 0, stream[i]>>>(outputDevPtr + i * size, inputDevPtr + i * size, size);
 for (int i = 0; i < 2; ++i)
     cudaMemcpyAsync(hostPtr + i * size, outputDevPtr + i * size, size, cudaMemcpyDeviceToHost, stream[i]);

可將stream[0]的內核啓動和stream[1]從主機端到設備端的數據拷貝重疊起來並行執行。

六、回調函數

        CUDA運行時提供了cudaStreamAddCallback()函數以在流中的任意位置插入一個回調函數點。回調函數運行於主機端，若是在默認流中插入回調函數，那麼它將等待全部其它流中的命令執行完成以後纔會開始執行。

        下面的代碼展現了回調函數技術的應用：

 void CUDART_CB MyCallback(cudaStream_t stream, cudaError_t status, void **data) {
     printf("Inside callback %d\n", (int)data);
 }
 ...
 for (int i = 0; i < 2; ++i) {
     cudaMemcpyAsync(devPtrIn[i], hostPtr[i], size, cudaMemcpyHostToDevice, stream[i]);
     MyKernel<<<100, 512, 0, stream[i]>>>(devPtrOut[i], devPtrIn[i], size);
     cudaMemcpyAsync(hostPtr[i], devPtrOut[i], size, cudaMemcpyDeviceToHost, stream[i]);
     cudaStreamAddCallback(stream[i], MyCallback, (void**)i, 0);
 }

上面的代碼定義了兩個流的操做，每一個流都完成一次主機端到設備端的數據拷貝，一次內核啓動，一次設備端到主機端的數據拷貝，最後增長了一個加入回調函數的操做。當設備端代碼運行到回調函數點的時候，設備將控制權交還給主機端，主機端運行完成之後再將控制權返還給設備端，而後設備端繼續運行。

        值得注意的是，在一個回調函數中，必定不能進行任何CUDA API的調用，直接的或者間接的都是不能夠的。