用CUDA实现数学计算
在《第一个Hello World程序》中，我们用CUDA实现了一个Hello World字符串打印程序，了解了CUDA编程的基本概念。但是CUDA的核心优势在于密集型的并行计算，本章将由浅入深的讲解CUDA在数学计算中的应用。
1. 简单的算术运算
我们从最简单的算术运算开始，如实现两个数的加法。完整需求如下：
用CUDA实现一个add的核函数：传入两个int32的参数，计算并获取两个参数相加的结果。
1.1. 完整代码
add.cu:
#include <stdio.h>
#include <cuda_runtime.h>
// 加法核函数 - 每个线程计算两个整数的和
__global__ void add(int a, int b, int* result) {
*result = a + b;
printf("[Device]: %d + %d = %d, thread %d in block %d\n",
a, b, *result, threadIdx.x, blockIdx.x);
}
int main() {
// 主机端数据准备
int h_a = 10; // 第一个加数
int h_b = 20; // 第二个加数
int h_result = 0; // 用于存储结果的变量
// 设备端指针
int* d_result;
// 分配设备内存
cudaMalloc((void**)&d_result, sizeof(int));
// 打印主机端数据
printf("[Host] before calculate, (a:%d, b:%d, c:%d)\n", h_a, h_b, h_result);
// 启动核函数
// 使用1个线程块，每个线程块1个线程
add<<<1, 1>>>(h_a, h_b, d_result);
// 等待GPU完成计算
cudaDeviceSynchronize();
// 将结果从设备拷贝回主机
cudaMemcpy(&h_result, d_result, sizeof(int), cudaMemcpyDeviceToHost);
// 打印最终结果
printf("[Host] after calculate, (a:%d, b:%d, c:%d)\n", h_a, h_b, h_result);
// 释放设备内存
cudaFree(d_result);
return 0;
}
1.2. 编译和运行
# 编译
nvcc ./add.cu
# 运行
./a.out
[Host] before calculate, (a:10, b:20, c:0)
[Device]: 10 + 20 = 30, thread 0 in block 0
[Host] after calculate, (a:10, b:20, c:30)
1.3. 源码分析
1.3.1. 参数传递机制
按值传递的参数 (a, b)
addKernel<<<1, 1>>>(h_a, h_b, d_result);
- 当调用核函数时，
h_a和h_b的值会被复制到GPU的寄存器或常量内存中。 - 每个线程都会获得这些参数的私有副本。
- 适合传递小型数据（如基本数据类型）。
按指针传递的参数 (result)
int* d_result;
cudaMalloc((void**)&d_result, sizeof(int));
d_result是一个指向设备内存的指针。- 必须在调用核函数前使用
cudaMalloc在设备上分配内存，并在核函数调用结束后调用cudaFree释放内存。 - 核函数通过解引用指针(
*result)来访问和修改设备内存中的值。
1.3.2. 执行配置
addKernel<<<1, 1>>>(h_a, h_b, d_result);
<<<1, 1>>>：执行配置，指定网格和块的维度- 第一个
1：使用1个线程块 - 第二个
1：每个线程块使用1个线程
- 第一个
- 对于这种简单的标量计算，只需要一个线程即可
1.3.3. 同步和数据传输
cudaDeviceSynchronize(); // 等待GPU完成计算
cudaMemcpy(&h_result, d_result, sizeof(int), cudaMemcpyDeviceToHost);
cudaDeviceSynchronize()：阻塞主机代码，直到GPU完成所有操作cudaMemcpy：将数据从设备内存复制回主机内存- 参数4：
cudaMemcpyDeviceToHost指定复制方向（设备→主机）
- 参数4：
2. 一维的向量运算
前面加法的算法只是一个Demo，让你先了解一下用CUDA实现数据计算的流程。因为他只是的一个极简单且单线程的计算，还不如直接在CPU计算效率高。
CUDA的优势在于大量数据的并行计算，现在以一个简单向量加法来了解一下CUDA并行计算的用法。
用CUDA实现一个向量加法运算的核函数：传入两个长度相同的int32的数组(向量)，针对两个数组的每一个元，素计算两数之后，并将结果保存至另一个数组，然后返回结果数组。
2.1. 完整代码
#include <stdio.h>
#include <cuda_runtime.h>
// 向量加法内核 - 每个线程处理一个元素
__global__ void vectorAdd(const int32_t *a, const int32_t *b, int32_t *c, int32_t n)
{
// 计算全局线程索引
int32_t global_id = blockIdx.x * blockDim.x + threadIdx.x;
// 确保不越界
if (global_id < n)
{
c[global_id] = a[global_id] + b[global_id];
}
}
int main()
{
const int32_t n = 1000000; // 100万个元素
int32_t size = n * sizeof(int32_t);
// 主机内存分配
int32_t *h_a = new int32_t[n];
int32_t *h_b = new int32_t[n];
int32_t *h_c = new int32_t[n];
// 初始化数据
for (int32_t i = 0; i < n; i++)
{
h_a[i] = i;
h_b[i] = i * 2;
h_c[i] = 0;
}
// 设备内存分配
int32_t *d_a, *d_b, *d_c;
cudaMalloc(&d_a, size);
cudaMalloc(&d_b, size);
cudaMalloc(&d_c, size);
// 拷贝数据到设备
cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);
// 配置线程结构
int32_t blockSize = 256;
int32_t gridSize = (n + blockSize - 1) / blockSize;
// 启动内核
vectorAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);
// 检查内核启动是否有错误
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("Kernel launch error: %s\n", cudaGetErrorString(err));
return -1;
}
// 等待所有线程完成计算
cudaDeviceSynchronize();
// 拷贝结果回主机
cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);
// 验证结果（只验证前10个数据）
for (int32_t i = 0; i < 10; i++)
{
printf("index_%d: %d + %d = %d\n", i, h_a[i], h_b[i], h_c[i]);
}
// 清理
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
delete[] h_a;
delete[] h_b;
delete[] h_c;
return 0;
}
2.2. 编译和运行
# 编译
nvcc ./vector.cu
# 运行
./a.out
index_0: 0 + 0 = 0
index_1: 1 + 2 = 3
index_2: 2 + 4 = 6
index_3: 3 + 6 = 9
index_4: 4 + 8 = 12
index_5: 5 + 10 = 15
index_6: 6 + 12 = 18
index_7: 7 + 14 = 21
index_8: 8 + 16 = 24
index_9: 9 + 18 = 27
2.3. 源码分析
2.3.1. 并行计算思想
向量加法是一个高度并行的问题：
- 每个元素的计算相互独立
- 不需要线程间的通信或同步
- 非常适合GPU的大规模并行架构
2.3.2. 配置线程结构
int32_t blockSize = 256;
int32_t gridSize = (n + blockSize - 1) / blockSize; // 向上取整
blockSize = 256：每个块有256个线程。块大小通常选择32的倍数，可以更好地利用GPU的warp调度机制，提高效率。gridSize = ceil(n / blockSize)：需要的块数（向上取整）。- 总线程数 =
gridSize × blockSize≥n。确保有足够的线程覆盖所有数据元素，多余的线程会被边界检查排除。
2.3.3. 错误处理
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("Kernel launch error: %s\n", cudaGetErrorString(err));
}
- 检查内核启动错误
- 检查内存分配和拷贝错误
- 使用
cudaDeviceSynchronize()确保内核执行完成
3. 二维的矩阵运算
CUDA的优势在于大量数据的并行计算，现在以一个二维的矩阵运算来了解一下CUDA并行计算的用法。需求如下：
用CUDA实现一个矩阵乘法运算的核函数：传入两个float类型的矩阵A和B，矩阵A的大小为
m*k，矩阵B的大小为k*n，对矩阵A和B进行矩阵乘法运算，生产矩阵C，然后返回矩阵C的运算结果。
3.1. 矩阵乘法原理
矩阵乘法的定义： 矩阵A（m×k）乘以矩阵B（kXn）得到矩阵C（m×n），其中C的第i行第j列的元素为A的第i行与B的第j列对应元素乘积的和。
矩阵乘法 C = A × B，其中：
- A 是 m × k 矩阵
- B 是 k × n 矩阵
- C 是 m × n 矩阵
C 中的每个元素 cᵢⱼ 是 A 的第 i 行和 B 的第 j 列的点积：
cᵢⱼ = Σ(aᵢₜ × bₜⱼ)，其中 t 从 0 到 k-1
3.2. 完整代码
#include <stdio.h>
#include <cuda_runtime.h>
#include <assert.h>
#include <stdlib.h>
#include <time.h>
// (用随机数)初始化矩阵
void initializeMatrix(float *matrix, int rows, int cols)
{
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
matrix[i * cols + j] = static_cast<float>(rand()) / RAND_MAX;
}
}
}
// 简单的矩阵乘法核函数 - 每个线程计算一个输出元素
__global__ void matrixMul(const float *A, const float *B, float *C,
int m, int n, int k)
{
// 计算当前线程对应的输出矩阵中的行和列
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
// 确保不越界
if (row < m && col < k)
{
float sum = 0.0f;
// 计算A的第row行和B的第col列的点积
for (int i = 0; i < n; i++)
{
sum += A[row * n + i] * B[i * k + col];
}
C[row * k + col] = sum;
}
}
void matrixMulFromGPU(const float *A, const float *B, float *C,
int M, int K, int N)
{
// 创建CUDA事件对象用于计时
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
// 记录开始时间
cudaEventRecord(start);
// 计算所需内存大小
size_t size_A = M * K * sizeof(float);
size_t size_B = K * N * sizeof(float);
size_t size_C = M * N * sizeof(float);
// 在设备上分配内存
float *d_A, *d_B, *d_C;
cudaMalloc(&d_A, size_A);
cudaMalloc(&d_B, size_B);
cudaMalloc(&d_C, size_C);
// 将输入数据从主机复制到设备
cudaMemcpy(d_A, A, size_A, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, B, size_B, cudaMemcpyHostToDevice);
// 配置线程结构
// 使用16x16的线程块（256个线程）
dim3 blockDim(16, 16);
// 计算网格大小（向上取整）
dim3 gridDim((N + blockDim.x - 1) / blockDim.x,
(M + blockDim.y - 1) / blockDim.y);
printf("Launching kernel with grid (%d, %d) and block (%d, %d)\n",
gridDim.x, gridDim.y, blockDim.x, blockDim.y);
// 计算矩阵的乘积
matrixMul<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
// 检查内核启动是否有错误
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess)
{
printf("Kernel launch error: %s\n", cudaGetErrorString(err));
// 清理CUDA事件
cudaEventDestroy(start);
cudaEventDestroy(stop);
return;
}
// 等待所有线程完成计算
cudaDeviceSynchronize();
// 将结果从设备复制回主机
cudaMemcpy(C, d_C, size_C, cudaMemcpyDeviceToHost);
// 释放设备内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
// 记录结束时间
cudaEventRecord(stop);
cudaEventSynchronize(stop);
// 计算执行时间
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
printf("GPU execution time: %.2f ms\n", milliseconds);
// 清理CUDA事件
cudaEventDestroy(start);
cudaEventDestroy(stop);
}
void matrixMulFromCPU(const float *A, const float *B, float *C,
int m, int n, int k)
{
// 在CPU上计算参考结果
clock_t cpu_start = clock();
for (int i = 0; i < m; i++)
{
for (int j = 0; j < k; j++)
{
float sum = 0.0f;
for (int t = 0; t < n; t++)
{
sum += A[i * n + t] * B[t * k + j];
}
C[i * k + j] = sum;
}
}
clock_t cpu_end = clock();
double cpu_time = static_cast<double>(cpu_end - cpu_start) / CLOCKS_PER_SEC * 1000.0;
printf("CPU execution time: %.2f ms\n", cpu_time);
}
// 打印矩阵（用于调试）
void printMatrix(const float *matrix, int rows, int cols, const char *name)
{
printf("%s:\n", name);
int R = std::min(5, rows);
int C = std::min(5, cols);
for (int i = 0; i < R; i++)
{
for (int j = 0; j < C; j++)
{
printf("%8.2f ", matrix[i * cols + j]);
}
printf(rows > 5 || cols > 5 ? "...\n" : "\n");
}
if (rows > 5 || cols > 5)
printf("...\n");
}
int main()
{
// 设置随机种子
srand(time(nullptr));
// 设置矩阵维度
// A: m x n, B: n x k, C: m x k
const int M = 512; // A的行数
const int K = 256; // A的列数/B的行数
const int N = 1024; // B的列数/C的列数
printf("Matrix multiplication: A(%dx%d) * B(%dx%d) = C(%dx%d)\n",
M, K, K, N, M, N);
// 在主机上分配和初始化矩阵
float *h_A = new float[M * K];
float *h_B = new float[K * N];
float *h_C = new float[M * N];
float *h_C_ref = new float[M * N]; // 用于参考验证
initializeMatrix(h_A, M, K);
initializeMatrix(h_B, K, N);
// 在GPU中执行计算
matrixMulFromGPU(h_A, h_B, h_C, M, K, N);
// 在CPU中执行计算
matrixMulFromCPU(h_A, h_B, h_C_ref, M, K, N);
// 打印部分结果作为示例
printMatrix(h_A, M, K, "Matrix A (first 5x5)");
printMatrix(h_B, K, N, "Matrix B (first 5x5)");
printMatrix(h_C, M, N, "Matrix C for GPU (first 5x5)");
printMatrix(h_C_ref, M, N, "Matrix C for CPU (first 5x5)");
// 释放主机内存
delete[] h_A;
delete[] h_B;
delete[] h_C;
delete[] h_C_ref;
return 0;
}
3.3. 编译和运行
# 编译
nvcc ./matrix.cu
# 运行
Matrix multiplication: A(512x256) * B(256x1024) = C(512x1024)
Launching kernel with grid (64, 32) and block (16, 16)
GPU execution time: 9.41 ms
CPU execution time: 580.94 ms
Matrix A (first 5x5):
0.38 0.60 0.95 0.21 0.00 ...
0.89 0.51 0.85 0.04 0.24 ...
0.81 0.19 0.33 0.66 0.48 ...
0.73 0.98 0.12 0.36 0.93 ...
0.26 0.40 0.66 0.98 0.74 ...
...
Matrix B (first 5x5):
0.32 0.90 0.48 0.26 0.74 ...
0.04 0.94 0.25 0.81 0.26 ...
0.25 0.75 0.01 0.39 0.94 ...
0.05 0.95 0.93 0.18 0.69 ...
0.38 0.81 0.25 0.86 0.20 ...
...
Matrix C for GPU (first 5x5):
61.79 71.57 64.29 65.21 69.29 ...
67.41 74.25 68.38 71.07 68.77 ...
60.84 69.20 62.10 63.84 64.24 ...
57.32 67.56 61.84 60.65 58.03 ...
64.66 75.07 69.87 68.77 67.60 ...
...
Matrix C for CPU (first 5x5):
61.79 71.57 64.29 65.21 69.29 ...
67.41 74.25 68.38 71.07 68.77 ...
60.84 69.20 62.10 63.84 64.24 ...
57.32 67.56 61.84 60.65 58.03 ...
64.66 75.07 69.87 68.77 67.60 ...
...
说明：
- 从打印结果可以看到
GPU版本使用时间是9.41ms，CPU版本使用时间是580.94ms。 - 通过GPU并行计算，实现了矩阵乘法计算的加速。速度提升了约62倍(580.94 / 9.41)。
- 上面的程序实现不是最优的(是最简单易懂的)，还可以做进步的优化，性能会更高。
3.4. 源码分析
3.4.1. 线程组织
对于矩阵乘法，我们使用二维网格和二维线程块：
dim3 blockDim(16, 16); // 16x16=256个线程/块
dim3 gridDim((k + 15) / 16, (m + 15) / 16); // 计算需要的块数
这种组织方式使得：
- 每个线程块处理输出矩阵的一个16×16子块
- 每个线程计算输出矩阵中的一个元素
3.4.2. 初始化矩阵
srand(time(nullptr)); // 第一行：随机数种子初始化
static_cast<float>(rand()) / RAND_MAX; // 第二行：生成[0, 1)范围内的随机浮点数
- 使用当前时间的时间戳(
time(nullptr))作为随机数生成器的种子，确保每次运行都有不同的随机数序列。 - 生成[0, 1)范围内的随机浮点数，用于初始化矩阵元素。
3.4.3. 性能对比
为了对比GPU并行计算的性能，我们写了一个功能相同的CPU版本的矩阵乘法程序，并对矩阵计算进行时间统计。
-
CUDA事件计时： CUDA事件使用GPU内部的高精度计时器，专门用于测量GPU操作的执行时间。
cudaEventCreate(&start); // 创建CUDA事件用于计时
cudaEventRecord(start); // 记录事件发生时间
cudaEventElapsedTime(&milliseconds, start, stop); // 计算两个事件之间的时间间隔（毫秒）
cudaEventDestroy(start); // 销毁CUDA事件对象 -
CPU时钟计计时器： 用于测量程序使用的CPU时间片的时间。注意：不包括时间偏轮转时的等待时间。
clock_t cpu_start = clock(); // 记录程序执行开始时间
clock_t cpu_end = clock(); // 记录程序执行结束时间
// 计算程序执行使用的CPU时间（毫秒）
double cpu_time = static_cast<double>(cpu_end - cpu_start) / CLOCKS_PER_SEC * 1000.0;