向量加法与点积的CUDA并行实现探索

向量加法与点积的CUDA并行实现探索

引言

在高性能计算领域，GPU（图形处理单元）因其强大的并行计算能力而备受青睐。CUDA（Compute Unified Device Architecture）作为NVIDIA提供的并行计算平台和编程模型，使得开发者能够充分利用GPU的并行性来加速计算密集型任务。本文将聚焦于两个基础但重要的线性代数操作——向量加法和点积，探讨它们在CUDA环境下的并行实现方法。

向量加法的CUDA实现

向量加法基础

向量加法是线性代数中最基本的操作之一，它涉及两个相同长度的向量，逐元素相加得到结果向量。例如，给定两个向量A和B，其和向量C的每个元素C[i] = A[i] + B[i]。

CUDA并行化思路

在CUDA中，向量加法的并行化相对直观。每个线程可以负责计算结果向量中的一个元素。假设向量长度为N，我们可以启动N个线程，每个线程根据其全局索引i来计算C[i] = A[i] + B[i]。

代码实现示例

#include <stdio.h>
#include <cuda_runtime.h>

#define N 1000000 // 向量长度

// CUDA核函数，执行向量加法
__global__ void vectorAdd(float *A, float *B, float *C) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < N) {
C[i] = A[i] + B[i];
}
}

int main() {
float *h_A, *h_B, *h_C; // 主机端向量
float *d_A, *d_B, *d_C; // 设备端向量

// 分配主机内存
h_A = (float*)malloc(N * sizeof(float));
h_B = (float*)malloc(N * sizeof(float));
h_C = (float*)malloc(N * sizeof(float));

// 初始化向量
for (int i = 0; i < N; i++) {
h_A[i] = 1.0f;
h_B[i] = 2.0f;
}

// 分配设备内存
cudaMalloc((void**)&d_A, N * sizeof(float));
cudaMalloc((void**)&d_B, N * sizeof(float));
cudaMalloc((void**)&d_C, N * sizeof(float));

// 将数据从主机复制到设备
cudaMemcpy(d_A, h_A, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, N * sizeof(float), cudaMemcpyHostToDevice);

// 定义线程块和网格大小
int blockSize = 256;
int gridSize = (N + blockSize – 1) / blockSize;

// 调用核函数
vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C);

// 将结果从设备复制回主机
cudaMemcpy(h_C, d_C, N * sizeof(float), cudaMemcpyDeviceToHost);

// 验证结果（可选）
for (int i = 0; i < N; i++) {
if (h_C[i] != 3.0f) {
printf("Error at index %d\\n", i);
break;
}
}

// 释放内存
free(h_A);
free(h_B);
free(h_C);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);

return 0;
}

点积的CUDA实现

点积基础

点积（或内积）是两个向量之间的一种运算，结果是一个标量。对于向量A和B，它们的点积定义为A·B = Σ(A[i] * B[i])，其中i从0到N-1。

CUDA并行化思路

点积的并行化稍微复杂一些，因为它涉及到一个求和操作。一种常见的方法是使用并行归约（parallel reduction）技术。每个线程块计算一部分元素的乘积之和，然后将这些部分和进一步归约以得到最终结果。

代码实现示例

#include <stdio.h>
#include <cuda_runtime.h>

#define N 1000000 // 向量长度
#define BLOCK_SIZE 256 // 线程块大小

// CUDA核函数，计算部分和
__global__ void dotProductPart(float *A, float *B, float *partialSums) {
__shared__ float sharedSum[BLOCK_SIZE];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x + threadIdx.x;

float sum = 0.0f;
if (i < N) {
sum = A[i] * B[i];
}

sharedSum[tid] = sum;
__syncthreads();

// 并行归约
for (int s = blockDim.x / 2; s > 0; s >>= 1) {
if (tid < s) {
sharedSum[tid] += sharedSum[tid + s];
}
__syncthreads();
}

if (tid == 0) {
partialSums[blockIdx.x] = sharedSum[0];
}
}

int main() {
float *h_A, *h_B; // 主机端向量
float *d_A, *d_B; // 设备端向量
float *d_partialSums, *h_partialSums; // 部分和
float result = 0.0f; // 最终点积结果

// 分配主机内存
h_A = (float*)malloc(N * sizeof(float));
h_B = (float*)malloc(N * sizeof(float));
h_partialSums = (float*)malloc((N + BLOCK_SIZE – 1) / BLOCK_SIZE * sizeof(float));

// 初始化向量
for (int i = 0; i < N; i++) {
h_A[i] = 1.0f;
h_B[i] = 2.0f;
}

// 分配设备内存
cudaMalloc((void**)&d_A, N * sizeof(float));
cudaMalloc((void**)&d_B, N * sizeof(float));
cudaMalloc((void**)&d_partialSums, (N + BLOCK_SIZE – 1) / BLOCK_SIZE * sizeof(float));

// 将数据从主机复制到设备
cudaMemcpy(d_A, h_A, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, N * sizeof(float), cudaMemcpyHostToDevice);

// 定义网格大小
int gridSize = (N + BLOCK_SIZE – 1) / BLOCK_SIZE;

// 调用核函数计算部分和
dotProductPart<<<gridSize, BLOCK_SIZE>>>(d_A, d_B, d_partialSums);

// 将部分和从设备复制回主机
cudaMemcpy(h_partialSums, d_partialSums, gridSize * sizeof(float), cudaMemcpyDeviceToHost);

// 在主机上完成最后的求和
for (int i = 0; i < gridSize; i++) {
result += h_partialSums[i];
}

printf("Dot product result: %f\\n", result);

// 释放内存
free(h_A);
free(h_B);
free(h_partialSums);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_partialSums);

return 0;
}