【CUDA並行程序設計系列（3）】CUDA線程模型

前言

CUDA並行程序設計系列是本人在學習CUDA時整理的資料，內容大都來源於對《CUDA並行程序設計：GPU編程指南》、《GPU高性能編程CUDA實戰》和CUDA Toolkit Documentation的整理。通過本系列整體介紹CUDA並行程序設計。內容包括GPU簡介、CUDA簡介、環境搭建、線程模型、內存、原子操作、同步、流和多GPU架構等。

本系列目錄：

• 【CUDA並行程序設計系列（1）】GPU技術簡介
• 【CUDA並行程序設計系列（2）】CUDA簡介及CUDA初步編程
• 【CUDA並行程序設計系列（3）】CUDA線程模型
• 【CUDA並行程序設計系列（4）】CUDA內存
• 【CUDA並行程序設計系列（5）】CUDA原子操作與同步
• 【CUDA並行程序設計系列（6）】CUDA流與多GPU
• 關於CUDA的一些學習資料

在前一章的代碼雖然是在GPU上執行，但並沒有涉及到並行運行。GPU計算的應用情景在很大程度上取決於能否從許多問題中發掘出大規模並行性。本文將介紹CUDA的線程模型，通過實例代碼（矢量加法和矩陣加法）演示CUDA下的並行程序設計。

矢量加法

我們通過比較在CPU和GPU的矢量加法的例子來看CUDA的並行執行，矢量加法CPU的代碼如下：

#include <stdio.h>
#define N 10
void vecAdd(int *a, int *b, int *c)
{
      int tid = 0;
      while (tid < N)
      {
            c[tid] = a[tid] + b[tid];
            ++tid;
      }
}
int main(void)
{
      int a[N], b[N], c[N];

      for (int i = 0; i < N; ++i)
      {
            a[i] = -i;
            b[i] = i * i;
      }
      vecAdd (a, b, c);
      for (int i = 0; i < N; ++i)
      {
            printf("%d + %d = %d\n", a[i], b[i], c[i]);
      }
      return 0;
}

GPU實現的代碼如下：

#include <stdio.h>
#define N 10
__global__ void vecAdd(int *a, int *b, int *c)
{
      int tid = threadIdx.x;
      if (tid < N)
      {
            c[tid] = a[tid] + b[tid];
      }
}
int main(void)
{
      int a[N], b[N], c[N];
      int *dev_a, *dev_b, *dev_c;
      for (int i = 0; i < N; ++i)
      {
            a[i] = -i;
            b[i] = i * i;
      }
      cudaMalloc( (void **)&dev_a, N * sizeof(int) );
      cudaMalloc( (void **)&dev_b, N * sizeof(int) );
      cudaMalloc( (void **)&dev_c, N * sizeof(int) );
      cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice );
      cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice );
      vecAdd<<<1, N>>>(dev_a, dev_b, dev_c);
      cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost );
      for (int i = 0; i < N; ++i)
      {
            printf("%d + %d = %d\n", a[i], b[i], c[i]);
      }
      cudaFree(dev_a);
      cudaFree(dev_b);
      cudaFree(dev_c);
      return 0;
}

關於cudaMalloc()、cudaMemcpy()在前一章已經講過，這里不再解釋。通過比較可以看到，這兩段代碼非常相似，主要區別有兩點：

• 調用核函數時使用vecAdd<<<1, N>>>(dev_a, dev_b, dev_c)。

• 核函數vecAdd中的for循環不見了，還使用了一個未定義的threadIdx。

使用<<<...>>>啟動核函數是告訴運行時如何啟動核函數，在這兩個參數中，第一個參數表示設備在執行核函數數使用的並行線程塊（block）的數量，第二個參數表示每個線程塊啟動線程（thread）的數量。顯然，本例中啟動了1個線程塊，每個線程塊啟動N個線程，每個線程將並行運行核函數。關於線程、線程塊將在下一節解釋。

這種運行方式引出了一個問題：既然GPU將運行核函數的N各副本，那如何在代碼中知道當前正在運行的是哪一個線程？這時就要使用內置的threadIdx，threadIdx.x代表當前線程的編號，這樣就可以計算位於這個索引出得數據，N個線程，每個線程處理一個索引，就完成了矢量加法。

代碼也可以改下成如下方式：

__global__ void vecAdd(int *a, int *b, int *c)
{
      int tid = blockIdx.x;
      if (tid < N)
      {
            c[tid] = a[tid] + b[tid];
      }
}

int main(void)
{
···
      vecAdd<<<N, 1>>>(dev_a, dev_b, dev_c);
···
}

這將啟動N個線程塊，每個線程塊啟動一個線程，內置的blockIdx.x代表當前線程塊的編號。

線程、線程塊以及線程網格

線程是並行程序的基本構建塊，和C語言的線程是一樣的。在前面已經看到，一個線程塊中可以有多個線程。事實上，線程塊又組成了一個線程網絡（grid）。線程網格（Grid）、線程塊（Block）和線程（Thread）的結構圖如下所示：

<<<...>>> 中的參數類型是int或dim3，dim3是CUDA定義的一種數據類型，代表一個三維的數據(x,y,z)，不過目前為止第三維只能是1。因此，可以用這個數據結構創建一個二維的線程塊和線程網絡。如上圖，每個線程塊的X軸方向上開啟了4個線程，Y軸方向上開啟了3個線程，每個線程塊開啟了4 * 3=12個線程。在線程網絡上，X軸方向上開啟了3個線程塊，Y軸方向上開啟了2個線程塊，線程網格共有3 * 2=6個線程塊。整個線程網絡線程總數為(4 * 3)*(3 * 2)=72，當然，這只是一個小的例子，在實際的程序的線程數可能是這個數量的1萬倍甚至更多。

由於程序中可能不止用到一個維度的線程索引，有可能用到X軸和Y軸兩個維度，因此，前面在核函數中使用threadIdx.x、blockIdx.x就不適用了，需要修改核函數，計算不同維度的索引。從上圖也可以看出，threadIdx.x、blockIdx.x成了相對索引。除了計算不同維度的相對索引外，有時可能還需要線性計算出相對整個線程網絡的絕對線程索引。因此，CUDA定義了一些新的數據結構來保存線程索引相關量：

• gridDim.x-線程網絡X維度上線程塊的數量
• gridDim.y-線程網絡Y維度上線程塊的數量
• blockDim.x-一個線程塊X維度上的線程數量
• blockDim.y-一個線程塊Y維度上的線程數量
• blockIdx.x-線程網絡X維度上的線程塊索引
• blockIdx.y-線程網絡Y維度上的線程塊索引
• threadIdx.x-線程塊X維度上的線程索引
• threadIdx.y-線程塊Y維度上的線程索引

圖中的各個對應的值為：

• gridDim.x=3
• gridDim.y=2
• blockDim.x=4
• blockDim.y=3

要計算當前線程相對於整個線程網絡的絕對索引，可使用如下代碼：

idx = (blockIdx.x * blockDim.x) + threadIdx.x;
idy = (blockIdx.y * blockDim.y) + threadIdx.y;
thread_idx = ((gridDim.x * blockDim.y) * idy) + idx;

矩陣加法

本節用一個計算矩陣加法的CUDA C代碼來展示二維線程模型，因為矩陣是二維的，所以非常適合使用二維線程來處理。我們將在每個線程塊中分別在X維度和Y維度啟動TPB個線程。

#include <stdio.h>

#define N 1024
#define TPB 16

__global__ void MatAdd(int A[N][N], int B[N][N], int 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 a[N][N], b[N][N], c[N][N];

int main()
{
      //int a[N][N], b[N][N], c[N][N];
     int (*dev_a)[N], (*dev_b)[N], (*dev_c)[N];
     for (int y = 0; y < N; ++y)
     {
           for (int x = 0; x < N; ++x)
           {
                 a[x][y] = x + y * N;
                 b[x][y] = x * x + y * N;
           }
     }
     cudaMalloc( &dev_a, N * N * sizeof(int) );
     cudaMalloc( &dev_b, N * N * sizeof(int) );
     cudaMalloc( &dev_c, N * N * sizeof(int) );
     cudaMemcpy(dev_a, a, N * N * sizeof(int), cudaMemcpyHostToDevice );
     cudaMemcpy(dev_b, b, N * N * sizeof(int), cudaMemcpyHostToDevice );
     dim3 threadsPerBlock(TPB, TPB);
     dim3 numBlocks( (N + TPB - 1) / threadsPerBlock.x, (N + TPB -1) / threadsPerBlock.y);
     MatAdd<<<numBlocks, threadsPerBlock>>>(dev_a, dev_b, dev_c);

     cudaMemcpy(c, dev_c, N * N * sizeof(int), cudaMemcpyDeviceToHost );

     int temp = 0;
     for (int y = 0; y < N; ++y)
     {
           for (int x = 0; x < N; ++x)
           {
                 temp = a[x][y] + b[x][y];
                 if (temp != c[x][y])
                 {
                        printf ("Failure at %d  %d\n", x, y);
                 }
           }
     }
     cudaFree(dev_a);
     cudaFree(dev_b);
     cudaFree(dev_c);
     return 0;
}

在這段代碼中，矩陣中每個索引使用一個線程來計算，通過

int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;

來計算當前線程的絕對坐標，該坐標正好對應於矩陣的位置。

需要注意的是，C語言使用太大的數組時，最好把數組定義為全局的，否則受棧的限制，可能會報錯：

Segmentation fault: 11

參考文獻

• 庫克. CUDA並行程序設計. 機械工業出版社, 2014.
• 桑德斯. GPU高性能編程CUDA實戰. 機械工業出版社, 2011.
• CUDA C Programming Guide
• CUDA Toolkit Documentation