绥化市网站建设_网站建设公司_Ruby_seo优化

假设我们要计算C=A×BC = A \times BC=A×B，其中AAA是M×KM \times KM×K矩阵，BBB是K×NK \times NK×N矩阵，CCC是M×NM \times NM×N矩阵。

1. 矩阵乘法回顾

矩阵CCC中任意元素Ci,jC_{i, j}Ci,j​的值，是通过将矩阵AAA的第iii行与矩阵BBB的第jjj列进行点积得到的：

Ci,j=∑k=0K−1Ai,k⋅Bk,jC_{i, j} = \sum_{k=0}^{K-1} A_{i, k} \cdot B_{k, j}Ci,j​=k=0∑K−1​Ai,k​⋅Bk,j​

2. 内存布局：行主序（Row-Major）

在 C/C++ 中，多维数组通常以**行主序（Row-Major Order）**存储在内存中，即一行接一行地存储。

• 对于M×NM \times NM×N矩阵XXX，元素Xi,jX_{i, j}Xi,j​的一维索引是：

  Index(i,j)=i×N+j\text{Index}(i, j) = i \times N + jIndex(i,j)=i×N+j

  其中NNN是矩阵的列数。

3. CUDA 编程步骤

步骤 1: 定义 Kernel 启动配置

我们需要启动M×NM \times NM×N个线程，让每个线程负责计算输出矩阵CCC中的一个元素Ci,jC_{i, j}Ci,j​。

步骤 2: 计算线程的全局索引(i,j)(i, j)(i,j)

Kernel 内部的关键是，将线程的xxx和yyy索引转换为矩阵CCC的行索引iii和列索引jjj。

• 行索引iii：对应于线程的全局yyy坐标。

• 列索引jjj：对应于线程的全局xxx坐标。

i=blockIdx.y×blockDim.y+threadIdx.yj=blockIdx.x×blockDim.x+threadIdx.xi = \text{blockIdx.y} \times \text{blockDim.y} + \text{threadIdx.y} \\ j = \text{blockIdx.x} \times \text{blockDim.x} + \text{threadIdx.x}i=blockIdx.y×blockDim.y+threadIdx.yj=blockIdx.x×blockDim.x+threadIdx.x

步骤 3: 实现 Kernel 逻辑

如果iii在MMM范围内，jjj在NNN范围内，线程执行AAA的第iii行和BBB的第jjj列的点积计算。

4. 完整的 CUDA C++ 代码实现

#include <stdio.h> #include <cuda_runtime.h> // 矩阵维度定义 (M x K) * (K x N) = (M x N) #define M 1024 #define K 512 #define N 1024 // 矩阵乘法 Kernel：每个线程计算 C 的一个元素 __global__ void matrixMul(const float* A, const float* B, float* C, int M, int K, int N) { // 1. 计算线程在输出矩阵 C 中的全局坐标 (i, j) // i 是行索引 (对应 Block Y 和 Thread Y) int i = blockIdx.y * blockDim.y + threadIdx.y; // j 是列索引 (对应 Block X 和 Thread X) int j = blockIdx.x * blockDim.x + threadIdx.x; // 2. 边界检查 (确保线程不超过矩阵 C 的维度 M x N) if (i < M && j < N) { float Cij = 0; // 3. 执行点积计算 (A 的第 i 行 和 B 的第 j 列) for (int k = 0; k < K; ++k) { // A[i, k] 的一维索引: i * K + k float Aik = A[i * K + k]; // B[k, j] 的一维索引: k * N + j float Bkj = B[k * N + j]; Cij += Aik * Bkj; } // 4. 将结果写回输出矩阵 C // C[i, j] 的一维索引: i * N + j C[i * N + j] = Cij; } } void setupAndLaunch() { // --- 1. 内存分配与初始化 --- size_t size_A = (size_t)M * K * sizeof(float); size_t size_B = (size_t)K * N * sizeof(float); size_t size_C = (size_t)M * N * sizeof(float); // Host 端内存 float *h_A, *h_B, *h_C; h_A = (float*)malloc(size_A); h_B = (float*)malloc(size_B); h_C = (float*)malloc(size_C); // 初始化 A 和 B (略，假设已初始化) // for (int i = 0; i < M*K; ++i) h_A[i] = ...; // Device 端内存 float *d_A, *d_B, *d_C; cudaMalloc((void**)&d_A, size_A); cudaMalloc((void**)&d_B, size_B); cudaMalloc((void**)&d_C, size_C); // --- 2. 数据传输：Host -> Device --- cudaMemcpy(d_A, h_A, size_A, cudaMemcpyHostToDevice); cudaMemcpy(d_B, h_B, size_B, cudaMemcpyHostToDevice); // --- 3. 配置启动参数与 Kernel 启动 --- // 定义 Thread Block 大小 (例如 16x16) int TILE_SIZE = 16; dim3 threadsPerBlock(TILE_SIZE, TILE_SIZE); // 计算 Grid 大小 // M / TILE_SIZE 向上取整，作为 Grid 的 Y 维度 int gridX = (N + TILE_SIZE - 1) / TILE_SIZE; // N / TILE_SIZE 向上取整，作为 Grid 的 X 维度 int gridY = (M + TILE_SIZE - 1) / TILE_SIZE; dim3 blocksPerGrid(gridX, gridY); // 启动 Kernel matrixMul<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, M, K, N); // 检查是否有 Kernel 启动错误 cudaError_t err = cudaGetLastError(); if (err != cudaSuccess) { fprintf(stderr, "CUDA launch error: %s\n", cudaGetErrorString(err)); return; } // 等待 GPU 完成 cudaDeviceSynchronize(); // --- 4. 数据传输：Device -> Host --- cudaMemcpy(h_C, d_C, size_C, cudaMemcpyDeviceToHost); // ... 验证结果 h_C ... // --- 5. 清理 --- cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); free(h_B); free(h_C); } int main() { setupAndLaunch(); return 0; }

5. 性能分析（该简单实现的局限性）

这个“从零开始”的实现虽然功能正确，但性能非常低效。

1. 全局内存带宽限制：计算Ci,jC_{i, j}Ci,j​时，每个线程都需要读取KKK个AAA元素和KKK个BBB元素。如果KKK很大，会导致大量的全局内存读取。

2. 内存访问不连续：

   • 读取Ai,kA_{i, k}Ai,k​（同一行）是连续的，对 GPU 有利。

   • 读取Bk,jB_{k, j}Bk,j​（同一列）在行主序存储中是不连续的（跨步访问）。这会导致严重的内存合并（Memory Coalescing）失败，极大降低内存吞吐量。

优化方向：使用共享内存（Shared Memory）

要实现高性能的矩阵乘法，必须采用 **Tile Blocking（分块）**技术，并将输入矩阵的块加载到 **共享内存（Shared Memory）**中，以实现数据重用和内存合并。

这是下一个层次的优化，涉及更复杂的 Kernel 逻辑和资源管理。