请求排查CUDA分块矩阵乘法核函数中的计算错误

开发者阿泽

2026-7-13

问题背景

我先编写了实现矩阵乘法的CPU代码，输入矩阵：

A = B = [[1 2 3 4],
         [5 6 7 8],
         [9 10 11 12],
         [13 14 15 16]]

正确的乘积结果为：

A*B = [[90 100 110 120], 
       [202 228 254 280], 
       [314 356 398 440], 
       [426 484 542 600]]

对应的CPU代码如下：

#include <cstdarg>
#include <fstream>
#include <iomanip>
#include <iostream>

using namespace std;

const int ROWS1 = 4;//1024;
const int COLS1 = 4;//1024;

const int ROWS2 = 4;//1024;
const int COLS2 = 4;//1024;

const int ROWS3 = ROWS1;
const int COLS3 = COLS2;

const int TILE_ROW_SIZE = 2;
const int TILE_COL_SIZE = 2;

#define IDX(tile_size, tile_i, relative_i) (tile_size * tile_i + relative_i)

void MultiplyAsSumOuterProductOfVectors(int *A, int *B, int *C,
int tile_row_size, int tile_col_size,
int cols1, int rows1, int cols2) {
  for (int tile_i = 0; tile_i < tile_col_size; tile_i++) {//x
      for (int tile_j = 0; tile_j < tile_row_size; tile_j++) {//y
          for (int tile_r = 0; tile_r < tile_col_size; tile_r++) {//x
              for (int cell_r = 0; cell_r < cols1; cell_r++) {//x
                  for (int cell_i = 0; cell_i < rows1; cell_i++) {//y
                      for (int cell_j = 0; cell_j < cols2; cell_j++) {//x
                          int r = IDX(TILE_COL_SIZE, tile_r, cell_r);
                          int i = IDX(TILE_ROW_SIZE, tile_i, cell_i);
                          int j = IDX(TILE_COL_SIZE, tile_j, cell_j);
                          C[i * COLS3 + j] += A[i * COLS1 + r] * B[r * COLS2 + j];
                      }
                  }
              }
          }
      }
  }
}

void printMatrix(int *mat, int rr, int cc) {
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      printf("%d ", mat[i * cc + j]);
    }
    printf("\n");
  }
  printf("\n");
}

void allocateMatrix(int *&a, int rows, int cols) {
  a = new int[rows * cols];
}

void freeMatrix(int *a) {
  delete[] a;
}

void initMatrix(int *mat, int rr, int cc) {
  int init = 1;
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      mat[i * cc + j] = init++;
    }
  }
}

void initMatrixZero(int *mat, int rr, int cc) {
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      mat[i * cc + j] = 0;
    }
  }
}

int main() {
  int *A, *B, *C;

  allocateMatrix(A, ROWS1, COLS1);
  initMatrix(A, ROWS1, COLS1);

  allocateMatrix(B, ROWS2, COLS2);
  initMatrix(B, ROWS2, COLS2);

  allocateMatrix(C, ROWS3, COLS3);
  initMatrixZero(C, ROWS3, COLS3);

  MultiplyAsSumOuterProductOfVectors(A, B, C, TILE_ROW_SIZE, TILE_COL_SIZE, COLS1 / TILE_COL_SIZE, ROWS1 / TILE_ROW_SIZE, COLS2 / TILE_COL_SIZE);

  printMatrix(C, ROWS3, COLS3);

  freeMatrix(A);
  freeMatrix(B);
  freeMatrix(C);

  return 0;
}

随后我将其改写为CUDA程序，但输出结果不符合预期：

66    116     86    136
   146    276    198    328
   226    436    310    520
   306    596    422    712

现提供完整CUDA代码，请求协助排查CUDA核函数中的错误：

#include <iostream>
#include <iomanip>

using namespace std;

const int ROWS1 = 4;//1024;
const int COLS1 = 4;//1024;

const int ROWS2 = 4;//1024;
const int COLS2 = 4;//1024;

const int ROWS3 = ROWS1;
const int COLS3 = COLS2;

const int TILE_ROW_SIZE = 2;//32;
const int TILE_COL_SIZE = 2;//32;

#define IDX(tile_size, tile_i, relative_i) (tile_size * tile_i + relative_i)


__global__ void MultiplyAsSumOuterProductOfVectors(int *A, int *B, int *C,
  int tile_row_size, int tile_col_size,
  int cols1, int rows1, int cols2) 
{
  int tile_i = blockIdx.y;
  int tile_j = blockIdx.x;

  int cell_i = threadIdx.y;
  int cell_j = threadIdx.x;

  for (int tile_r = 0; tile_r < tile_col_size; tile_r++)
  {
    int r = IDX(TILE_COL_SIZE, tile_r, cell_j);

    __shared__ int subA[TILE_ROW_SIZE][TILE_COL_SIZE];
    __shared__ int subB[TILE_ROW_SIZE][TILE_COL_SIZE];

    subA[cell_i][cell_j] = A[IDX(cols1, (tile_i * TILE_ROW_SIZE + cell_i), r)];
    subB[cell_i][cell_j] = B[IDX(cols2, r, (tile_j * TILE_COL_SIZE + cell_j))];

    __syncthreads();

    for (int cell_r = 0; cell_r < TILE_ROW_SIZE; cell_r++) 
    {
      int c_i = tile_i * TILE_ROW_SIZE + cell_i;
      int c_j = tile_j * TILE_COL_SIZE + cell_j;
      
      if (c_i < rows1 && c_j < cols2)
      {
        C[IDX(COLS3, c_i, c_j)] += subA[cell_i][cell_r] * subB[cell_r][cell_j];
      }
    }

    __syncthreads();
  }
}

void printMatrix(int *mat, int rr, int cc) {
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      cout << setw(6) << mat[i * cc + j] << " ";
    }
    cout << endl;
  }
  cout << endl;
}

void allocateMatrix(int *&a, int rows, int cols) {
  a = new int[rows * cols];
}

void freeMatrix(int *a) {
  delete[] a;
}

void initMatrix(int *mat, int rr, int cc) {
  int init = 1;
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      mat[i * cc + j] = init++;
    }
  }
}

void initMatrixZero(int *mat, int rr, int cc) {
  for (int i = 0; i < rr; i++) {
    for (int j = 0; j < cc; j++) {
      mat[i * cc + j] = 0;
    }
  }
}

int main() {
  int *A, *B, *C;
  int *d_A, *d_B, *d_C;

  allocateMatrix(A, ROWS1, COLS1);
  initMatrix(A, ROWS1, COLS1);

  allocateMatrix(B, ROWS2, COLS2);
  initMatrix(B, ROWS2, COLS2);

  allocateMatrix(C, ROWS3, COLS3);
  initMatrixZero(C, ROWS3, COLS3);

  // Allocate device memory
  cudaMalloc((void **)&d_A, ROWS1 * COLS1 * sizeof(int));
  cudaMalloc((void **)&d_B, ROWS2 * COLS2 * sizeof(int));
  cudaMalloc((void **)&d_C, ROWS3 * COLS3 * sizeof(int));

  // Copy input matrices from host to device
  cudaMemcpy(d_A, A, ROWS1 * COLS1 * sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_B, B, ROWS2 * COLS2 * sizeof(int), cudaMemcpyHostToDevice);

  // Set grid and block dimensions
  dim3 gridSize(COLS3 / TILE_COL_SIZE, ROWS3 / TILE_ROW_SIZE);
  dim3 blockSize(TILE_COL_SIZE, TILE_ROW_SIZE);

  // Launch the kernel
  MultiplyAsSumOuterProductOfVectors<<<gridSize, blockSize>>>(d_A, d_B, d_C, TILE_ROW_SIZE, TILE_COL_SIZE, COLS1, ROWS1, COLS2);

  // Copy result matrix from device to host
  cudaMemcpy(C, d_C, ROWS3 * COLS3 * sizeof(int), cudaMemcpyDeviceToHost);

  // Print the result matrix
  cout << "Result Matrix:" << endl;
  printMatrix(C, ROWS3, COLS3);

  // Free device memory
  cudaFree(d_A);
  cudaFree(d_B);
  cudaFree(d_C);

  // Free host memory
  freeMatrix(A);
  freeMatrix(B);
  freeMatrix(C);

  return 0;
}

错误分析与修正

1. 共享内存定义位置错误

共享内存subA和subB被定义在tile_r循环内部，每次循环都会重新分配共享内存，破坏线程同步逻辑，还会降低效率。需将共享内存定义在核函数最外层，循环之外。

2. 矩阵索引计算错误

原代码中使用IDX宏计算矩阵索引时逻辑混乱，导致读取的矩阵元素错误：

读取A矩阵时，正确索引应为(tile_i * TILE_ROW_SIZE + cell_i) * cols1 + (tile_r * TILE_COL_SIZE + cell_j)
读取B矩阵时，正确索引应为(tile_r * TILE_ROW_SIZE + cell_i) * cols2 + (tile_j * TILE_COL_SIZE + cell_j)

3. 循环逻辑错误

循环条件tile_r < tile_col_size错误，应改为tile_r < cols1 / tile_col_size，确保遍历完所有需要的tile块（A的列数/B的行数对应的tile数量）
累加过程直接写入全局内存C可能导致线程竞争，改用局部变量sum完成累加后一次性写入，更安全高效

修正后的CUDA核函数

__global__ void MultiplyAsSumOuterProductOfVectors(int *A, int *B, int *C,
  int tile_row_size, int tile_col_size,
  int cols1, int rows1, int cols2) 
{
  int tile_i = blockIdx.y;
  int tile_j = blockIdx.x;

  int cell_i = threadIdx.y;
  int cell_j = threadIdx.x;

  // 共享内存定义在循环外
  __shared__ int subA[TILE_ROW_SIZE][TILE_COL_SIZE];
  __shared__ int subB[TILE_ROW_SIZE][TILE_COL_SIZE];

  // 当前线程负责的C矩阵元素位置
  int c_i = tile_i * TILE_ROW_SIZE + cell_i;
  int c_j = tile_j * TILE_COL_SIZE + cell_j;

  int sum = 0;

  // 遍历所有需要的tile块
  for (int tile_r = 0; tile_r < cols1 / tile_col_size; tile_r++)
  {
    // 加载A的当前tile块到共享内存
    int a_col = tile_r * tile_col_size + cell_j;
    subA[cell_i][cell_j] = A[c_i * cols1 + a_col];

    // 加载B的当前tile块到共享内存
    int b_row = tile_r * tile_row_size + cell_i;
    subB[cell_i][cell_j] = B[b_row * cols2 + c_j];

    // 等待所有线程完成共享内存加载
    __syncthreads();

    // 块内累加计算
    for (int k = 0; k < tile_col_size; k++) 
    {
      sum += subA[cell_i][k] * subB[k][cell_j];
    }

    // 等待所有线程完成当前轮累加，再进行下一轮tile加载
    __syncthreads();
  }

  // 将最终结果写入C矩阵
  if (c_i < rows1 && c_j < cols2)
  {
    C[c_i * COLS3 + c_j] = sum;
  }
}

额外说明

修正后的代码会输出正确的矩阵乘积结果，同时优化了内存访问模式和线程同步逻辑，适合扩展到更大尺寸的矩阵计算。

内容的提问来源于stack exchange，提问作者user366312

请求排查CUDA分块矩阵乘法核函数中的计算错误

错误分析与修正 #

1. 共享内存定义位置错误 #

2. 矩阵索引计算错误 #

3. 循环逻辑错误 #

修正后的CUDA核函数 #