计算矩阵的乘法

使用ChatGPT辅助写作。作为一个土狗，尝试使用CUDA计算

A_{1024 \times 1024} \times B_{1024 \times 1024}

A_{1024 \times 1024} = B_{1024 \times 1024} = J_{1024 \times 1024}

其中 $J$ 表示一个 $全1矩阵$

安装CUDA

使用Windows 11的操作系统，Visual Studio作为集成开发环境。安装CUDA的时候要注意VS的版本与CUDA的版本有一致性的要求。例如我使用的Visual Studio Community 2022，其中MVSC的版本是 $19.42.34436$ ，我之前安装的CUDA 12.1就不能用了，从官网重新下载了12.8的版本。

用VS2022新建CUDA 12.8的工程，修改kernel.cu的内容。

初始化两个全1矩阵

 size_t bytes = N * N * sizeof(int);
 int* h_A = (int*)malloc(bytes);
 int* h_B = (int*)malloc(bytes);
 int* h_C = (int*)malloc(bytes);
 std::fill(h_A, h_A + N * N, 1);
 std::fill(h_B, h_B + N * N, 1);
 std::fill(h_C, h_C + N * N, 0);

其中h_C为存放结果的数组。为了方便内存管理，我们把 $1024 \times 1024$ 的矩阵，按行映射为 $1024 \times 1024$ 大小的数组h_A和h_B，分别代表两个矩阵。也就是说我们想要矩阵A的i行j列，就是 $h_A[i*1024 + j]$

获取CUDA的存储空间

有点类似我们使用opengl的时候，要把主机内存的数据拷贝到显卡的内存里，再使用GPU计算，CUDA同样需要如此。

 int* d_A, * d_B, * d_C;

 cudaMalloc(&d_A, bytes);
 cudaMalloc(&d_B, bytes);
 cudaMalloc(&d_C, bytes);

 cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
 cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);

获取GPU的计算资源

 dim3 threadsPerBlock(16, 16);

 dim3 numBlocks(64, 64);

threadBlock是GPU的计算单元，每个threadBlock有一定数量的GPU线程。

每个GPU线程的计算方法

我们有 $1024 \times 1024$ 个元素需要计算，也就需要对应数来你的线程数我们定义有 $64 \times 64$ 个threadBlock，每个block有 $16 \times 16$ 个线程。 $16 \times 16 \times 64 \times 64 = 1024 \times 1024$ 每个线程负责矩阵C的一个元素。

每个threadBlock有

blockIdx
blockDim
threadIdx

我们已经把所有的线程映射到了C的每个元素上，根据矩阵的乘法定义:

C[i][j] = \sum_{i=0}^{n} A[i][k] \times B[k][j]

i = blockIdx.y \times blockDim.y + threadIdx.y

j = blockIdx.x \times blockDim.x + threadIdx.x

代码如下:

__global__ void matrixMultiplication(int *c, const int *a, const int *b)
{
 int row = blockIdx.y * blockDim.y + threadIdx.y;
 int col = blockIdx.x * blockDim.x + threadIdx.x;
 int sum = 0;
    for (int i = 0; i < N; i++)
    {
  sum += a[row * N + i] * b[i * N + col];
    }
 c[row * N + col] = sum;
}

CUDA代码的调用方式如下

matrixMultiplication <<<numBlocks, threadsPerBlock >>> (d_C, d_A, d_B );

完整代码如下


#include <numeric>

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <stdlib.h>

#define N 1024*8

__global__ void matrixMultiplication(int *c, const int *a, const int *b)
{
 int row = blockIdx.y * blockDim.y + threadIdx.y;
 int col = blockIdx.x * blockDim.x + threadIdx.x;
 int sum = 0;
    for (int i = 0; i < N; i++)
    {
  sum += a[row * N + i] * b[i * N + col];
    }
 c[row * N + col] = sum;
}

int genMatrix()
{
 return 0;
}

int main()
{
 size_t bytes = N * N * sizeof(int);
 int* h_A = (int*)malloc(bytes);
 int* h_B = (int*)malloc(bytes);
 int* h_C = (int*)malloc(bytes);
 std::fill(h_A, h_A + N * N, 1);
 std::fill(h_B, h_B + N * N, 1);
 std::fill(h_C, h_C + N * N, 0);

 int* d_A, * d_B, * d_C;

 cudaMalloc(&d_A, bytes);
 cudaMalloc(&d_B, bytes);
 cudaMalloc(&d_C, bytes);

 cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
 cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);

 dim3 threadsPerBlock(16*8, 16*8);

 dim3 numBlocks(64, 64);

 matrixMultiplication <<<numBlocks, threadsPerBlock >>> (d_C, d_A, d_B );

 cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);

 cudaFree(d_A);
 cudaFree(d_B);
 cudaFree(d_C);
 free(h_A);
 free(h_B);
 free(h_C);
 return 0;
}