为何矩阵乘法代码中B矩阵必须预先转置？

问题背景 #

实验要求将以下CPU矩阵乘法代码转换为CUDA程序，实现GPU矩阵乘法并进行耗时测试：

void Multiply(float* A, float* B, float* C, int N)
{
    // B is already transposed
    int loc_a, loc_b, loc_c;
    printf("In Multiply\n");
    for (int i =0; i< N; i++) {
        for (int j = 0; j<N;j++) {
            loc_c = i*N+j;
            loc_a = i*N;
            loc_b = j*N;
            C[loc_c] = 0.0f;
            for (int k=0;k<N;k++) {
                C[loc_c] += A[loc_a]*B[loc_b];
                loc_a++;
                loc_b++;
            }
        }
    }
}

我对这段代码存在疑惑：为何进行矩阵乘法时，B矩阵必须预先转置？

核心原因：优化内存访问连续性 #

不管是CPU还是GPU，硬件缓存对连续内存访问的优化力度极大——连续读取的数据会被缓存起来，后续访问直接从高速缓存取，比从主存快几个数量级。

看原代码的计算逻辑：
计算C[i][j]时，A矩阵取的是第i行的连续元素（loc_a从i*N开始逐次+1），属于连续内存访问，能充分利用缓存。如果B矩阵不转置，原本计算时需要取B的第j列元素，也就是B[j + k*N]，这种访问是跨步式的（步长为N），内存地址不连续，缓存根本存不下这些分散的数据，每次读取都要从主存加载，速度慢得离谱。

把B转置之后，原来的第j列变成了第j行，此时取B的元素就是B[j*N + k]，和A一样是连续内存访问，缓存命中率拉满，整体乘法的执行效率会显著提升。

逻辑正确性验证 #

原代码的乘法逻辑等价于标准矩阵乘法C = A × B_original，其中B是B_original的转置。因为标准矩阵乘法中C[i][j] = Σ(A[i][k] × B_original[k][j])，转置后B[j][k] = B_original[k][j]，所以代码里的Σ(A[loc_a] × B[loc_b])就是标准乘法的计算逻辑，结果完全正确。

完整CPU版本源代码 #

#ifndef CUDA_LAB_MATRIX_H
#define CUDA_LAB_MATRIX_H

#include <stdio.h>
#include <stdlib.h>
#include <chrono>
#include <iostream>
#include <cmath>

#include <cstdarg> // for va_list, va_start, va_arg, va_end

void CreateMatrix(float* &A, int count, ...) {
    A = (float*)malloc(count * sizeof(float));
    va_list args;
    va_start(args, count); // start the variable argument list

    for (int i = 0; i < count; i++) {
        A[i] = va_arg(args, int); // get the next argument from the list
    }

    va_end(args); // end the variable argument list
}

void Transpose(float *A, float* At, int N){
    for (int i = 0; i<N;i++ ) {
        int loc = i*N;
        int loc_t;
        for (int j=0;j<N;j++) {
            loc_t = j*N+i;
            At[loc_t] = A[loc];
            loc++;
        }
    }
}

void Multiply(float* A, float* B, float* C, int N)
{
    C = (float*)malloc(N * sizeof(float));
    // B is transposed
    int loc_a, loc_b, loc_c;
    printf("In Multiply\n");
    for (int i =0; i< N; i++) {
        for (int j = 0; j<N;j++) {
            loc_c = i*N+j;
            loc_a = i*N;
            loc_b = j*N;
            C[loc_c] = 0.0f;
            for (int k=0;k<N;k++) {
                float temp = A[loc_a]*B[loc_b];
                C[loc_c] += temp;
                loc_a++;
                loc_b++;
            }
        }
    }
}

void PrintMat(float* A, int row, int ext_row, int col, int ext_col, int N)
{
    int cur_row;
    int loc;
    cur_row = row;
    for (int i = 0; i< ext_row; i++)
    {
        loc = cur_row*N +col;
        for (int j=0; j< ext_col;j++)
        {
            printf("%f  ",A[loc+j]);
        }
        printf("\n");
        cur_row++;
    }
}

void CompareMatrices(float* A, float *B, int N)
{
    int count =0;
    float Sum=0.0f;
    float diff;
    int loc =0;

    for (int i =0;i<N;i++)
    {
        for (int j=0;j<N;j++)
        {
            if (A[loc]!=B[loc])
            {
                Sum += fabs(A[loc]-B[loc]);
                count++;
            }
            loc++;
        }
    }
    printf("Difference: %f\n",Sum);
    printf("Count: %d\n",count);
}

#endif //CUDA_LAB_MATRIX_H

CUDA实现参考思路（满足实验要求） #

要转成CUDA程序，核心是延续内存连续访问的优化思路，同时利用GPU的并行能力：

• 用线程网格分配任务：每个线程负责计算C中的一个元素C[i][j]，通过blockIdx和threadIdx计算对应的i和j
• 预先转置B矩阵（可以在CPU上完成后拷贝到GPU，也可以在GPU上执行转置核函数）
• 使用CUDA事件API（cudaEventCreate、cudaEventRecord、cudaEventElapsedTime）测试GPU乘法的耗时
• 严格遵循CUDA内存流程：cudaMalloc分配GPU内存 → cudaMemcpy把CPU数据拷贝到GPU → 启动核函数计算 → cudaMemcpy把结果拷贝回CPU → cudaFree释放GPU内存

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