leetGPU解题笔记（2）

1.题面

矩阵乘法

简单
编写一个程序，在GPU上实现两个32位浮点数矩阵的乘法。给定一个维度为$M\times K$的矩阵$A$和一个维度为$K\times N$的矩阵$B$，计算乘积矩阵$C$，其维度为$M\times N$。所有矩阵均以行优先格式存储。

实现要求

• 仅使用原生特性（不允许使用外部库）
• solve函数签名必须保持不变
• 最终结果必须存储在矩阵$C$中

示例1：

输入：
矩阵$A$（$2\times 3$）：
$\left[\begin{array}{ccc}1.0& 2.0& 3.0\\ 4.0& 5.0& 6.0\end{array}\right]$

矩阵$B$（$3\times 2$）：
$\left[\begin{array}{cc}7.0& 8.0\\ 9.0& 10.0\\ 11.0& 12.0\end{array}\right]$

输出：
矩阵$C$（$2\times 2$）：
$\left[\begin{array}{cc}58.0& 64.0\\ 139.0& 154.0\end{array}\right]$

示例2：

输入：
矩阵$A$（$1\times 2$）：
$\left[\begin{array}{cc}2.5& 3.5\end{array}\right]$

矩阵$B$（$2\times 3$）：
$\left[\begin{array}{ccc}1.0& 2.0& 3.0\\ 4.0& 5.0& 6.0\end{array}\right]$

输出：
矩阵$C$（$1\times 3$）：
$\left[\begin{array}{ccc}16.5& 22.5& 28.5\end{array}\right]$

约束条件

• $1\le M,N,K\le 8192$
• 性能测试使用$M=8192$、$N=6144$、$K=4096$的矩阵

2. 给出的代码

// A, B, C are device pointers (i.e. pointers to memory on the GPU)
void solve(const float* A, const float* B, float* C, int M, int N, int K) {
    dim3 threadsPerBlock(16, 16);
    dim3 blocksPerGrid((K + threadsPerBlock.x - 1) / threadsPerBlock.x,
                       (M + threadsPerBlock.y - 1) / threadsPerBlock.y);
    
    matrix_multiplication_kernel<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, M, N, K);
    cudaDeviceSynchronize();
}

这段代码是使用CUDA在GPU上实现矩阵乘法的核心调度部分。以下是对其关键组件的详细解释：

1. 函数参数与目的

void solve(const float* A, const float* B, float* C, int M, int N, int K)

• 参数含义：
  • A、B、C：指向GPU显存的矩阵指针。
  • M、N、K：矩阵维度。A为M×K，B为K×N，C为M×N。
• 目标：计算矩阵乘积 C = A × B。

2. 线程块与网格配置

dim3 threadsPerBlock(16, 16);
dim3 blocksPerGrid((M + threadsPerBlock.x - 1) / threadsPerBlock.x,
                   (K + threadsPerBlock.y - 1) / threadsPerBlock.y);

注意，现在是对一个M*K的大矩形，使用16_16的小正方形进行的划分。

二维线程块 (threadsPerBlock)

• 每个线程块包含 16×16=256 个线程，以二维矩阵形式组织。
• 这是CUDA处理矩阵运算的典型配置，便于并行计算子矩阵。

二维网格 (blocksPerGrid)

• 网格维度：(blocksPerGrid.x, blocksPerGrid.y)。
• 计算方式：
• 结果矩阵是M_N dot N_K ->M_K
  • blocksPerGrid.x = (M + 15) / 16：这是结果的第一维度，和X做对应
  • blocksPerGrid.y = (K + 15) / 16：这是结果的第二维度，和M做对应。
• 作用：将整个矩阵计算任务划分为多个子任务，每个线程块处理一个子矩阵。

3. 核函数调用

matrix_multiplication_kernel<<<blocksPerGrid, threadsPerBlock>>>(A, B, C, M, N, K);

核函数执行配置

• 网格：二维结构，包含 blocksPerGrid.x × blocksPerGrid.y 个线程块。
• 线程块：二维结构，每个块包含 16×16 个线程。

4. 同步操作

cudaDeviceSynchronize();

• 作用：确保核函数执行完成后才返回，避免后续操作访问未完成的结果。
• 必要性：CUDA核函数调用是异步的，必须显式同步以保证数据一致性。

QA环节

其实输入的也还是一维数组，完全是靠着MN和K划分的逻辑矩阵，传入的是一个float指针。

Q1: dim3 是什么意思，举几个简单的例子？

A1: dim3 的定义与作用

dim3 是 CUDA 中用于定义三维线程块/网格结构的内置数据类型，本质是 unsigned int 的三元组 (x, y, z)。
它简化了线程和线程块的组织方式，支持一维、二维或三维的索引系统，常用于：

1. 定义线程块维度（每个线程块包含的线程数）。
2. 定义网格维度（整个网格包含的线程块数）。

示例1：一维线程块与网格

dim3 threadsPerBlock(256);         // 一维线程块，包含256个线程
dim3 blocksPerGrid(10);            // 一维网格，包含10个线程块

• 线程块结构：threadsPerBlock.x = 256，y=1，z=1。
• 网格结构：blocksPerGrid.x = 10，y=1，z=1。
• 总线程数：10 × 256 = 2560 个线程，适用于一维数据（如向量）。

示例2：二维线程块与网格（矩阵计算）

dim3 threadsPerBlock(16, 16);      // 二维线程块，16×16=256个线程
dim3 blocksPerGrid(4, 4);          // 二维网格，4×4=16个线程块

• 线程块结构：x=16，y=16，z=1，线程按二维矩阵排列。
• 网格结构：x=4，y=4，z=1，线程块按二维矩阵排列。
• 总线程数：16 × 16 × 4 × 4 = 4096 个线程，适用于二维数据（如矩阵）。

示例3：三维线程块（图像处理）

dim3 threadsPerBlock(8, 8, 4);     // 三维线程块，8×8×4=256个线程
dim3 blocksPerGrid(10, 10, 1);     // 三维网格，10×10×1=100个线程块

• 线程块结构：x=8，y=8，z=4，适用于处理三维数据（如图像的RGB通道）。
• 网格结构：x=10，y=10，z=1。
• 总线程数：256 × 100 = 25600 个线程，适用于三维数据并行处理。

Q1.1：为什么需要三维结构？

三维结构（x, y, z）能更自然地映射多维数据：

• 一维：处理向量、数组。
• 二维：处理矩阵、图像（行列）。
• 三维：处理立体数据、视频（行列+时间/通道）。

Q1.2：如何计算线程全局索引？

对于二维网格和二维线程块，全局索引计算方式为：

int globalRow = blockIdx.y * blockDim.y + threadIdx.y;  // 全局行索引
int globalCol = blockIdx.x * blockDim.x + threadIdx.x;  // 全局列索引

例如，在矩阵乘法中，每个线程可通过此公式定位到结果矩阵 C 的对应位置。

Q1.3：维度上限是多少？

• 线程块维度：通常 x,y,z 最大为 1024、1024、64（总线程数 ≤1024）。
• 网格维度：通常 x,y,z 最大为 2^31-1、65535、65535。
  （具体限制取决于 GPU 架构，可通过 cudaDeviceProp 查询）。

总结

dim3 是 CUDA 中组织并行计算的基础工具，通过简洁的三维结构支持灵活的线程布局，尤其适合处理多维数据（如矩阵、图像）。合理使用 dim3 能大幅提升 GPU 计算效率。

题目解答

#include "solve.h"
#include <cuda_runtime.h>

__global__ void matrix_multiplication_kernel(const float* A, const float* B, float* C, int M, int N, int K) {
    int x=threadIdx.x+blockDim.x*blockIdx.x;
    int y=threadIdx.y+blockIdx.y*blockDim.y;
    if(x<M&&y<K){
        float sum=0;
        for(int i=0;i<N;i++){
            sum+=A[x*N+i]*B[i*K+y];
        }
        C[x*K+y]=sum;

    
   }

核心思路就是找到需要的那一个结果的位置，然后进行一个for循环。位置翻译正确即可。
注意的是，这里传入的是一维float指针，不支持二维数组直接使用a[x][y],需要迂回一下计算一维坐标a[x*Y+y]。

其实是非常慢的，击败了4.3%的对手。
针对矩阵乘法的优化浩如烟海，我们暂时不急着研究了，这个复杂度依然是n方的。
本次实验主要追求正确性，以入门为主要追求。