CUDA C/C++ 教程一：加速應用程式

1. CUDA 基礎

1.1. CUDA 簡介

GPU 加速計算正在逐步取代 CPU 計算，近年來加速計算帶來了越來越多的突破性進展，各類應用程式對加速計算日益增長地需求、便捷地撰寫加速計算的程式的需求以及不斷改進的支持加速計算的硬體設施，所有這一切都在推動著計算方式從 CPU 計算過渡到 GPU 加速計算，

無論是從出色的性能還是易用性來看，CUDA 計算平臺均是加速計算的重要實作方式，CUDA 提供了一種可擴展于 C、C++、Python 和 Fortran 等語言的編碼介面，并行化后的代碼能夠在 NVIDIA GPU 上運行，以大幅加速應用程式，它包含有 DNN、BLAS、圖形分析 和 FFT 等等庫，并且還附帶功能強大的命令列和可視化分析器，

CUDA 支持許多領域的超性能計算應用程式：計算流體動力學、分子動力學、量子化學、物理學 和高性能計算 (HPC)等等，

學習 CUDA 將能幫你加速自己的應用程式，應用程式加速后的執行速度會遠遠超過原本在 CPU 上的執行速度，使那些在 CPU 上性能受限的計算得以進行下去，在本教程中, 你將學習使用 CUDA 的 C/C++ 介面作為加速應用程式編程的入門知識，這些入門知識足以讓你加速自己的 CPU 應用程式，以獲得性能上的巨大提升并幫你邁入全新的計算領域，

1.2. 學習前的準備作業

如要充分利用本教程學習CUDA，那么你應該要先有如下知識儲備：

• 在 C++/C 中宣告變數、撰寫回圈并使用 if/else 陳述句，
• 在 C++/C 中定義和呼叫函式，
• 在 C++/C 中分配陣列，

說白了就是要有C或C++語言的基礎，此外不需要事先知道任何關于 CUDA 的知識，當你在本教程完成學習后，你就可以做到：

• 撰寫、編譯及運行既可呼叫 CPU 函式也可啟動 GPU 核函式的 C/C++ 程式，
• 通過配置引數控制并行執行緒的層次結構，
• 重構串行回圈以在 GPU 上并行執行其迭代，
• 分配和釋放可用于 CPU 和 GPU 的記憶體，
• 處理 CUDA 代碼產生的錯誤，
• 加速 CPU 應用程式，

1.3. 加速系統的硬體設施

帶有GPU的計算機系統稱為加速系統（又稱異構系統，即指包含CPU和GPU的系統），在一個包含 NVIDIA GPU 的加速系統的實驗環境上，可以使用 nvidia-smi 命令查詢有關此 GPU 的資訊，例如：

nvidia-smi

按回車之后，將輸出該機器上的GPU資訊

需要注意的是，加速系統在運行程式時首先會運行 CPU 程式，在運行到需要GPU進行大規模并行計算的函式時，再將對應函式載入GPU執行，

也就是說，由GPU加速的依然還是純CPU的應用程式，只是某些模塊在運行時調入了GPU中，該模塊在同步完畢后將會重新回到CPU中執行主程式的后續代碼：

2. 撰寫在GPU運行的代碼

CUDA 為許多編程語言提供了擴展介面，而在本教程用CUDA為 C/C++ 提供的介面來展示，對編程語言的擴展可以讓開發人員在 GPU 上更加方便的運行 CUDA 庫的函式，

以下是一個 .cu 檔案（.cu 是 CUDA 加速程式的檔案擴展名，實際上.cu檔案只是含有CUDA代碼的.cpp檔案，沒有別的特殊之處），其中包含兩個函式，第一個函式 CPUFunction() 將在 CPU 上運行，第二個函式 GPUFunction() 將在 GPU 上運行：

// 在CPU上運行的函式
void CPUFunction() {
  printf("This function is defined to run on the CPU.\n");
}

// 在GPU上運行的函式
__global__ void GPUFunction() {
  printf("This function is defined to run on the GPU.\n");
}

int main() {
  CPUFunction(); // 呼叫CPU函式

  GPUFunction<<<1, 1>>>(); // 呼叫GPU函式
  cudaDeviceSynchronize(); // 同步
}

根據上面的代碼，我們來講解一些需要特別注意的重要代碼行，以及加速計算中使用的一些其他常用術語：
__global__ void GPUFunction()：

• __global__ 關鍵字表明該函式將在 GPU 上運行并可全域呼叫（ 既可以由CPU ，也可以由 GPU 呼叫）；
• 通常，我們將在 CPU 上執行的代碼稱為 Host （主機）代碼，而將在 GPU 上運行的代碼稱為 Device （設備）代碼；
• 注意回傳型別為 void，使用 __global__ 關鍵字定義的函式回傳值需為 void 型別，

GPUFunction<<<1, 1>>>()：

• 通常，我們把要運行在 GPU 上的函式稱為 kernel （核）函式;
• 啟動核(kernel)函式時，我們必須事先配置GPU引數，使用 <<< ... >>> 語法向核函式傳遞兩個必要的引數;
• 在 <<< ... >>> 中傳遞的引數用于為核函式設定執行緒的層次結構，第一個引數定義執行緒塊(Block)的數量，第二個引數定義Block中含有的執行緒(Thread)數量，例如本例中的核函式 GPUFunction() 將在包含 1 個執行緒（第二個配置引數）的 1 個執行緒塊（第一個執行配置引數）上運行，

cudaDeviceSynchronize()：

• 與其他并行化的代碼類似，核函式啟動方式為異步，即 CPU 代碼將繼續執行而不會等待核函式執行完成；
• 呼叫 CUDA 提供的函式 cudaDeviceSynchronize 可以讓Host 代碼(CPU) 等待 Device 代碼(GPU) 執行完畢，再在CPU上繼續執行，

2.1. 撰寫運行一個 Hello GPU 核函式

#include <stdio.h>

void helloCPU() {
  printf("Hello from the CPU.\n");
}

// __global__ 表明這是一個全域GPU核函式.
__global__ void helloGPU() {
  printf("Hello from the GPU.\n");
}

int main() {
  helloCPU(); // 呼叫CPU函式

   /* 使用 <<<...>>> 配置核函式的GPU引數，
   * 第一個1表示1個執行緒塊，第二個1表示每個執行緒塊1個執行緒，*/
  helloGPU<<<1, 1>>>(); // 呼叫GPU函式
  cudaDeviceSynchronize(); // `cudaDeviceSynchronize` 同步CPU和GPU
}

現在來編譯并運行加速后的CUDA代碼，將上述檔案命名為hello-gpu.cu，執行命令：

nvcc hello-gpu.cu -o hello-gpu
./hello-gpu

得到結果：

3. CUDA 執行緒的層次結構

從上面的圖中可以看出，CUDA執行緒的層次結構分為三層：Thread（執行緒）、Block（塊）、Grid（網格），網格由塊組成，塊由執行緒組成，

3.1. 運行核函式

我們可以通過配置引數指定核函式如何在 GPU 的多個執行緒中并行運行，具體來說，就可以配置 Block 的數量以及每個 Block 中所包含 Thread 的數量，配置引數的語法如下：

<<< Block 數, 每個Block中的 Thread 數>>>

啟動核函式時，核函式代碼由我們自行配置的 Block 中的每個 Thread 執行，因此，如果假設已定義一個名為 someKernel 的核函式，則GPU執行緒可以配置為下列情況：

• someKernel<<<1, 1>>() 在GPU中為該核函式分配1個具有1個執行緒的執行緒塊，核函式中的代碼將只運行1次；
• someKernel<<<1, 10>>() 在GPU中為該核函式分配1個具有10個執行緒的執行緒塊，核函式中的代碼將運行10次；
• someKernel<<<10, 1>>() 在GPU中為該核函式分配10個具有1個執行緒的執行緒塊，核函式中的代碼將運行10次；
• someKernel<<<10, 10>>() 在GPU中為該核函式分配10個具有10個執行緒的執行緒塊，核函式中的代碼將運行100次；

啟動并行運行的核函式示例：

#include <stdio.h>

__global__ void firstParallel() {
  printf("This is running in parallel.\n");
}

int main() {
  firstParallel<<<5, 5>>>(); // 在GPU中為核函式分配5個具有5個執行緒的執行緒塊，將運行25次；
  cudaDeviceSynchronize(); // 同步
}

將上述代碼命名為basic-parallel.cu，然后編譯運行：

nvcc basic-parallel.cu -o basic-parallel
./basic-parallel

結果如下，數了一下，確實是25次：

3.2. 執行緒和塊的索引

如圖所示，每個執行緒在其執行緒塊的內部都會被分配一個索引，從 0 開始，此外，每個執行緒塊也會被分配一個索引，也是從 0 開始，正如執行緒組成執行緒塊，執行緒塊又會組成網格(Grid)，而網格是 CUDA 執行緒層次結構中級別最高的物體，它沒有索引，

簡言之，CUDA 核函式在由一個或多個執行緒塊組成的網格中執行，且每個執行緒塊中均包含相同數量的一個或多個執行緒（每個執行緒塊中的執行緒數量相同），

在核函式中，可以通過兩個變數來獲取到索引： threadIdx.x （執行緒索引）和 blockIdx.x（執行緒塊索引），

現在讓我們來使用索引控制特定的執行緒和塊：

#include <stdio.h>

// 核函式
__global__ void printSuccessForCorrectExecutionConfiguration() {
  // 當執行到第255個執行緒塊的第1023個執行緒時，才輸出
  if(threadIdx.x == 1023 && blockIdx.x == 255) {
    printf("Success!\n"); // 輸出 Success！
    printf("threadIdx.x: %d\n", threadIdx.x); // 輸出執行緒ID
    printf("blockIdx.x: %d\n", blockIdx.x); // 輸出執行緒塊ID
  }
}

int main() {
  // 配置該核函式由256個含有1024個執行緒的執行緒塊中執行
  printSuccessForCorrectExecutionConfiguration<<<256, 1024>>>();
  cudaDeviceSynchronize(); // 同步
}

將上述代碼命名為thread-and-block-idx.cu，然后編譯運行：

nvcc thread-and-block-idx.cu -o thread-and-block-idx
./thread-and-block-idx

輸出：

3.3. 用 CUDA 加速 For 回圈

到此為止，加速 for 回圈就是一個可行的操作了，在加速計算中，for 回圈不再順序執行每次迭代，而是讓每次迭代都在不同的執行緒中并行執行，

例如，現在有以下在 CPU 中執行的 for 回圈：

int N = 10;
for (int i = 0; i < N; ++i) {
  printf("%d\n", i);
}

如要并行此回圈，必須執行以下 2 個步驟：

• 撰寫用于執行單次迭代作業的核函式，
• 呼叫核函式時為它配置執行引數，即并行的執行緒數，每個執行緒執行一次迭代，

如下例程式：

#include <stdio.h>

// 核函式
__global__ void loop() {
  // 輸出每一個執行緒的執行緒號（0~9）
  printf("This is iteration number %d\n", threadIdx.x);
}

int main() {
  loop<<<1, 10>>>(); // 執行核函式
  cudaDeviceSynchronize();
}

將上述代碼命名為single-block-loop.cu，然后編譯運行：

nvcc single-block-loop.cu -o single-block-loop
./single-block-loop

輸出：

3.4. 管理不同塊之間的執行緒

之前提到過，一個執行緒塊可以包含多個執行緒，那么我們就可以調整執行緒塊的大小以實作更多型別的并行化，執行緒塊包含的執行緒具有數量限制：確切地說是 1024 個（即每個塊中的執行緒數量 <= 1024），通常為了增加加速應用程式中的并行量，我們需要利用多個執行緒塊，并在它們之間進行協調，

CUDA 核函式中，記錄了每個塊中執行緒數的變數是 blockDim.x（一個執行緒塊中包含的執行緒數量，每個塊中包含的執行緒數都是一樣的），通過將此變數與 blockIdx.x 和 threadIdx.x 變數結合使用，并借助運算式 threadIdx.x + blockIdx.x * blockDim.x 計算執行緒ID，該運算式可以用C++中訪問二維陣列的索引計算來類比看待，以增強理解，

以下是詳細示例：

配置引數 <<<10, 10>>> 將啟動共計擁有 100 個執行緒的網格，該網格又分為由 10 個執行緒組成的 10 個執行緒塊（即一個執行緒塊中含有10個執行緒，blockDim.x=10），這時候，就可以利用運算式 threadIdx.x + blockIdx.x * blockDim.x 來計算某個執行緒的唯一索引（0 至 99 之間）了，

• 如果執行緒塊 blockIdx.x 索引為 0，則 blockIdx.x * blockDim.x 為 0，以 0 為起始索引加上可能的 threadIdx.x 值（0 至 9），便可在網格中找到索引為 0 至 9 的執行緒，
• 如果執行緒塊 blockIdx.x 索引為 1，則 blockIdx.x * blockDim.x 為 10，以 10 為起始索引加上可能的 threadIdx.x 值（0 至 9），便可在網格中找到索引為 10 至 19 的執行緒，
• 如果執行緒塊 blockIdx.x 索引為 5，則 blockIdx.x * blockDim.x 為 50，以 50 為起始索引加上可能的 threadIdx.x 值（0 至 9），便可在網格中找到索引為 50 至 59 的執行緒，
• 如果執行緒塊 blockIdx.x 索引為 9，則 blockIdx.x * blockDim.x 為 90，以 90 為起始索引加上可能的 threadIdx.x 值（0 至 9），便可在網格中找到索引為 90 至 99 的執行緒，

現在我們來加速具有多個執行緒塊的For回圈：

#include <stdio.h>

__global__ void loop()
{
  // 在Grid中遍歷所有thread
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  printf("%d\n", i);
}

int main()
{
  /*
   * 配置引數還可以試試其他的，例如：
   * <<<5, 2>>>
   * <<<10, 1>>>
   */
  loop<<<2, 5>>>();
  cudaDeviceSynchronize();
}

將上述代碼命名為multi-block-loop.cu，然后編譯運行：

nvcc multi-block-loop.cu -o multi-block-loop
./multi-block-loop

輸出：

4. 分配可同時被GPU和CPU訪問的記憶體

CUDA 的最新版本（版本 6 和更高版本）可以便捷地分配和釋放既可用于 Host 也可被 Device 訪問的記憶體，

在 Host（CPU）中，我們一般適用malloc 和 free 來分配和釋放記憶體，但這樣分配的記憶體無法直接被Device（GPU）訪問，所以在這里我們用cudaMallocManaged 和 cudaFree 兩個函式來分配和釋放同時可被 Host 和 Device 訪問的記憶體，如下例所示：

// CPU
int N = 10;
size_t size = N * sizeof(int);
int *a;

a = (int *)malloc(size); // 分配CPU記憶體
free(a); // 釋放CPU記憶體

// GPU
int N = 10;
size_t size = N * sizeof(int);
int *a;

cudaMallocManaged(&a, size);// 為a分配CPU和GPU記憶體
cudaFree(a); // 釋放GPU記憶體

實際上，cudaMallocManaged在統一記憶體中創建了一個托管記憶體池（CPU上有，GPU上也有），記憶體池中已分配的空間可以通過相同的指標直接被CPU和GPU訪問，底層系統在統一的記憶體空間中自動地在設備和主機間進行傳輸，資料傳輸對應用來說是透明的，大大簡化了代碼，

現在讓我們來看看如何利用GPU來執行陣列元素的乘法操作：

#include <stdio.h>

// 初始化陣列
void init(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    a[i] = i;
  }
}

// CUDA 核函式，所有元素乘2
__global__ void doubleElements(int *a, int N) {
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < N) {
    a[i] *= 2;
  }
}

// 檢查陣列內所有元素的值是否均為復數
bool checkElementsAreDoubled(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    if (a[i] != i*2) return false;
  }
  return true;
}

int main() {
  int N = 1000;
  int *a;

  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size); // 為a分配CPU和GPU空間

  init(a, N); // 為陣列a賦值
  size_t threads_per_block = 256; // 定義每個block的thread數量
  size_t number_of_blocks = (N + threads_per_block - 1) / threads_per_block; // 定義block的數量

  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N); // 執行核函式
  cudaDeviceSynchronize(); // 同步

  bool areDoubled = checkElementsAreDoubled(a, N); // 檢查元素是否為復數
  printf("All elements were doubled? %s\n", areDoubled ? "TRUE" : "FALSE");

  cudaFree(a); // 釋放由cudaMallocManaged
}

將上述代碼命名為double-elements.cu，然后編譯運行：

nvcc double-elements.cu -o double-elements
./double-elements

輸出：

5. 網格大小與實際并行作業量不匹配

5.1. 網格大于作業量

鑒于 GPU 的硬體特性，執行緒塊中的執行緒數最好配置為 32 的倍數，但是在實際作業中，很可能會出現這樣的情況，我們手動配置引數所創建的執行緒數無法匹配為實作并行回圈所需的執行緒數，比如實際上需要執行1230次回圈，但是你卻配置了2048個執行緒，

我們不可能每次配置引數的時候都手動去算一遍最佳配置，更何況并不是所有的數都是 32 的倍數，不過這個問題現在已經可以通過以下三個步驟輕松地解決：

• 首先，設定配置引數，使執行緒總數超過實際作業所需的執行緒數，
• 然后，在向核函式傳遞引數時傳遞一個用于表示要處理的資料集總大小或完成作業所需的總執行緒數 N，
• 最后，計算網格內的執行緒索引后（使用 threadIdx + blockIdx*blockDim），判斷該索引是否超過 N，只在不超過的情況下執行與核函式相關的作業，

以下是一種可選的配置方式，適用于 作業總量 N 和執行緒塊中的執行緒數已知的情況，如此一來，便可確保網格中至少始終能執行 N 次任務，且最多只浪費 1 個執行緒塊的執行緒數量：

// 假設N是已知的
int N = 100000;

// 把每個block中的thread數設為256
size_t threads_per_block = 256;

// 根據N和thread數量配置Block數量
size_t number_of_blocks = (N + threads_per_block - 1) / threads_per_block;

some_kernel<<<number_of_blocks, threads_per_block>>>(N);

由于上述執行配置致使網格中的執行緒數超過 N，因此需要注意 some_kernel 定義中的內容，以確保 some_kernel 在由其中一個額外的（大于N的）執行緒執行時不會嘗試訪問超出范圍的資料元素，也就是要加個判斷：

__global__ some_kernel(int N) {
  int idx = threadIdx.x + blockIdx.x * blockDim.x;

  if (idx < N) { // 保證執行緒ID小于元素數量N
    // 并行代碼
  }

使用不匹配的配置引數來加速 For 回圈

#include <stdio.h>

__global__ void initializeElementsTo(int initialValue, int *a, int N) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  if (i < N) {
    a[i] = initialValue;
  }
}

int main() {
  int N = 1000;

  int *a;
  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size);

  size_t threads_per_block = 256;
  // 這是慣用的CUDA語法
  // 為 number_of_blocks 分配一個值，以確保執行緒數至少與指標 a 中可供訪問的元素數同樣多，
  size_t number_of_blocks = (N + threads_per_block - 1) / threads_per_block;

  int initialValue = 6; // 初始化的值
  initializeElementsTo<<<number_of_blocks, threads_per_block>>>(initialValue, a, N);
  cudaDeviceSynchronize();

  // 檢查元素值是否被初始化
  for (int i = 0; i < N; ++i) {
    if(a[i] != initialValue) {
      printf("FAILURE: target value: %d\t a[%d]: %d\n", initialValue, i, a[i]);
      exit(1);
    }
  }
  printf("SUCCESS!\n");

  cudaFree(a);
}

將上述代碼命名為mismatched-config-loop.cu，然后編譯運行：

nvcc mismatched-config-loop.cu -o mismatched-config-loop
./mismatched-config-loop

輸出：

5.2. 網格小于作業量

有時，作業量比網格大，或者出于某種原因，一個網格中的執行緒數量可能會小于實際作業量的大小，請思考一下包含 1000 個元素的陣列和包含 250 個執行緒的網格（此處使用極小的規模以便于說明），此網格中的每個執行緒將需使用 4 次，如要實作此操作，一種常用方法便是在核函式中使用跨網格回圈，

在跨網格回圈中，每個執行緒將在網格內使用 threadIdx + blockIdx*blockDim 計算自身唯一的索引，并對陣列內該索引的元素執行相應運算，然后用網格中的執行緒數加上自身索引值，并重復此操作，直至超出陣列范圍，

例如，對于包含 500 個元素的陣列 a 和包含 250 個執行緒的網格，網格中索引為 20 的執行緒將執行如下操作：

• 對 a[20] 執行相應運算；
• 將執行緒索引增加 250，使網格的大小達到 270
• 對a[270] 執行相應運算；
• 將執行緒索引增加 250，使網格的大小達到 520
• 由于 520 現已超出陣列范圍，因此執行緒將停止作業，

CUDA 提供一個記錄了網格中執行緒塊數的變數：gridDim.x，然后可以利用它來計算網格中的總執行緒數，即網格中的執行緒塊數乘以每個執行緒塊中的執行緒數：gridDim.x * blockDim.x，現在來看看以下核函式中網格跨度回圈的示例：

__global void kernel(int *a, int N)
{
  int indexWithinTheGrid = threadIdx.x + blockIdx.x * blockDim.x;
  int gridStride = gridDim.x * blockDim.x; // grid 的一個跨步

  for (int i = indexWithinTheGrid; i < N; i += gridStride) {
    // 對 a[i] 的操作;
  }
}

上面是一個簡單的例子，現在我們來看看一個更詳細的例子，使用了跨網格回圈來處理比網格更大的陣列：

#include <stdio.h>

// 初始化陣列a
void init(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    a[i] = i;
  }
}

__global__ void doubleElements(int *a, int N) {

  // 使用grid-stride回圈，這樣每個執行緒可以處理陣列中的多個元素，
  int idx = blockIdx.x * blockDim.x + threadIdx.x;
  int stride = gridDim.x * blockDim.x; // grid 的一個跨步

  for (int i = idx; i < N; i += stride) {
    a[i] *= 2;
  }
}

// 檢查陣列內所有元素的值是否均為復數
bool checkElementsAreDoubled(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    if (a[i] != i*2) return false;
  }
  return true;
}

int main() {
  int N = 10000;
  int *a;
  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size);

  init(a, N); // 初始化陣列a

  size_t threads_per_block = 256; // 每個block的thread數量
  size_t number_of_blocks = 32; // block數量

  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N);
  cudaDeviceSynchronize();

  bool areDoubled = checkElementsAreDoubled(a, N); 
// 檢查陣列內所有元素的值是否均為復數
  printf("All elements were doubled? %s\n", areDoubled ? "TRUE" : "FALSE");

  cudaFree(a);
}

將上述代碼命名為grid-stride-double.cu，然后編譯運行：

nvcc grid-stride-double.cu -o grid-stride-double
./grid-stride-double

輸出：

6. 錯誤處理

CUDA 函式發生錯誤時會回傳一個型別為 cudaError_t 的變數，該變數可用于檢查呼叫函式時是否發生錯誤，以下是對呼叫 cudaMallocManaged 函式執行錯誤處理的示例：

cudaError_t err;
err = cudaMallocManaged(&a, N)                    // 假設a和N已經被定義

if (err != cudaSuccess) { // `cudaSuccess` 是一個 CUDA 變數.
  printf("Error: %s\n", cudaGetErrorString(err)); // `cudaGetErrorString` 是一個 CUDA 函式.
}

但是，核函式并不會回傳型別為 cudaError_t 的值（因為核函式的回傳值為void），為檢查執行核函式時是否發生錯誤（例如配置錯誤），CUDA 提供了 cudaGetLastError 函式，可以用于檢查核函式執行期間發生的錯誤，

// 這段程式中的核函式會出一個CUDA錯誤，但是核函式本身無法捕獲該錯誤
someKernel<<<1, -1>>>();  // 執行緒數不能為-1

cudaError_t err;
err = cudaGetLastError(); // `cudaGetLastError` 會捕獲上面代碼中的最近的一個錯誤
if (err != cudaSuccess) {
  printf("Error: %s\n", cudaGetErrorString(err));
}

另一個要注意的點是，為了捕捉在異步核函式執行期間發生的錯誤，一定要檢查后續同步 CPU 與 GPU 時 API 呼叫所回傳的狀態（例如 cudaDeviceSynchronize）；如果之前執行的某一個核函式失敗了，則將會回傳錯誤，

添加錯誤處理的示例：

#include <stdio.h>

// 初始化陣列a
void init(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    a[i] = i;
  }
}

// CUDA 核函式 陣列元素值乘2
__global__ void doubleElements(int *a, int N) {
  int idx = blockIdx.x * blockDim.x + threadIdx.x;
  int stride = gridDim.x * blockDim.x;

  // for (int i = idx; i < N; i += stride) {
  // 這里出現一個數值越界錯誤
  for (int i = idx; i < N + stride; i += stride) {
    a[i] *= 2;
  }
}

// 檢查陣列元素是否均為復數
bool checkElementsAreDoubled(int *a, int N) {
  int i;
  for (i = 0; i < N; ++i) {
    if (a[i] != i*2) return false;
  }
  return true;
}

int main() {
  int N = 10000;
  int *a;

  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size);
  init(a, N);

  cudaError_t syncErr, asyncErr; // 定義錯誤處理變數
  
  // size_t threads_per_block = 1024;
  // 執行緒數大于1024（前面說過每個block的執行緒數不能超過1024）
  size_t threads_per_block = 2048;
  size_t number_of_blocks = 32;
  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N); // 執行核函式

  syncErr = cudaGetLastError(); // 捕獲核函式執行期間發生的錯誤
  asyncErr = cudaDeviceSynchronize(); // 同步，并捕獲同步期間發生的錯誤

  // 輸出錯誤 說明：兩個錯誤需分別設定（即每次運行時只保留一個錯誤）
  if (syncErr != cudaSuccess) printf("Error: %s\n", cudaGetErrorString(syncErr));
  if (asyncErr != cudaSuccess) printf("Error: %s\n", cudaGetErrorString(asyncErr));

  bool areDoubled = checkElementsAreDoubled(a, N); // 驗證陣列元素值是否均為復數
  printf("All elements were doubled? %s\n", areDoubled ? "TRUE" : "FALSE");

  cudaFree(a);
}

將上述代碼命名為add-error-handling.cu，然后編譯運行：

nvcc add-error-handling.cu -o add-error-handling
./add-error-handling

輸出：

6.1. 定制一個 CUDA 錯誤處理宏

創建一個包裝 CUDA 函式呼叫的宏對于檢查錯誤十分有用，以下是一個宏示例，我們可以在其他的 CUDA 代碼中隨時使用：

#include <stdio.h>
#include <assert.h>

// CUDA 錯誤處理宏
inline cudaError_t checkCuda(cudaError_t result)
{
  if (result != cudaSuccess) {
    fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
    assert(result == cudaSuccess);
  }
  return result;
}

int main() {
  // checkCuda 宏可以回傳 CUDA 函式回傳的錯誤型別`cudaError_t`的值
  checkCuda( cudaDeviceSynchronize() )
}

7. 總結

至此，我們已經完成了我們預期的學習目標：

• 撰寫、編譯及運行既可呼叫 CPU 函式也可啟動GPU核函式的 C/C++ 程式，
• 使用執行配置控制并行執行緒層次結構，
• 重構串行回圈以在 GPU 上并行執行其迭代，
• 分配和釋放可用于 CPU 和 GPU 的記憶體，
• 處理 CUDA 代碼生成的錯誤，

現在，加速 CPU 應用程式進行是可行的了，

7.1 用 CUDA 實作向量加法

為了展示一下如何綜合運用本篇教程提到的內容，我們通過一個向量與向量加分的案例來串用以上知識：

#include <stdio.h>
#include <assert.h>

// CUDA 錯誤處理宏
inline cudaError_t checkCuda(cudaError_t result)
{
  if (result != cudaSuccess) {
    fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
    assert(result == cudaSuccess);
  }
  return result;
}

// 初始化陣列 a
void initWith(float num, float *a, int N) {
  for(int i = 0; i < N; ++i) {
    a[i] = num;
  }
}

// 向量加法核函式
__global__ void addVectorsInto(float *result, float *a, float *b, int N) {
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  int stride = blockDim.x * gridDim.x;

  for(int i = index; i < N; i += stride) {
    result[i] = a[i] + b[i]; // 元素a[i] + 元素 b[i]
  }
}

// 檢查 CUDA 向量加分是否計算正確
void checkElementsAre(float target, float *array, int N) {
  for(int i = 0; i < N; i++) {
    if(array[i] != target) {
      printf("FAIL: array[%d] - %0.0f does not equal %0.0f\n", i, array[i], target);
      exit(1);
    }
  }
  printf("SUCCESS! All values added correctly.\n");
}

int main() {
  const int N = 10;
  size_t size = N * sizeof(float);

  float *a;
  float *b;
  float *c;

  // 分配記憶體，且檢查執行期間發生的錯誤
  checkCuda( cudaMallocManaged(&a, size) );
  checkCuda( cudaMallocManaged(&b, size) );
  checkCuda( cudaMallocManaged(&c, size) );

  initWith(3, a, N); // 將陣列a中所有的元素初始化為3
  initWith(4, b, N); // 將陣列b中所有的元素初始化為4
  initWith(0, c, N); // 將陣列c中所有的元素初始化為0，陣列c是結果向量

  // 配置引數
  size_t threadsPerBlock = 256;
  size_t numberOfBlocks = (N + threadsPerBlock - 1) / threadsPerBlock;
  addVectorsInto<<<numberOfBlocks, threadsPerBlock>>>(c, a, b, N); // 執行核函式

  checkCuda( cudaGetLastError() ); // 檢查核函式執行期間發生的錯誤
  checkCuda( cudaDeviceSynchronize() ); // 同步，且檢查執行期間發生的錯誤

  checkElementsAre(7, c, N);  // 檢查向量加的結果是否正確

  // 釋放記憶體，且檢查執行期間發生的錯誤
  checkCuda( cudaFree(a) );
  checkCuda( cudaFree(b) );
  checkCuda( cudaFree(c) );
}

7.2. 二維和三維的網格和塊

網格和執行緒塊最多可以定義有 3 個維度，使用多個維度定義網格和執行緒塊在處理具有多個維度的資料時可能很有效，例如二維矩陣，如果要定義二維或三維的網格或執行緒塊，可以使用 CUDA 的 dim3 關鍵字來定義多維網格或塊，即如下所示：

dim3 threads_per_block(16, 16, 1);
dim3 number_of_blocks(16, 16, 1);
someKernel<<<number_of_blocks, threads_per_block>>>();

鑒于以上示例，someKernel 內部的變數 gridDim.x、gridDim.y、blockDim.x 和 blockDim.y 均將等于 16，

7.3 用 CUDA 實作矩陣乘法

#include <stdio.h>
#define N  64

// GPU 矩陣乘法
__global__ void matrixMulGPU( int * a, int * b, int * c ) {
  int val = 0;

  int row = blockIdx.x * blockDim.x + threadIdx.x;
  int col = blockIdx.y * blockDim.y + threadIdx.y;

  if (row < N && col < N) {
    for ( int k = 0; k < N; ++k )
      val += a[row * N + k] * b[k * N + col];
    c[row * N + col] = val;
  }
}

// CPU矩陣乘法
void matrixMulCPU( int * a, int * b, int * c ) {
  int val = 0;

  for( int row = 0; row < N; ++row )
    for( int col = 0; col < N; ++col ) {
      val = 0;
      for ( int k = 0; k < N; ++k )
        val += a[row * N + k] * b[k * N + col];
      c[row * N + col] = val;
    }
}

int main() {
  int *a, *b, *c_cpu, *c_gpu;
  int size = N * N * sizeof (int); // Number of bytes of an N x N matrix

  // 分配記憶體
  cudaMallocManaged (&a, size);
  cudaMallocManaged (&b, size);
  cudaMallocManaged (&c_cpu, size);
  cudaMallocManaged (&c_gpu, size);

  // 初始化陣列
  for( int row = 0; row < N; ++row )
    for( int col = 0; col < N; ++col )
    {
      a[row * N + col] = row;
      b[row * N + col] = col + 2;
      c_cpu[row * N + col] = 0;
      c_gpu[row * N + col] = 0;
    }

  dim3 threads_per_block (16, 16, 1); // 一個 16 * 16 的執行緒陣
  dim3 number_of_blocks ((N / threads_per_block.x) + 1, (N / threads_per_block.y) + 1, 1);
  matrixMulGPU <<< number_of_blocks, threads_per_block >>> ( a, b, c_gpu ); // 執行核函式

  cudaDeviceSynchronize(); // 同步

  matrixMulCPU( a, b, c_cpu ); // 執行 CPU 版本的矩陣乘法

  // 比較 CPU 和 GPU 兩種方法的計算結果是否一致
  bool error = false;
  for( int row = 0; row < N && !error; ++row )
    for( int col = 0; col < N && !error; ++col )
      if (c_cpu[row * N + col] != c_gpu[row * N + col]) {
        printf("FOUND ERROR at c[%d][%d]\n", row, col);
        error = true;
        break;
      }
  if (!error)
    printf("Success!\n");

  // 釋放記憶體
  cudaFree(a); cudaFree(b);
  cudaFree( c_cpu ); cudaFree( c_gpu );
}