本文是 CUDA C 和 C ++的一個系列，它是 CUDA 并行計算平臺的 C / C ++接口。本系列文章假定您熟悉 C 語言編程。我們將針對 Fortran 程序員運行一系列關(guān)于 CUDA Fortran 的文章。這兩個系列將介紹 CUDA 平臺上并行計算的基本概念。從這里起，除非我另有說明，我將用“ CUDA C ”作為“ CUDA C 和 C ++”的速記。 CUDA C 本質(zhì)上是 C / C ++，具有幾個擴(kuò)展，允許使用并行的多個線程在 GPU 上執(zhí)行函數(shù)。

CUDA 編程模型基礎(chǔ)

在我們跳轉(zhuǎn)到 CUDA C 代碼之前， CUDA 新手將從 CUDA 編程模型的基本描述和使用的一些術(shù)語中受益。

CUDA 編程模型是一個異構(gòu)模型，其中使用了 CPU 和 GPU 。在 CUDA 中， host 指的是 CPU 及其存儲器， device 是指 GPU 及其存儲器。在主機(jī)上運行的代碼可以管理主機(jī)和設(shè)備上的內(nèi)存，還可以啟動在設(shè)備上執(zhí)行的函數(shù) kernels 。這些內(nèi)核由許多 GPU 線程并行執(zhí)行。

鑒于 CUDA 編程模型的異構(gòu)性， CUDA C 程序的典型操作序列是：

聲明并分配主機(jī)和設(shè)備內(nèi)存。

初始化主機(jī)數(shù)據(jù)。

將數(shù)據(jù)從主機(jī)傳輸?shù)皆O(shè)備。

執(zhí)行一個或多個內(nèi)核。

將結(jié)果從設(shè)備傳輸?shù)街鳈C(jī)。

記住這個操作序列，讓我們看一個 CUDA C 示例。

第一個 CUDA C 程序

在最近的一篇文章中，我演示了 薩克斯比的六種方法 ，其中包括一個 CUDA C 版本。 SAXPY 代表“單精度 A * X + Y ”，是并行計算的一個很好的“ hello world ”示例。在這篇文章中，我將剖析 CUDA C SAXPY 的一個更完整的版本，詳細(xì)解釋它的作用和原因。完整的 SAXPY 代碼是：

#include 

__global__
void saxpy(int n, float a, float *x, float *y)
{
 int i = blockIdx.x*blockDim.x + threadIdx.x;
 if (i < n) y[i] = a*x[i] + y[i];
}

int main(void)
{
  int N = 1<<20;
  float *x, *y, *d_x, *d_y;
  x = (float*)malloc(N*sizeof(float));
  y = (float*)malloc(N*sizeof(float));

  cudaMalloc(&d_x, N*sizeof(float));?
  cudaMalloc(&d_y, N*sizeof(float));

  for (int i = 0; i < N; i++) {
    x[i] = 1.0f;
    y[i] = 2.0f;
  }

  cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

  // Perform SAXPY on 1M elements
  saxpy<<<(N+255)/256, 256>>>(N, 2.0f, d_x, d_y);

  cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

  float maxError = 0.0f;
  for (int i = 0; i < N; i++)
    maxError = max(maxError, abs(y[i]-4.0f));
  printf("Max error: %f
", maxError);

  cudaFree(d_x);
  cudaFree(d_y);
  free(x);
  free(y);
}

函數(shù)saxpy是在 GPU 上并行運行的內(nèi)核，main函數(shù)是宿主代碼。讓我們從宿主代碼開始討論這個程序。

主機(jī)代碼

main 函數(shù)聲明兩對數(shù)組。

  float *x, *y, *d_x, *d_y;
  x = (float*)malloc(N*sizeof(float));
  y = (float*)malloc(N*sizeof(float));

  cudaMalloc(&d_x, N*sizeof(float));
  cudaMalloc(&d_y, N*sizeof(float));

指針x和y指向以典型方式使用malloc分配的主機(jī)陣列，d_x和d_y數(shù)組指向從CUDA運行時API使用cudaMalloc函數(shù)分配的設(shè)備數(shù)組。CUDA中的主機(jī)和設(shè)備有獨立的內(nèi)存空間，這兩個空間都可以從主機(jī)代碼進(jìn)行管理（CUDAC內(nèi)核也可以在支持它的設(shè)備上分配設(shè)備內(nèi)存）。

然后，主機(jī)代碼初始化主機(jī)數(shù)組。在這里，我們設(shè)置了一個 1 數(shù)組，以及一個 2 數(shù)組。

  for (int i = 0; i < N; i++) {
    x[i] = 1.0f;
    y[i] = 2.0f;
  }

為了初始化設(shè)備數(shù)組，我們只需使用cudaMemcpy將數(shù)據(jù)從x和y復(fù)制到相應(yīng)的設(shè)備數(shù)組d_x和d_y，它的工作方式與標(biāo)準(zhǔn)的 Cmemcpy函數(shù)一樣，只是它采用了第四個參數(shù)，指定了復(fù)制的方向。在本例中，我們使用cudaMemcpyHostToDevice指定第一個（目標(biāo)）參數(shù)是設(shè)備指針，第二個（源）參數(shù)是主機(jī)指針。

  cudaMemcpy(d_x, x, N*sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(d_y, y, N*sizeof(float), cudaMemcpyHostToDevice);

在運行內(nèi)核之后，為了將結(jié)果返回到主機(jī)，我們使用cudaMemcpy和cudaMemcpyDeviceToHost，從d_y指向的設(shè)備數(shù)組復(fù)制到y指向的主機(jī)數(shù)組。

cudaMemcpy(y, d_y, N*sizeof(float), cudaMemcpyDeviceToHost);

啟動內(nèi)核

cord [EZX13 內(nèi)核由以下語句啟動：

saxpy<<<(N+255)/256, 256>>>(N, 2.0, d_x, d_y);

三個 V 形符號之間的信息是 執(zhí)行配置 ，它指示有多少設(shè)備線程并行執(zhí)行內(nèi)核。在 CUDA 中，軟件中有一個線程層次結(jié)構(gòu)，它模仿線程處理器在 GPU 上的分組方式。在 CUDA 編程模型中，我們談到啟動一個 grid 為 螺紋塊 的內(nèi)核。執(zhí)行配置中的第一個參數(shù)指定網(wǎng)格中線程塊的數(shù)量，第二個參數(shù)指定線程塊中的線程數(shù)。

線程塊和網(wǎng)格可以通過為這些參數(shù)傳遞 dim3 （一個由 CUDA 用 x 、 y 和 z 成員定義的簡單結(jié)構(gòu)）值來生成一維、二維或三維的線程塊和網(wǎng)格，但是對于這個簡單的示例，我們只需要一維，所以我們只傳遞整數(shù)。在本例中，我們使用包含 256 個線程的線程塊啟動內(nèi)核，并使用整數(shù)算術(shù)來確定處理數(shù)組（ （N+255）/256 ）的所有 N 元素所需的線程塊數(shù)。

對于數(shù)組中的元素數(shù)不能被線程塊大小平均整除的情況，內(nèi)核代碼必須檢查內(nèi)存訪問是否越界。

清理

完成后，我們應(yīng)該釋放所有分配的內(nèi)存。對于使用 cudaMalloc（） 分配的設(shè)備內(nèi)存，只需調(diào)用 cudaFree（） 。對于主機(jī)內(nèi)存，請像往常一樣使用 free（） 。

cudaFree(d_x);
  cudaFree(d_y);
  free(x);
  free(y);

設(shè)備代碼

現(xiàn)在我們繼續(xù)討論內(nèi)核代碼。

__global__
void saxpy(int n, float a, float *x, float *y)
{
 int i = blockIdx.x*blockDim.x + threadIdx.x;
 if (i < n) y[i] = a*x[i] + y[i];
}

在 CUDA 中，我們使用 __global__ de __global__ 說明符定義諸如 Clara 這樣的內(nèi)核。設(shè)備代碼中定義的變量不需要指定為設(shè)備變量，因為假定它們駐留在設(shè)備上。在這種情況下， n 、 a 和 i 變量將由每個線程存儲在寄存器中，指針 x 和 y 必須是指向設(shè)備內(nèi)存地址空間的指針。這確實是真的，因為當(dāng)我們從宿主代碼啟動內(nèi)核時，我們將 d_x 和 d_y 傳遞給了內(nèi)核。但是，前兩個參數(shù) n 和 a 沒有在主機(jī)代碼中顯式傳輸?shù)皆O(shè)備。因為函數(shù)參數(shù)在 C / C ++中是默認(rèn)通過值傳遞的，所以 CUDA 運行時可以自動處理這些值到設(shè)備的傳輸。 CUDA 運行時 API 的這一特性使得在 GPU 上啟動內(nèi)核變得非常自然和簡單——這幾乎與調(diào)用 C 函數(shù)一樣。

在我們的 saxpy 內(nèi)核中只有兩行。如前所述，內(nèi)核由多個線程并行執(zhí)行。如果我們希望每個線程處理結(jié)果數(shù)組的一個元素，那么我們需要一種區(qū)分和標(biāo)識每個線程的方法。 CUDA 定義變量 blockDim 、 blockIdx 和 threadIdx 。這些預(yù)定義變量的類型為 dim3 ，類似于主機(jī)代碼中的執(zhí)行配置參數(shù)。預(yù)定義變量 blockDim 包含在內(nèi)核啟動的第二個執(zhí)行配置參數(shù)中指定的每個線程塊的維度。預(yù)定義變量 threadIdx 和 blockIdx 分別包含線程塊中線程的索引和網(wǎng)格中的線程塊的索引。表達(dá)式：

    int i = blockDim.x * blockIdx.x + threadIdx.x

生成用于訪問數(shù)組元素的全局索引。我們在這個例子中沒有使用它，但是還有一個 gridDim ，它包含在啟動的第一個執(zhí)行配置參數(shù)中指定的網(wǎng)格維度。

在使用該索引訪問數(shù)組元素之前，將根據(jù)元素的數(shù)量 n 檢查其值，以確保沒有越界內(nèi)存訪問。如果一個數(shù)組中的元素數(shù)不能被線程塊大小平均整除，并且結(jié)果內(nèi)核啟動的線程數(shù)大于數(shù)組大小，則需要進(jìn)行此檢查。內(nèi)核的第二行執(zhí)行 SAXPY 的元素級工作，除了邊界檢查之外，它與 SAXPY 主機(jī)實現(xiàn)的內(nèi)部循環(huán)相同。

if (i < n) y[i] = a*x[i] + y[i];

編譯和運行代碼

CUDA C 編譯器 nvcc 是 NVIDIA CUDA 工具箱 的一部分。為了編譯我們的 SAXPY 示例，我們將代碼保存在一個擴(kuò)展名為。 cu 的文件中，比如說 saxpy.cu 。然后我們可以用 nvcc 編譯它。

nvcc -o saxpy saxpy.cu

然后我們可以運行代碼：

% ./saxpy
Max error: 0.000000

總結(jié)與結(jié)論

通過對 SAXPY 的一個簡單的 CUDA C 實現(xiàn)的演練，您現(xiàn)在了解了編程 CUDA C 的基本知識。將 C 代碼“移植”到 CUDA C 只需要幾個 C 擴(kuò)展：設(shè)備內(nèi)核函數(shù)的 __global__ de Clara 說明符；啟動內(nèi)核時使用的執(zhí)行配置；內(nèi)置的設(shè)備變量 blockDim 、 blockIdx 和 threadIdx 用來識別和區(qū)分并行執(zhí)行內(nèi)核的 GPU 線程。

異類 CUDA 編程模型的一個優(yōu)點是，將現(xiàn)有代碼從 C 移植到 CUDA C 可以逐步完成，一次只能移植一個內(nèi)核。

在本系列的下一篇文章中，我們將研究一些性能度量和度量。

關(guān)于作者

Mark Harris 是 NVIDIA 杰出的工程師，致力于 RAPIDS 。 Mark 擁有超過 20 年的 GPUs 軟件開發(fā)經(jīng)驗，從圖形和游戲到基于物理的模擬，到并行算法和高性能計算。當(dāng)他還是北卡羅來納大學(xué)的博士生時，他意識到了一種新生的趨勢，并為此創(chuàng)造了一個名字： GPGPU （圖形處理單元上的通用計算）。

審核編輯：郭婷