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