深入浅出CUDA

“CUDA 是 NVIDIA 的 GPGPU 模型，它使用 C 语言为基础，可以直接以大多数 人熟悉的 C 语言，写出在显示芯片上执行的程序，而不需要去学习特定的显示 芯片的指令或是特殊的结构。” CUDA 是什么？能吃吗？ 编者注：NVIDIA 的 GeFoce 8800GTX 发布后，它的通用计算架构 CUDA 经过一年多的推广后，现 在已经在有相当多的论文发表，在商业应用软件等方面也初步出现了视频编解码、金融、地质勘 探、科学计算等领域的产品，是时候让我们对其作更深一步的了解。为了让大家更容易了解 CUDA， 我们征得 Hotball 的本人同意，发表他最近亲自撰写的本文。这篇文章的特点是深入浅出，也包 含了 hotball 本人编写一些简单 CUDA 程序的亲身体验，对于希望了解 CUDA 的读者来说是非常不 错的入门文章，PCINLIFE 对本文的发表没有作任何的删减，主要是把一些台湾的词汇转换成大 陆的词汇以及作了若干"编者注"的注释。 现代的显示芯片已经具有高度的可程序化能力，由于显示芯片通常具有相当高的内存带宽， 以及大量的执行单元，因此开始有利用显示芯片来帮助进行一些计算工作的想法，即 GPGPU。 CUDA 即是 NVIDIA 的 GPGPU 模型。 NVIDIA 的新一代显示芯片，包括 GeForce 8 系列及更新的显示芯片都支持 CUDA。NVIDIA 免费提供 CUDA 的开发工具（包括 Windows 版本和 Linux 版本）、程序范例、文件等等， 可以在 CUDA Zone 下载。 GPGPU 的优缺点 使用显示芯片来进行运算工作，和使用 CPU 相比，主要有几个好处： 1. 显示芯片通常具有更大的内存带宽。例如，NVIDIA 的 GeForce 8800GTX 具有 超过 50GB/s 的内存带宽，而目前高阶 CPU 的内存带宽则在 10GB/s 左右。 2. 显示芯片具有更大量的执行单元。例如 GeForce 8800GTX 具有 128 个 "stream processors"，频率为 1.35GHz。CPU 频率通常较高，但是执行单元的数目 则要少得多。 3. 和高阶 CPU 相比，显卡的价格较为低廉。例如目前一张 GeForce 8800GT 包 括 512MB 内存的价格，和一颗 2.4GHz 四核心 CPU 的价格相若。 当然，使用显示芯片也有它的一些缺点： 1. 显示芯片的运算单元数量很多，因此对于不能高度并行化的工作，所能带来的 帮助就不大。 2. 显示芯片目前通常只支持 32 bits 浮点数，且多半不能完全支持 IEEE 754 规 格， 有些运算的精确度可能较低。目前许多显示芯片并没有分开的整数运算单元， 因此整数运算的效率较差。 3. 显示芯片通常不具有分支预测等复杂的流程控制单元，因此对于具有高度分支 的程序，效率会比较差。 4. 目前 GPGPU 的程序模型仍不成熟，也还没有公认的标准。例如 NVIDIA 和 AMD/ATI 就有各自不同的程序模型。 整体来说，显示芯片的性质类似 stream processor，适合一次进行大量相同的工作。CPU 则 比较有弹性，能同时进行变化较多的工作。 CUDA 架构 CUDA 是 NVIDIA 的 GPGPU 模型，它使用 C 语言为基础，可以直接以大多数人熟悉的 C 语 言，写出在显示芯片上执行的程序，而不需要去学习特定的显示芯片的指令或是特殊的结构。 在 CUDA 的架构下，一个程序分为两个部份：host 端和 device 端。Host 端是指在 CPU 上 执行的部份，而 device 端则是在显示芯片上执行的部份。Device 端的程序又称为 "kernel"。通常 host 端程序会将数据准备好后，复制到显卡的内存中，再由显示芯片执行 device 端程序，完成后再由 host 端程序将结果从显卡的内存中取回。 由于 CPU 存取显卡内存时只能透过 PCI Express 接口，因此速度较慢（PCI Express x16 的 理论带宽是双向各 4GB/s），因此不能太常进行这类动作，以免降低效率。 在 CUDA 架构下，显示芯片执行时的最小单位是 thread。数个 thread 可以组成一个 block。 一个 block 中的 thread 能存取同一块共享的内存，而且可以快速进行同步的动作。 每一个 block 所能包含的 thread 数目是有限的。不过，执行相同程序的 block，可以组 成 grid。不同 block 中的 thread 无法存取同一个共享的内存，因此无法直接互通或进行 同步。因此，不同 block 中的 thread 能合作的程度是比较低的。不过，利用这个模式， 可以让程序不用担心显示芯片实际上能同时执行的 thread 数目限制。例如，一个具有很少 量执行单元的显示芯片，可能会把各个 block 中的 thread 顺序执行，而非同时执行。不 同的 grid 则可以执行不同的程序（即 kernel）。 Grid、block 和 thread 的关系，如下图所示： 每个 thread 都有自己的一份 register 和 local memory 的空间。同一个 block 中的每 个 thread 则有共享的一份 share memory。此外，所有的 thread（包括不同 block 的 thread）都共享一份 global memory、constant memory、和 texture memory。不同的 grid 则有各自的 global memory、constant memory 和 texture memory。这些不同的内存的差 别，会在之后讨论。 执行模式 由于显示芯片大量并行计算的特性，它处理一些问题的方式，和一般 CPU 是不同的。主要 的特点包括： 1. 内存存取 latency 的问题：CPU 通常使用 cache 来减少存取主内存的次数， 以避免内存 latency 影响到执行效率。显示芯片则多半没有 cache（或很小），而 利用并行化执行的方式来隐藏内存的 latency（即，当第一个 thread 需要等待内 存读取结果时，则开始执行第二个 thread，依此类推）。 2. 分支指令的问题：CPU 通常利用分支预测等方式来减少分支指令造成的 pipeline bubble。显示芯片则多半使用类似处理内存 latency 的方式。不过，通 常显示芯片处理分支的效率会比较差。 因此，最适合利用 CUDA 处理的问题，是可以大量并行化的问题，才能有效隐藏内存的 latency，并有效利用显示芯片上的大量执行单元。使用 CUDA 时，同时有上千个 thread 在 执行是很正常的。因此，如果不能大量并行化的问题，使用 CUDA 就没办法达到最好的效率 了。 CUDA Toolkit 的安装 目前 NVIDIA 提供的 CUDA Toolkit（可从这里下载）支持 Windows （32 bits 及 64 bits 版本）及许多不同的 Linux 版本。 CUDA Toolkit 需要配合 C/C++ compiler。在 Windows 下，目前只支持 Visual Studio 7.x 及 Visual Studio 8（包括免费的 Visual Studio C++ 2005 Express）。Visual Studio 6 和 gcc 在 Windows 下是不支援的。在 Linux 下则只支援 gcc。 这里简单介绍一下在 Windows 下设定并使用 CUDA 的方式。 下载及安装 在 Windows 下，CUDA Toolkit 和 CUDA SDK 都是由安装程序的形式安装的。CUDA Toolkit 包括 CUDA 的基本工具，而 CUDA SDK 则包括许多范例程序以及链接库。基本上要写 CUDA 的程序，只需要安装 CUDA Toolkit 即可。不过 CUDA SDK 仍值得安装，因为里面的许多范 例程序和链接库都相当有用。 CUDA Toolkit 安装完后，预设会安装在 C:\CUDA 目录里。其中包括几个目录： • bin -- 工具程序及动态链接库 • doc -- 文件 • include -- header 檔 • lib -- 链接库档案 • open64 -- 基于 Open64 的 CUDA compiler • src -- 一些原始码 安装程序也会设定一些环境变量，包括： • CUDA_BIN_PATH -- 工具程序的目录，默认为 C:\CUDA\bin • CUDA_INC_PATH -- header 文件的目录，默认为 C:\CUDA\inc • CUDA_LIB_PATH -- 链接库文件的目录，默认为 C:\CUDA\lib 在 Visual Studio 中使用 CUDA CUDA 的主要工具是 nvcc，它会执行所需要的程序，将 CUDA 程序代码编译成执行档 (或 object 檔) 。在 Visual Studio 下，我们透过设定 custom build tool 的方式，让 Visual Studio 会自动执行 nvcc。 这里以 Visual Studio 2005 为例： 1. 首先，建立一个 Win32 Console 模式的 project（在 Application Settings 中记得勾选 Empty project），并新增一个档案，例如 main.cu。 2. 在 main.cu 上右键单击，并选择 Properties。点选 General，确定 Tool 的 部份是选择 Custom Build Tool。 3. 选择 Custom Build Step，在 Command Line 使用以下设定： o Release 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccbin "$(VCInstallDir)bin" -c -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)\$(InputName).obj $(InputFileName) o Debug 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccbin "$(VCInstallDir)bin" -c -D_DEBUG -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)\$(InputName).obj $(InputFileName) 4. 如果想要使用软件仿真的模式，可以新增两个额外的设定： o EmuRelease 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccbin "$(VCInstallDir)bin" -deviceemu -c -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)\$(InputName).obj $(InputFileName) o EmuDebug 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccbin "$(VCInstallDir)bin" -deviceemu -c -D_DEBUG -DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o $(ConfigurationName)\$(InputName).obj $(InputFileName) 5. 对所有的配置文件，在 Custom Build Step 的 Outputs 中加入 $(ConfigurationName)\$(InputName).obj。 6. 选择 project，右键单击选择 Properties，再点选 Linker。对所有的配置文 件修改以下设定： o General/Enable Incremental Linking：No o General/Additional Library Directories：$(CUDA_LIB_PATH) o Input/Additional Dependencies：cudart.lib 这样应该就可以直接在 Visual Studio 的 IDE 中，编辑 CUDA 程序后，直接 build 以及 执行程序了。 第一个 CUDA 程序 CUDA 目前有两种不同的 API：Runtime API 和 Driver API，两种 API 各有其适用的范围。 由于 runtime API 较容易使用，一开始我们会以 runetime API 为主。 CUDA 的初始化 首先，先建立一个档案 first_cuda.cu。如果是使用 Visual Studio 的话，则请先按照这 里的设定方式设定 project。 要使用 runtime API 的时候，需要 include cuda_runtime.h。所以，在程序的最前面， 加上 #include #include 接下来是一个 InitCUDA 函式，会呼叫 runtime API 中，有关初始化 CUDA 的功能： bool InitCUDA() { int count; cudaGetDeviceCount(&count); if(count == 0) { fprintf(stderr, "There is no device.\n"); return false; } int i; for(i = 0; i < count; i++) { cudaDeviceProp prop; if(cudaGetDeviceProperties(∝, i) == cudaSuccess) { if(prop.major >= 1) { break; } } } if(i == count) { fprintf(stderr, "There is no device supporting CUDA 1.x.\n"); return false; } cudaSetDevice(i); return true; } 这个函式会先呼叫 cudaGetDeviceCount 函式，取得支持 CUDA 的装置的数目。如果系 统上没有支持 CUDA 的装置，则它会传回 1，而 device 0 会是一个仿真的装置，但不支持 CUDA 1.0 以上的功能。所以，要确定系统上是否有支持 CUDA 的装置，需要对每个 device 呼叫 cudaGetDeviceProperties 函式，取得装置的各项数据，并判断装置支持的 CUDA 版本（prop.major 和 prop.minor 分别代表装置支持的版本号码，例如 1.0 则 prop.major 为 1 而 prop.minor 为 0）。 透过 cudaGetDeviceProperties 函式可以取得许多数据，除了装置支持的 CUDA 版本 之外，还有装置的名称、内存的大小、最大的 thread 数目、执行单元的频率等等。详情可 参考 NVIDIA 的 CUDA Programming Guide。 在找到支持 CUDA 1.0 以上的装置之后，就可以呼叫 cudaSetDevice 函式，把它设为目 前要使用的装置。 最后是 main 函式。在 main 函式中我们直接呼叫刚才的 InitCUDA 函式，并显示适当的讯 息： int main() { if(!InitCUDA()) { return 0; } printf("CUDA initialized.\n"); return 0; } 这样就可以利用 nvcc 来 compile 这个程序了。使用 Visual Studio 的话，若按照先前的 设定方式，可以直接 Build Project 并执行。 nvcc 是 CUDA 的 compile 工具，它会将 .cu 檔拆解出在 GPU 上执行的部份，及在 host 上执行的部份，并呼叫适当的程序进行 compile 动作。在 GPU 执行的部份会透过 NVIDIA 提供的 compiler 编译成中介码，而 host 执行的部份则会透过系统上的 C++ compiler 编 译（在 Windows 上使用 Visual C++ 而在 Linux 上使用 gcc）。 编译后的程序，执行时如果系统上有支持 CUDA 的装置，应该会显示 CUDA initialized. 的 讯息，否则会显示相关的错误讯息。 利用 CUDA 进行运算 到目前为止，我们的程序并没有做什么有用的工作。所以，现在我们加入一个简单的动作， 就是把一大堆数字，计算出它的平方和。 首先，把程序最前面的 include 部份改成： #include #include #include #define DATA_SIZE 1048576 int data[DATA_SIZE]; 并加入一个新函式 GenerateNumbers： void GenerateNumbers(int *number, int size) { for(int i = 0; i < size; i++) { number[i] = rand() % 10; } } 这个函式会产生一大堆 0 ~ 9 之间的随机数。 要利用 CUDA 进行计算之前，要先把数据复制到显卡内存中，才能让显示芯片使用。因此， 需要取得一块适当大小的显卡内存，再把产生好的数据复制进去。在 main 函式中加入： GenerateNumbers(data, DATA_SIZE); int* gpudata, *result; cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE); cudaMalloc((void**) &result, sizeof(int)); cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice); 上面这段程序会先呼叫 GenerateNumbers 产生随机数，并呼叫 cudaMalloc 取得一块显 卡内存（result 则是用来存取计算结果，在稍后会用到），并透过 cudaMemcpy将产生 的随机数复制到显卡内存中。cudaMalloc 和 cudaMemcpy 的用法和一般的 malloc 及 memcpy 类似，不过 cudaMemcpy 则多出一个参数，指示复制内存的方向。在这里因为是从 主内存复制到显卡内存，所以使用 cudaMemcpyHostToDevice。如果是从显卡内存到主 内存，则使用 cudaMemcpyDeviceToHost。这在之后会用到。 接下来是要写在显示芯片上执行的程序。在 CUDA 中，在函式前面加上 __global__表示 这个函式是要在显示芯片上执行的。因此，加入以下的函式： __global__ static void sumOfSquares(int *num, int* result) { int sum = 0; int i; for(i = 0; i < DATA_SIZE; i++) { sum += num[i] * num[i]; } *result = sum; } 在显示芯片上执行的程序有一些限制，例如它不能有传回值。其它的限制会在之后提到。 接下来是要让 CUDA 执行这个函式。在 CUDA 中，要执行一个函式，使用以下的语法： 函式名称<<>>(参数...); 呼叫完后，还要把结果从显示芯片复制回主内存上。在 main 函式中加入以下的程序： sumOfSquares<<<1, 1, 0>>>(gpudata, result); int sum; cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); printf("sum: %d\n", sum); 因为这个程序只使用一个 thread，所以 block 数目、thread 数目都是 1。我们也没有使 用到任何 shared memory，所以设为 0。编译后执行，应该可以看到执行的结果。 为了确定执行的结果正确，我们可以加上一段以 CPU 执行的程序代码，来验证结果： sum = 0; for(int i = 0; i < DATA_SIZE; i++) { sum += data[i] * data[i]; } printf("sum (CPU): %d\n", sum); 编译后执行，确认两个结果相同。 计算运行时间 CUDA 提供了一个 clock 函式，可以取得目前的 timestamp，很适合用来判断一段程序执行 所花费的时间（单位为 GPU 执行单元的频率）。这对程序的优化也相当有用。要在我们的 程序中记录时间，把 sumOfSquares 函式改成： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { int sum = 0; int i; clock_t start = clock(); for(i = 0; i < DATA_SIZE; i++) { sum += num[i] * num[i]; } *result = sum; *time = clock() - start; } 把 main 函式中间部份改成： int* gpudata, *result; clock_t* time; cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE); cudaMalloc((void**) &result, sizeof(int)); cudaMalloc((void**) &time, sizeof(clock_t)); cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice); sumOfSquares<<<1, 1, 0>>>(gpudata, result, time); int sum; clock_t time_used; cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost); cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); printf("sum: %d time: %d\n", sum, time_used); 编译后执行，就可以看到执行所花费的时间了。 如果计算实际运行时间的话，可能会注意到它的执行效率并不好。这是因为我们的程序并没 有利用到 CUDA 的主要的优势，即并行化执行。在下一段文章中，会讨论如何进行优化的动 作。 改良第一个 CUDA 程序 在上一篇文章中，我们做了一个计算一大堆数字的平方和的程序。不过，我们也提到这个程 序的执行效率并不理想。当然，实际上来说，如果只是要做计算平方和的动作，用 CPU 做 会比用 GPU 快得多。这是因为平方和的计算并不需要太多运算能力，所以几乎都是被内存 带宽所限制。因此，光是把数据复制到显卡内存上的这个动作，所需要的时间，可能已经和 直接在 CPU 上进行计算差不多了。 不过，如果进行平方和的计算，只是一个更复杂的计算过程的一部份的话，那么当然在 GPU 上计算还是有它的好处的。而且，如果数据已经在显卡内存上（例如在 GPU 上透过某种算 法产生），那么，使用 GPU 进行这样的运算，还是会比较快的。 刚才也提到了，由于这个计算的主要瓶颈是内存带宽，所以，理论上显卡的内存带宽是相当 大的。这里我们就来看看，倒底我们的第一个程序，能利用到多少内存带宽。 程序的并行化 我们的第一个程序，并没有利用到任何并行化的功能。整个程序只有一个 thread。在 GeForce 8800GT 上面，在 GPU 上执行的部份（称为 "kernel"）大约花费 640M 个频率。 GeForce 8800GT 的执行单元的频率是 1.5GHz，因此这表示它花费了约 0.43 秒的时间。1M 个 32 bits 数字的数据量是 4MB，因此，这个程序实际上使用的内存带宽，只有 9.3MB/s 左 右！这是非常糟糕的表现。 为什么会有这样差的表现呢？这是因为 GPU 的架构特性所造成的。在 CUDA 中，一般的数 据复制到的显卡内存的部份，称为 global memory。这些内存是没有 cache 的，而且，存 取 global memory 所需要的时间（即 latency）是非常长的，通常是数百个 cycles。由于 我们的程序只有一个 thread，所以每次它读取 global memory 的内容，就要等到实际读取 到数据、累加到 sum 之后，才能进行下一步。这就是为什么它的表现会这么的差。 由于 global memory 并没有 cache，所以要避开巨大的 latency 的方法，就是要利用大量 的 threads。假设现在有大量的 threads 在同时执行，那么当一个 thread 读取内存，开 始等待结果的时候，GPU 就可以立刻切换到下一个 thread，并读取下一个内存位置。因此， 理想上当 thread 的数目够多的时候，就可以完全把 global memory 的巨大 latency 隐藏 起来了。 要怎么把计算平方和的程序并行化呢？最简单的方法，似乎就是把数字分成若干组，把各组 数字分别计算平方和后，最后再把每组的和加总起来就可以了。一开始，我们可以把最后加 总的动作，由 CPU 来进行。 首先，在 first_cuda.cu 中，在 #define DATA_SIZE 的后面增加一个 #define，设定 thread 的数目： #define DATA_SIZE 1048576 #define THREAD_NUM 256 接着，把 kernel 程序改成： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { const int tid = threadIdx.x; const int size = DATA_SIZE / THREAD_NUM; int sum = 0; int i; clock_t start; if(tid == 0) start = clock(); for(i = tid * size; i < (tid + 1) * size; i++) { sum += num[i] * num[i]; } result[tid] = sum; if(tid == 0) *time = clock() - start; } 程序里的 threadIdx 是 CUDA 的一个内建的变量，表示目前的 thread 是第几个 thread （由 0 开始计算）。以我们的例子来说，会有 256 个 threads，所以同时会有 256 个 sumOfSquares 函式在执行，但每一个的 threadIdx.x 则分别会是 0 ~ 255。利用这 个变量，我们就可以让每一份函式执行时，对整个数据不同的部份计算平方和。另外，我们 也让计算时间的动作，只在 thread 0（即 threadIdx.x = 0 的时候）进行。 同样的，由于会有 256 个计算结果，所以原来存放 result 的内存位置也要扩大。把 main 函式中的中间部份改成： int* gpudata, *result; clock_t* time; cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE); cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM); cudaMalloc((void**) &time, sizeof(clock_t)); cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice); sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time); int sum[THREAD_NUM]; clock_t time_used; cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM, cudaMemcpyDeviceToHost); cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); cudaFree(time); 可以注意到我们在呼叫 sumOfSquares 函式时，指定 THREAD_NUM 为 thread 的数目。 最后，在 CPU 端把计算好的各组数据的平方和进行加总： int final_sum = 0; for(int i = 0; i < THREAD_NUM; i++) { final_sum += sum[i]; } printf("sum: %d time: %d\n", final_sum, time_used); final_sum = 0; for(int i = 0; i < DATA_SIZE; i++) { sum += data[i] * data[i]; } printf("sum (CPU): %d\n", final_sum); 编译后执行，确认结果和原来相同。 这个版本的程序，在 GeForce 8800GT 上执行，只需要约 8.3M cycles，比前一版程序快了 77 倍！这就是透过大量 thread 来隐藏 latency 所带来的效果。 不过，如果计算一下它使用的内存带宽，就会发现其实仍不是很理想，大约只有 723MB/s 而 已。这和 GeForce 8800GT 所具有的内存带宽是很大的差距。为什么会这样呢？ 内存的存取模式 显卡上的内存是 DRAM，因此最有效率的存取方式，是以连续的方式存取。前面的程序，虽 然看起来是连续存取内存位置（每个 thread 对一块连续的数字计算平方和），但是我们要 考虑到实际上 thread 的执行方式。前面提过，当一个 thread 在等待内存的数据时，GPU 会 切换到下一个 thread。也就是说，实际上执行的顺序是类似 thread 0 -> thread 1 -> thread 2 -> ... 因此，在同一个 thread 中连续存取内存，在实际执行时反而不是连续了。要让实际执行结 果是连续的存取，我们应该要让 thread 0 读取第一个数字，thread 1 读取第二个数字… 依此类推。所以，我们可以把 kernel 程序改成如下： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { const int tid = threadIdx.x; int sum = 0; int i; clock_t start; if(tid == 0) start = clock(); for(i = tid; i < DATA_SIZE; i += THREAD_NUM) { sum += num[i] * num[i]; } result[tid] = sum; if(tid == 0) *time = clock() - start; } 编译后执行，确认结果相同。 仅仅是这样简单的修改，实际执行的效率就有很大的差别。在 GeForce 8800GT 上，上面的 程序执行需要的频率是 2.6M cycles，又比前一版程序快了三倍。不过，这样仍只有 2.3GB/s 的带宽而已。 这是因为我们使用的 thread 数目还是不够多的原因。理论上 256 个 threads 最多只能隐 藏 256 cycles 的 latency。但是 GPU 存取 global memory 时的 latency 可能高达 500 cycles 以上。如果增加 thread 数目，就可以看到更好的效率。例如，可以把 THREAD_NUM 改成 512。在 GeForce 8800GT 上，这可以让执行花费的时间减少到 1.95M cycles。有些改进，但是仍不够大。不幸的是，目前 GeForce 8800GT 一个 block 最多只 能有 512 个 threads，所以不能再增加了，而且，如果 thread 数目增加太多，那么在 CPU 端要做的最后加总工作也会变多。 更多的并行化 前面提到了 block。在之前介绍呼叫 CUDA 函式时，也有提到 "block 数目" 这个参数。到 目前为止，我们都只使用一个 block。究竟 block 是什么呢？ 在 CUDA 中，thread 是可以分组的，也就是 block。一个 block 中的 thread，具有一个 共享的 shared memory，也可以进行同步工作。不同 block 之间的 thread 则不行。在我 们的程序中，其实不太需要进行 thread 的同步动作，因此我们可以使用多个 block 来进 一步增加 thread 的数目。 首先，在 #define DATA_SIZE 的地方，改成如下： #define DATA_SIZE 1048576 #define BLOCK_NUM 32 #define THREAD_NUM 256 这表示我们会建立 32 个 blocks，每个 blocks 有 256 个 threads，总共有 32*256 = 8192 个 threads。 接着，我们把 kernel 部份改成： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { const int tid = threadIdx.x; const int bid = blockIdx.x; int sum = 0; int i; if(tid == 0) time[bid] = clock(); for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) { sum += num[i] * num[i]; } result[bid * THREAD_NUM + tid] = sum; if(tid == 0) time[bid + BLOCK_NUM] = clock(); } blockIdx.x 和 threadIdx.x 一样是 CUDA 内建的变量，它表示的是目前的 block 编 号。另外，注意到我们把计算时间的方式改成每个 block 都会记录开始时间及结束时间。 main 函式部份，修改成： int* gpudata, *result; clock_t* time; cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE); cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM * BLOCK_NUM); cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2); cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice); sumOfSquares<<>>(gpudata, result, time); int sum[THREAD_NUM * BLOCK_NUM]; clock_t time_used[BLOCK_NUM * 2]; cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM * BLOCK_NUM, cudaMemcpyDeviceToHost); cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); cudaFree(time); int final_sum = 0; for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) { final_sum += sum[i]; } clock_t min_start, max_end; min_start = time_used[0]; max_end = time_used[BLOCK_NUM]; for(int i = 1; i < BLOCK_NUM; i++) { if(min_start > time_used[i]) min_start = time_used[i]; if(max_end < time_used[i + BLOCK_NUM]) max_end = time_used[i + BLOCK_NUM]; } printf("sum: %d time: %d\n", final_sum, max_end - min_start); 基本上我们只是把 result 的大小变大，并修改计算时间的方式，把每个 block 最早的开 始时间，和最晚的结束时间相减，取得总运行时间。 这个版本的程序，执行的时间减少很多，在 GeForce 8800GT 上只需要约 150K cycles，相 当于 40GB/s 左右的带宽。不过，它在 CPU 上执行的部份，需要的时间加长了（因为 CPU 现 在需要加总 8192 个数字）。为了避免这个问题，我们可以让每个 block 把自己的每个 thread 的计算结果进行加总。 Thread 的同步 前面提过，一个 block 内的 thread 可以有共享的内存，也可以进行同步。我们可以利用 这一点，让每个 block 内的所有 thread 把自己计算的结果加总起来。把 kernel 改成如 下： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { extern __shared__ int shared[]; const int tid = threadIdx.x; const int bid = blockIdx.x; int i; if(tid == 0) time[bid] = clock(); shared[tid] = 0; for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) { shared[tid] += num[i] * num[i]; } __syncthreads(); if(tid == 0) { for(i = 1; i < THREAD_NUM; i++) { shared[0] += shared[i]; } result[bid] = shared[0]; } if(tid == 0) time[bid + BLOCK_NUM] = clock(); } 利用 __shared__ 声明的变量表示这是 shared memory，是一个 block 中每个 thread 都 共享的内存。它会使用在 GPU 上的内存，所以存取的速度相当快，不需要担心 latency 的 问题。 __syncthreads() 是一个 CUDA 的内部函数，表示 block 中所有的 thread 都要同步到 这个点，才能继续执行。在我们的例子中，由于之后要把所有 thread 计算的结果进行加总， 所以我们需要确定每个 thread 都已经把结果写到 shared[tid] 里面了。 接下来，把 main 函式的一部份改成： int* gpudata, *result; clock_t* time; cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE); cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM); cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2); cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice); sumOfSquares<<>>(gpudata, result, time); int sum[BLOCK_NUM]; clock_t time_used[BLOCK_NUM * 2]; cudaMemcpy(∑, result, sizeof(int) * BLOCK_NUM, cudaMemcpyDeviceToHost); cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); cudaFree(time); int final_sum = 0; for(int i = 0; i < BLOCK_NUM; i++) { final_sum += sum[i]; } 可以注意到，现在 CPU 只需要加总 BLOCK_NUM 也就是 3