VS2015 cuda c++ matrixMul矩阵乘法经典案例详解

文章目录

• • 1.源码演示
  • 2.头文件和宏定义部分（1-20行）
  • 3.CUDA Kernel函数（24-90行）
  • • 索引计算（30-37行）
    • 矩阵A的遍历参数（40-46行）
    • 矩阵B的遍历参数（49-52行）
    • 核心计算循环（59-86行）
    • 结果写回（89-90行）
  • 4.主机端辅助函数（92-98行）
  • 5.主测试函数 MatrixMultiply（101-271行）
  • • 主机内存分配（104-115行）
    • GPU内存分配（126-130行）
    • 性能计时设置（135-138行）
    • 异步内存传输（142-143行）
    • 内核启动配置（146-148行）
    • 预热和性能测试（155-180行）
    • 性能计算（189-200行）
    • 结果验证（211-228行）
  • 6.main函数（273-321行）
  • • 命令行参数解析（278-287行）
    • 默认配置（299-302行）
    • 维度验证（312-316行）
  • 7.关键优化技术总结

1.源码演示

/**
* Matrix multiplication: C = A * B.
* Host code.
*
* This sample implements matrix multiplication which makes use of shared memory
* to ensure data reuse, the matrix multiplication is done using tiling approach.
* It has been written for clarity of exposition to illustrate various CUDA programming
* principles, not with the goal of providing the most performant generic kernel for matrix multiplication.
* See also:
* V. Volkov and J. Demmel, "Benchmarking GPUs to tune dense linear algebra,"
* in Proc. 2008 ACM/IEEE Conf. on Supercomputing (SC '08),
* Piscataway, NJ: IEEE Press, 2008, pp. Art. 31:1-11.
*/

// System includes
#include <stdio.h>
#include <assert.h>

// CUDA runtime
#include <cuda_runtime.h>

// Helper functions and utilities to work with CUDA
#include <helper_functions.h>
#include <helper_cuda.h>

/**
* Matrix multiplication (CUDA Kernel) on the device: C = A * B
* wA is A's width and wB is B's width
*/
template <int BLOCK_SIZE> __global__ void MatrixMulCUDA(float *C, float *A,
float *B, int wA,
int wB) {
// Block index
int bx = blockIdx.x;
int by = blockIdx.y;

// Thread index
int tx = threadIdx.x;
int ty = threadIdx.y;

// Index of the first sub-matrix of A processed by the block
int aBegin = wA * BLOCK_SIZE * by;

// Index of the last sub-matrix of A processed by the block
int aEnd = aBegin + wA – 1;

// Step size used to iterate through the sub-matrices of A
int aStep = BLOCK_SIZE;

// Index of the first sub-matrix of B processed by the block
int bBegin = BLOCK_SIZE * bx;

// Step size used to iterate through the sub-matrices of B
int bStep = BLOCK_SIZE * wB;

// Csub is used to store the element of the block sub-matrix
// that is computed by the thread
float Csub = 0;

// Loop over all the sub-matrices of A and B
// required to compute the block sub-matrix
for (int a = aBegin, b = bBegin;
a <= aEnd;
a += aStep, b += bStep) {
// Declaration of the shared memory array As used to
// store the sub-matrix of A
__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];

// Declaration of the shared memory array Bs used to
// store the sub-matrix of B
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

// Load the matrices from device memory
// to shared memory; each thread loads
// one element of each matrix
As[ty][tx] = A[a + wA * ty + tx];
Bs[ty][tx] = B[b + wB * ty + tx];

// Synchronize to make sure the matrices are loaded
__syncthreads();

// Multiply the two matrices together;
// each thread computes one element
// of the block sub-matrix
#pragma unroll

for (int k = 0; k < BLOCK_SIZE; ++k) {
Csub += As[ty][k] * Bs[k][tx];
}

// Synchronize to make sure that the preceding
// computation is done before loading two new
// sub-matrices of A and B in the next iteration
__syncthreads();
}

// Write the block sub-matrix to device memory;
// each thread writes one element
int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
C[c + wB * ty + tx] = Csub;
}

void ConstantInit(float *data, int size, float val) {
for (int i = 0; i < size; ++i) {
data[i] = val;
}
}

/**
* Run a simple test of matrix multiplication using CUDA
*/
int MatrixMultiply(int argc, char **argv,
int block_size, const dim3 &dimsA,
const dim3 &dimsB) {
// Allocate host memory for matrices A and B
unsigned int size_A = dimsA.x * dimsA.y;
unsigned int mem_size_A = sizeof(float) * size_A;
float *h_A = reinterpret_cast<float *>(malloc(mem_size_A));
unsigned int size_B = dimsB.x * dimsB.y;
unsigned int mem_size_B = sizeof(float) * size_B;
float *h_B = reinterpret_cast<float *>(malloc(mem_size_B));
cudaStream_t stream;

// Initialize host memory
const float valB = 0.01f;
ConstantInit(h_A, size_A, 1.0f);
ConstantInit(h_B, size_B, valB);

// Allocate device memory
float *d_A, *d_B, *d_C;

// Allocate host matrix C
dim3 dimsC(dimsB.x, dimsA.y, 1);
unsigned int mem_size_C = dimsC.x * dimsC.y * sizeof(float);
float *h_C = reinterpret_cast<float *>(malloc(mem_size_C));

if (h_C == NULL) {
fprintf(stderr, "Failed to allocate host matrix C!\\n");
exit(EXIT_FAILURE);
}

checkCudaErrors(cudaMalloc(reinterpret_cast<void **>(&d_A), mem_size_A));
checkCudaErrors(cudaMalloc(reinterpret_cast<void **>(&d_B), mem_size_B));
checkCudaErrors(cudaMalloc(reinterpret_cast<void **>(&d_C), mem_size_C));
// Allocate CUDA events that we'll use for timing
cudaEvent_t start, stop;
checkCudaErrors(cudaEventCreate(&start));
checkCudaErrors(cudaEventCreate(&stop));

checkCudaErrors(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));

// copy host memory to device
checkCudaErrors(cudaMemcpyAsync(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice, stream));
checkCudaErrors(cudaMemcpyAsync(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice, stream));

// Setup execution parameters
dim3 threads(block_size, block_size);
dim3 grid(dimsB.x / threads.x, dimsA.y / threads.y);

// Create and start timer
printf("Computing result using CUDA Kernel…\\n");

// Performs warmup operation using matrixMul CUDA kernel
if (block_size == 16) {
MatrixMulCUDA<16> <<< grid, threads, 0, stream>>>(d_C, d_A, d_B,
dimsA.x, dimsB.x);
} else {
MatrixMulCUDA<32> <<< grid, threads, 0, stream>>>(d_C, d_A, d_B,
dimsA.x, dimsB.x);
}

printf("done\\n");
checkCudaErrors(cudaStreamSynchronize(stream));

// Record the start event
checkCudaErrors(cudaEventRecord(start, stream));

// Execute the kernel
int nIter = 300;

for (int j = 0; j < nIter; j++) {
if (block_size == 16) {
MatrixMulCUDA<16> <<<grid, threads, 0, stream>>>(d_C, d_A, d_B,
dimsA.x, dimsB.x);
} else {
MatrixMulCUDA<32> <<<grid, threads, 0, stream>>>(d_C, d_A, d_B,
dimsA.x, dimsB.x);
}
}

// Record the stop event
checkCudaErrors(cudaEventRecord(stop, stream));

// Wait for the stop event to complete
checkCudaErrors(cudaEventSynchronize(stop));

float msecTotal = 0.0f;
checkCudaErrors(cudaEventElapsedTime(&msecTotal, start, stop));

// Compute and print the performance
float msecPerMatrixMul = msecTotal / nIter;
double flopsPerMatrixMul = 2.0 * static_cast<double>(dimsA.x) *
static_cast<double>(dimsA.y) *
static_cast<double>(dimsB.x);
double gigaFlops = (flopsPerMatrixMul * 1.0e-9f) /
(msecPerMatrixMul / 1000.0f);
printf(
"Performance= %.2f GFlop/s, Time= %.3f msec, Size= %.0f Ops," \\
" WorkgroupSize= %u threads/block\\n",
gigaFlops,
msecPerMatrixMul,
flopsPerMatrixMul,
threads.x * threads.y);

// Copy result from device to host
checkCudaErrors(cudaMemcpyAsync(h_C, d_C, mem_size_C, cudaMemcpyDeviceToHost, stream));
checkCudaErrors(cudaStreamSynchronize(stream));

printf("Checking computed result for correctness: ");
bool correct = true;

// test relative error by the formula
// |<x, y>_cpu – <x,y>_gpu|/<|x|, |y|> < eps
double eps = 1.e-6; // machine zero

for (int i = 0; i < static_cast<int>(dimsC.x * dimsC.y); i++) {
double abs_err = fabs(h_C[i] – (dimsA.x * valB));
double dot_length = dimsA.x;
double abs_val = fabs(h_C[i]);
double rel_err = abs_err / abs_val / dot_length;

if (rel_err > eps) {
printf("Error! Matrix[%05d]=%.8f, ref=%.8f error term is > %E\\n",
i, h_C[i], dimsA.x * valB, eps);
correct = false;
}
}

printf("%s\\n", correct ? "Result = PASS" : "Result = FAIL");

// Clean up memory
free(h_A);
free(h_B);
free(h_C);
checkCudaErrors(cudaFree(d_A));
checkCudaErrors(cudaFree(d_B));
checkCudaErrors(cudaFree(d_C));
checkCudaErrors(cudaEventDestroy(start));
checkCudaErrors(cudaEventDestroy(stop));
printf("\\nNOTE: The CUDA Samples are not meant for performance"\\
"measurements. Results may vary when GPU Boost is enabled.\\n");

if (correct) {
return EXIT_SUCCESS;
} else {
return EXIT_FAILURE;
}
}

/**
* Program main
*/
int main(int argc, char **argv) {
printf("[Matrix Multiply Using CUDA] – Starting…\\n");

if (checkCmdLineFlag(argc, (const char **)argv, "help") ||
checkCmdLineFlag(argc, (const char **)argv, "?")) {
printf("Usage -device=n (n >= 0 for deviceID)\\n");
printf(" -wA=WidthA -hA=HeightA (Width x Height of Matrix A)\\n");
printf(" -wB=WidthB -hB=HeightB (Width x Height of Matrix B)\\n");
printf(" Note: Outer matrix dimensions of A & B matrices" \\
" must be equal.\\n");

exit(EXIT_SUCCESS);
}

// This will pick the best possible CUDA capable device, otherwise
// override the device ID based on input provided at the command line
int dev = findCudaDevice(argc, (const char **)argv);

int block_size = 32;

dim3 dimsA(5 * 2 * block_size, 5 * 2 * block_size, 1);
dim3 dimsB(5 * 4 * block_size, 5 * 2 * block_size, 1);

// width of Matrix A
if (checkCmdLineFlag(argc, (const char **)argv, "wA")) {
dimsA.x = getCmdLineArgumentInt(argc, (const char **)argv, "wA");
}

// height of Matrix A
if (checkCmdLineFlag(argc, (const char **)argv, "hA")) {
dimsA.y = getCmdLineArgumentInt(argc, (const char **)argv, "hA");
}

// width of Matrix B
if (checkCmdLineFlag(argc, (const char **)argv, "wB")) {
dimsB.x = getCmdLineArgumentInt(argc, (const char **)argv, "wB");
}

// height of Matrix B
if (checkCmdLineFlag(argc, (const char **)argv, "hB")) {
dimsB.y = getCmdLineArgumentInt(argc, (const char **)argv, "hB");
}

if (dimsA.x != dimsB.y) {
printf("Error: outer matrix dimensions must be equal. (%d != %d)\\n",
dimsA.x, dimsB.y);
exit(EXIT_FAILURE);
}

printf("MatrixA(%d,%d), MatrixB(%d,%d)\\n", dimsA.x, dimsA.y,
dimsB.x, dimsB.y);

int matrix_result = MatrixMultiply(argc, argv, block_size, dimsA, dimsB);

exit(matrix_result);
}

这是一个CUDA C++实现的矩阵乘法示例，使用分块(tiling)和共享内存(shared memory)来优化性能。将逐行解析代码：

2.头文件和宏定义部分（1-20行）

// System includes
#include <stdio.h>
#include <assert.h>

// CUDA runtime
#include <cuda_runtime.h>

// Helper functions and utilities to work with CUDA
#include <helper_functions.h>
#include <helper_cuda.h>

• 包含标准C库和CUDA运行时库
• helper_functions.h和helper_cuda.h是CUDA Samples提供的辅助函数，用于设备选择、错误检查等

3.CUDA Kernel函数（24-90行）

template <int BLOCK_SIZE> __global__ void MatrixMulCUDA(float *C, float *A,
float *B, int wA,
int wB)

• __global__：声明这是一个CUDA kernel函数，在GPU上执行，从CPU调用
• template <int BLOCK_SIZE>：编译时确定块大小（16或32），允许编译器优化

索引计算（30-37行）

int bx = blockIdx.x; // 块在x方向上的索引
int by = blockIdx.y; // 块在y方向上的索引
int tx = threadIdx.x; // 线程在块内的x索引
int ty = threadIdx.y; // 线程在块内的y索引

矩阵A的遍历参数（40-46行）

int aBegin = wA * BLOCK_SIZE * by; // 当前块处理的A的起始行索引
int aEnd = aBegin + wA – 1; // 结束条件
int aStep = BLOCK_SIZE; // 步长（每次移动一个块的大小）

• 矩阵A是M×K大小（M = dimsA.y, K = dimsA.x）
• wA是矩阵A的宽度（即K）
• 每个块负责计算一个BLOCK_SIZE × BLOCK_SIZE的子矩阵

矩阵B的遍历参数（49-52行）

int bBegin = BLOCK_SIZE * bx; // 当前块处理的B的起始列索引
int bStep = BLOCK_SIZE * wB; // 步长（跳过整个子矩阵的宽度）

• 矩阵B是K×N大小（K = dimsB.y, N = dimsB.x）
• wB是矩阵B的宽度（即N）

核心计算循环（59-86行）

float Csub = 0; // 每个线程累加的结果
for (int a = aBegin, b = bBegin; a <= aEnd; a += aStep, b += bStep)

• 外层循环遍历所有需要计算的子矩阵对

共享内存声明（64-70行）

__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

• __shared__：声明共享内存，同一块内的线程共享
• 用于缓存当前需要计算的子矩阵，减少全局内存访问

数据加载（73-74行）

As[ty][tx] = A[a + wA * ty + tx];
Bs[ty][tx] = B[b + wB * ty + tx];

• 每个线程加载一个元素到共享内存
• 地址计算：a + wA * ty + tx是A矩阵中元素的行优先索引

同步（77行）

__syncthreads();

• 块内同步，确保所有线程完成数据加载后再进行计算

计算部分（81-83行）

#pragma unroll // 告诉编译器展开循环
for (int k = 0; k < BLOCK_SIZE; ++k) {
Csub += As[ty][k] * Bs[k][tx];
}

• 每个线程计算一个输出元素
• 从共享内存中读取As和Bs进行计算

再次同步（86行）

__syncthreads();

• 确保所有线程完成当前子矩阵计算后，再加载下一组子矩阵

结果写回（89-90行）

int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
C[c + wB * ty + tx] = Csub;

• 将计算结果写回全局内存

4.主机端辅助函数（92-98行）

void ConstantInit(float *data, int size, float val) {
for (int i = 0; i < size; ++i) {
data[i] = val;
}
}

• 用常数值初始化数组

5.主测试函数 MatrixMultiply（101-271行）

主机内存分配（104-115行）

unsigned int size_A = dimsA.x * dimsA.y;
float *h_A = reinterpret_cast<float *>(malloc(mem_size_A));

• h_A表示host端的矩阵A

GPU内存分配（126-130行）

checkCudaErrors(cudaMalloc(reinterpret_cast<void **>(&d_A), mem_size_A));

• d_A表示device端的矩阵A
• checkCudaErrors宏用于检查CUDA API调用是否成功

性能计时设置（135-138行）

cudaEvent_t start, stop;
checkCudaErrors(cudaEventCreate(&start));
checkCudaErrors(cudaEventCreate(&stop));

• CUDA事件用于高精度计时

异步内存传输（142-143行）

checkCudaErrors(cudaMemcpyAsync(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice, stream));

• 异步传输，与GPU计算可以重叠

内核启动配置（146-148行）

dim3 threads(block_size, block_size); // 每个块的线程数（16×16或32×32）
dim3 grid(dimsB.x / threads.x, dimsA.y / threads.y); // 网格中的块数

预热和性能测试（155-180行）

// 预热运行
MatrixMulCUDA<16> <<< grid, threads, 0, stream>>>(…);

// 正式测试，运行300次取平均
int nIter = 300;
for (int j = 0; j < nIter; j++) {
// 执行kernel
}

性能计算（189-200行）

float msecPerMatrixMul = msecTotal / nIter;
double flopsPerMatrixMul = 2.0 * dimsA.x * dimsA.y * dimsB.x;
double gigaFlops = (flopsPerMatrixMul * 1.0e-9f) / (msecPerMatrixMul / 1000.0f);

• FLOPs计算：矩阵乘法C = A × B需要2×M×K×N次浮点运算
• 除以时间得到GFLOPS性能指标

结果验证（211-228行）

double eps = 1.e-6;
for (int i = 0; i < static_cast<int>(dimsC.x * dimsC.y); i++) {
double rel_err = abs_err / abs_val / dot_length;
if (rel_err > eps) {
correct = false;
}
}

• 计算相对误差，验证结果的正确性

6.main函数（273-321行）

命令行参数解析（278-287行）

if (checkCmdLineFlag(argc, (const char **)argv, "help")) {
// 显示帮助信息
}

• 支持命令行指定矩阵维度和GPU设备

默认配置（299-302行）

int block_size = 32; // 默认块大小32
dim3 dimsA(5 * 2 * block_size, 5 * 2 * block_size, 1); // A: 320×320
dim3 dimsB(5 * 4 * block_size, 5 * 2 * block_size, 1); // B: 640×320

维度验证（312-316行）

if (dimsA.x != dimsB.y) {
printf("Error: outer matrix dimensions must be equal.\\n");
}

• 确保矩阵维度匹配：A的列数必须等于B的行数

7.关键优化技术总结

这个实现展示了CUDA编程的核心概念，是学习GPU并行计算的经典示例。