为了更好地说明如何在Linux实现C++与CUDA的混合编程，我接下来将以实现对一个矩阵的每一个元素的取模运算。
1. 头文件和文件形式
要在C++中编写CUDA代码，需要引入头文件：
#include "cuda_runtime.h" #include "cublas_v2.h" #include "device_launch_parameters.h"
然后将文件后缀由.cpp改为.cu
2. 搭建程序框架
一个CUDA程序的框架类似于普通C++程序，下面我们来展示其中一种简单的形式：
// mod_test.cu #include "cuda_runtime.h" #include "cublas_v2.h" #include "device_launch_parameters.h" // 自己实现的CUDA取模运算函数 __global__ void mod_kernel(args) { // ... } // 进行矩阵元素取模运算 void mat_mod_cu(args){ // ... // 调用自己实现的CUDA取模运算函数 mod_kernel <<< gridSize, blockSize >>> (args); // ... } int main(void){ // 定义cuda消息处理器 cublasHandle_t cuHandle; cublasStatus_t status = cublasCreate(&cuHandle); // ... // 进行矩阵元素取模运算 mat_mod_cu(args); // ... // 销毁cuda消息处理器 cublasDestroy(cuHandle); return 0; }
在上面的框架中：
__global__ void mod_kernel(args)
是自己编写的要在GPU中执行的代码，__global__指明该函数是kernel函数，调用时用<<< >>>语法配置并行化参数；返回值是void的；mod_kernel是CUDA函数名；args是参数。
接下来我们基于此框架实现CUDA取模运算。
2.1. main() 函数
为了获得良好的封装性，main()函数需要向mat_mod()传入需要取模的矩阵、其行列数、取模值、cuda控制符，那么可以在main()函数中编写如下代码：
int main(void){ // 1.定义cuda消息处理器 cublasHandle_t cuHandle; cublasStatus_t status = cublasCreate(&cuHandle); // 2.定义要取模的矩阵 long q = 3; // 取模数 long rowSize = 4; // 矩阵行数 long colSize = 6; // 矩阵列数 long** M1 = uniformMat(rowSize, colSize, 5, 9); // 生成一个rowSize*colSize的随机矩阵，元素值在[5, 9]内随机选取 // 3.进行矩阵元素取模运算 long** M2 = mat_mod_cu(M1, rowSize, colSize, q, cuHandle); // 4.输出原矩阵M1和取模后的矩阵M2 cout << "M1: " << endl; for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { cout << M1[i][j] << " "; } cout << endl; } cout << "M2: " << endl; for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { cout << M2[i][j] << " "; } cout << endl; } // 5.销毁cuda消息处理器 cublasDestroy(cuHandle); return 0; }
uniformMat()函数是用于生成随机矩阵的函数，它的实现如下：
long** uniformMat(long rowSize, long colSize, long minValue, long maxValue) { long** mat = new long* [rowSize]; for (long i = 0; i < rowSize; i++) mat[i] = new long[colSize]; srand((unsigned)time(NULL)); for (long i = 0; i < rowSize; i++) { for (long j = 0; j < colSize; j++) { mat[i][j] = (long)(rand() % (maxValue - minValue + 1)) + minValue; } } return mat; }
2.2. mat_mod_cu() 函数
mat_mod_cu()函数主要复制将矩阵存入GPU内存中，然后调用cuda函数进行运算，再将运算结果从GPU中调出。它的定义如下：
long** mat_mod_cu(long** M1, long rowSize, long colSize, long q, cublasHandle_t cuHandle) { // 1.定义结果矩阵，用于返回 long** M2 = new long* [rowSize]; for (long i = 0; i < rowSize; i++) M2[i] = new long[colSize]; // 2.分配CPU资源 double* h_M1 = (double*)malloc(rowSize * colSize * sizeof(double)); double* h_M2 = (double*)malloc(rowSize * colSize * sizeof(double)); // 初始化h_M1 for (long i = 0; i < rowSize; i++) { for (long j = 0; j < colSize; j++) { h_M1[i * colSize + j] = (double)M1[i][j]; } } // 3.分配GPU资源 double* d_M1; double* d_M2; cudaMalloc((void**)&d_M1, rowSize * colSize * sizeof(double)); cudaMalloc((void**)&d_M2, rowSize * colSize * sizeof(double)); // 将CPU数据拷贝到GPU上 cudaMemcpy(d_M1, h_M1, rowSize * colSize * sizeof(double), cudaMemcpyHostToDevice); // 4.定义kernel的执行配置 int blockSize; int minGridSize; int gridSize; // 获取GPU的信息，并配置最优参数 cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, mod_kernel, 0, rowSize * colSize); gridSize = (rowSize * colSize + blockSize - 1) / blockSize; // 5.执行核函数 // 取模 mod_kernel <<< gridSize, blockSize >>> (d_M1, d_M2, q, rowSize*colSize); // 6.将GPU数据拷贝到CPU上 cudaMemcpy(h_M2, d_M2, rowSize * colSize * sizeof(double), cudaMemcpyDeviceToHost); // 7.赋值给结果矩阵 for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { M2[i][j] = static_cast<long>(h_M2[i * colSize + j]); } } // 8.清理掉使用过的内存 free(h_M1); free(h_M2); cudaFree(d_M1); cudaFree(d_M2); return M2; }
CUDA对矩阵的运算，需要先把矩阵展平为一维向量，因此也就有了第二步，分配CPU资源。然后在把这个一维向量传入GPU中，也就是第三步。之后再对位于GPU中的矩阵进行运算，在运算之前需要指定一些用于并行计算的参数，在第四步中已经设置了自动获取，因此不再需要手动配置。这样就可以执行核函数进行取模运算了。运算完毕之后从GPU内取出数据到CPU，再展开为二维矩阵，返回，就行了。
2.3. mod_kernel() 函数
mod_kernel() 是最终在GPU中执行计算的函数，它的定义如下：
__global__ void mod_kernel(double* d_M1, double* d_M2, long q, int n) { // index int idx = blockIdx.x * blockDim.x + threadIdx.x; long x; if (idx < n){ x = d_M1[idx]; d_M2[idx] = x % q; } }
由代码可见，核函数并不一定需要for循环来遍历执行整个矩阵，这是因为在调用该核函数时指定了gridSize和blockSize，使得矩阵在GPU内是以并行的形式执行的。只需要明确遍历完矩阵内每一个元素即可。
如果想在核函数中输出信息，那么可以用printf，而不是cout。
3. 完整代码及运算结果
完整代码：
/* mod_test.cu */ #include "cuda_runtime.h" #include "cublas_v2.h" #include "device_launch_parameters.h" #include <iostream> #include <stdio.h> using namespace std; __global__ void mod_kernel(double* d_M1, double* d_M2, long q, int n) { // index int idx = blockIdx.x * blockDim.x + threadIdx.x; long x; if (idx < n){ x = d_M1[idx]; // printf("x: %d\n", x) d_M2[idx] = x % q; } } long** mat_mod_cu(long** M1, long rowSize, long colSize, long q, cublasHandle_t cuHandle) { // 1.定义结果矩阵，用于返回 long** M2 = new long* [rowSize]; for (long i = 0; i < rowSize; i++) M2[i] = new long[colSize]; // 2.分配CPU资源 double* h_M1 = (double*)malloc(rowSize * colSize * sizeof(double)); double* h_M2 = (double*)malloc(rowSize * colSize * sizeof(double)); // 初始化h_M1 for (long i = 0; i < rowSize; i++) { for (long j = 0; j < colSize; j++) { h_M1[i * colSize + j] = (double)M1[i][j]; } } // 3.分配GPU资源 double* d_M1; double* d_M2; cudaMalloc((void**)&d_M1, rowSize * colSize * sizeof(double)); cudaMalloc((void**)&d_M2, rowSize * colSize * sizeof(double)); // 将CPU数据拷贝到GPU上 cudaMemcpy(d_M1, h_M1, rowSize * colSize * sizeof(double), cudaMemcpyHostToDevice); // 4.定义kernel的执行配置 int blockSize; int minGridSize; int gridSize; // 获取GPU的信息，并配置最优参数 cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize, mod_kernel, 0, rowSize * colSize); gridSize = (rowSize * colSize + blockSize - 1) / blockSize; // 5.执行核函数 // 取模 mod_kernel <<< gridSize, blockSize >>> (d_M1, d_M2, q, rowSize*colSize); // 6.将GPU数据拷贝到CPU上 cudaMemcpy(h_M2, d_M2, rowSize * colSize * sizeof(double), cudaMemcpyDeviceToHost); // 7.赋值给结果矩阵 for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { M2[i][j] = static_cast<long>(h_M2[i * colSize + j]); } } // 8.清理掉使用过的内存 free(h_M1); free(h_M2); cudaFree(d_M1); cudaFree(d_M2); return M2; } long** uniformMat(long rowSize, long colSize, long minValue, long maxValue) { long** mat = new long* [rowSize]; for (long i = 0; i < rowSize; i++) mat[i] = new long[colSize]; srand((unsigned)time(NULL)); for (long i = 0; i < rowSize; i++) { for (long j = 0; j < colSize; j++) { mat[i][j] = (long)(rand() % (maxValue - minValue + 1)) + minValue; } } return mat; } int main(void){ // 1.定义cuda消息处理器 cublasHandle_t cuHandle; cublasStatus_t status = cublasCreate(&cuHandle); // 2.定义要取模的矩阵 long q = 3; // 取模数 long rowSize = 4; // 矩阵行数 long colSize = 6; // 矩阵列数 long** M1 = uniformMat(rowSize, colSize, 5, 9); // 生成一个rowSize*colSize的随机矩阵，元素值在[5, 9]内随机选取 // 3.进行矩阵元素取模运算 long** M2 = mat_mod_cu(M1, rowSize, colSize, q, cuHandle); // 4.输出原矩阵M1和取模后的矩阵M2 cout << "M1: " << endl; for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { cout << M1[i][j] << " "; } cout << endl; } cout << "M2: " << endl; for (int i = 0; i < rowSize; i++) { for (int j = 0; j < colSize; j++) { cout << M2[i][j] << " "; } cout << endl; } // 5.销毁cuda消息处理器 cublasDestroy(cuHandle); return 0; }
在命令行中执行：
nvcc -lcublas mod_test.cu -o mt ./mt 1 2 运算结果： M1: 9 6 5 9 9 8 9 8 6 5 7 9 5 6 6 6 8 7 5 6 5 9 9 5 M2: 0 0 2 0 0 2 0 2 0 2 1 0 2 0 0 0 2 1 2 0 2 0 0 2
成功！
