CUDA矩阵乘法算法分析

CUDA矩阵乘法算法分析

矩阵乘法是科学计算的基本构建块。此外，矩阵乘法的算法模式具有代表性。许多其他算法与矩阵乘法共享类似的优化技术。因此，矩阵乘法是学习并行编程中最重要的例子之一。

CUDA 矩阵乘法的源代码可在 gitlab 上找到。建议使用 git 获取源代码，它允许提取可能提供的任何更新，并提供一种在试验代码时，对所做的任何更改进行版本控制的方法。 可以使用以下命令，获取代码并执行示例：

git clone https://gitlab.com/ecatue/gpu_matrixmul_cuda.git

cd gpu_matrixmul_cuda

make

./matrixmul

下面关于如何在 GPU 上的 CUDA 中，优化矩阵乘法的例子。

矩阵乘法matrixMul问题

给定一个 M x K 矩阵 A 和一个 K x N 矩阵 B，将 A 乘以 B 并将结果存储到 M x N 矩阵 C 中。

现在以matrixMul 为示例，将展示几种在 GPU 上优化矩阵乘法的技术。它们中的大多数是通用的，可以应用于其他应用程序。这些技术是：

1. 分块计算
2. 内存合并
3. 避免内存库冲突
4. 计算优化。
5. 循环展开
6. 预取

这些优化技术的性能如下图所示。注意：这些优化针对矩阵大小为 8800 x 4096 的 NVIDIA 4096 GT GPU 进行了调整，对于其他 GPU和其他矩阵大小，可能不是最佳的。

将从在CPU上运行的简单串行代码开始，然后逐步进行这些优化。

在 CPU 上的简化计算

void main(){

  define A, B, C

  for i = 0 to M do

    for j = 0 to N do

      /* compute element C(i,j) */

      for k = 0 to K do

        C(i,j) <= C(i,j) + A(i,k) * B(k,j)

      end

    end

  end

}

为了简化，图中使用了平方矩阵 （M==N==K）。下图显示了计算元素 C3,11（红色）的内存占用量。这可以看作是一行 A（蓝色）和一列 B（绿色）的内积。

GPU上的简化计算

/* Codes running on CPU */

 

void main(){

 

    define A_cpu, B_cpu, C_cpu in the CPU memory

    define A_gpu, B_gpu, C_gpu in the GPU memory

 

    memcopy A_cpu to A_gpu

    memcopy B_cpu to B_gpu

 

    dim3 dimBlock(16, 16)

    dim3 dimGrid(N/dimBlock.x, M/dimBlock.y)

 

    matrixMul<<<dimGrid, dimBlock>>>(A_gpu,B_gpu,C_gpu,K)

 

    memcopy C_gpu to C_cpu

 

}

/* Codes running on GPU */

 

__global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

 

    temp <= 0

 

    i <= blockIdx.y * blockDim.y + threadIdx.y    // Row i of matrix C

    j <= blockIdx.x * blockDim.x + threadIdx.x    // Column j of matrix C

 

    for k = 0 to K-1 do

        accu <= accu + A_gpu(i,k) * B_gpu(k,j)

    end

 

    C_gpu(i,j) <= accu

 

}

GPU 上的简化计算分配一个线程来计算矩阵 C 的一个元素。每个线程从全局内存中加载一行矩阵 A 和一列矩阵 B做内积，并将结果存储回全局内存中的矩阵 C。该图显示了存储矩阵A、B 和 C 的，全局内存上的一个线程的内存占用量。

在简化计算中，计算量为 2 x M x N x K flop，而全局内存访问量为 2 x M x N x K字节。计算与内存比率约为 1/4（flop/byte）。因此，简化计算是带宽受限的。

通过平铺增加计算与内存比率

为了增加计算与内存比率，可以应用平铺矩阵乘法。一个线程块计算矩阵 C 的一个图块。线程块中的一个线程计算分块的一个元素。该图显示了一个 32 x 32 矩阵，分为四个 16 x 16 图块。为了完成这种计算这，可以创建四个线程块，每个线程块具有 16 x 16 个线程。

GPU 内核在多次迭代中计算 C。在每次迭代中，一个线程块将 A 的一个分块和 B 的一个分块，从全局内存加载到共享内存，执行计算，并将 C 的临时结果存储在寄存器中。完成所有迭代后，线程块将一个图块C存储到全局内存中。例如，一个线程块可以在两次迭代中，计算 C0,0：C0,0 = A0,0 B0,0 + A0,1 B1,0。

在第一次迭代中，线程块将分块 A0,0 和图块 B0,0，从全局内存加载到共享内存中。每个线程执行内积，以产生一个 C 元素。C的这个元素存储在寄存器中，将在下一次迭代中累积。

在第二次迭代中，线程块将图块 A0,1 和图块 B1,0，从全局内存加载到共享内存中。每个线程执行内积，以产生一个 C 元素，该元素与先前的值一起累积。如果这是最后一次迭代，则将寄存器文件中的 C 元素存储回全局内存中。

CPU 代码保持不变。这里只显示 GPU 内核。

/* Codes running on GPU */

 

__global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

 

    __shared__ float A_tile(blockDim.y, blockDim.x)

    __shared__ float B_tile(blockDim.x, blockDim.y)

 

    accu <= 0

 

    /* Accumulate C tile by tile. */

 

    for tileIdx = 0 to (K/blockDim.x - 1) do

 

        /* Load one tile of A and one tile of B into shared mem */

 

        // Row i of matrix A

        i <= blockIdx.y * blockDim.y + threadIdx.y

        // Column j of matrix A

        j <= tileIdx * blockDim.x + threadIdx.x

        // Load A(i,j) to shared mem

        A_tile(threadIdx.y, threadIdx.x) <= A_gpu(i,j)

        // Load B(j,i) to shared mem

        B_tile(threadIdx.x, threadIdx.y) <= B_gpu(j,i) // Global Mem Not coalesced

        // Synchronize before computation

        __sync()

 

        /* Accumulate one tile of C from tiles of A and B in shared mem */

 

        for k = 0 to threadDim.x do

            // Accumulate for matrix C

            accu <= accu + A_tile(threadIdx.y,k) * B_tile(k,threadIdx.x)

        end

        // Synchronize

        __sync()

 

    end

 

    // Row i of matrix C

    i <= blockIdx.y * blockDim.y + threadIdx.y

    // Column j of matrix C

    j <= blockIdx.x * blockDim.x + threadIdx.x

    // Store accumulated value to C(i,j)

    C_gpu(i,j) <= accu

 

}

在平铺实现中，计算量仍然是 2 x M x N x K flop。但是，使用切片大小B，全局内存访问量为2 x M x N x K / B字节。计算与内存比率约为 B/4（flop/字节）。现在可以通过更改图块大小 B 来调整计算与内存比率。

全局内存合并

C/C++ 中的二维数组是行主数组。在上面的平铺实现中，相邻线程已合并到矩阵A，但没有合并到矩阵B。在列主要语言（如 Fortran）中，问题恰恰相反。一个明显的解决方案是，在将矩阵 B 卸载到 GPU 内存之前，由 CPU 转置矩阵 B。

/* Codes running on GPU */

 

__global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

 

    __shared__ float A_tile(blockDim.y, blockDim.x)

    __shared__ float B_tile(blockDim.x, blockDim.y)

 

    accu <= 0

 

    /* Accumulate C tile by tile. */

 

    for tileIdx = 0 to (K/blockDim.x - 1) do

 

        /* Load one tile of A and one tile of B into shared mem */

 

        // Row i of matrix A

        i <= blockIdx.y * blockDim.y + threadIdx.y

        // Column j of matrix A

        j <= tileIdx * blockDim.x + threadIdx.x

        // Load A(i,j) to shared mem

        A_tile(threadIdx.y, threadIdx.x) <= A_gpu(i,j)

        // Load B(i,j) to shared mem

        B_tile(threadIdx.x, threadIdx.y) <= B_gpu(i,j) // Global Mem Coalesced

        // Synchronize before computation

        __sync()

 

        /* Accumulate one tile of C from tiles of A and B in shared mem */

 

        for k = 0 to threadDim.x do

            // Accumulate for matrix C    // Shared Mem Bank conflict

            accu <= accu + A_tile(threadIdx.y,k) * B_tile(threadIdx.x,k)

        end

        // Synchronize

        __sync()

 

    end

 

    // Row i of matrix C

    i <= blockIdx.y * blockDim.y + threadIdx.y

    // Column j of matrix C

    j <= blockIdx.x * blockDim.x + threadIdx.x

    // Store accumulated value to C(i,j)

    C_gpu(i,j) <= accu

 

}

避免共享内存库冲突

/* Codes running on GPU */

 

__global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

 

    __shared__ float A_tile(blockDim.y, blockDim.x)

    __shared__ float B_tile(blockDim.x, blockDim.y)

 

    accu <= 0

 

    /* Accumulate C tile by tile. */

 

    for tileIdx = 0 to (K/blockDim.x - 1) do

 

        /* Load one tile of A and one tile of B into shared mem */

 

        // Row i of matrix A

        i <= blockIdx.y * blockDim.y + threadIdx.y

        // Column j of matrix A

        j <= tileIdx * blockDim.x + threadIdx.x

        // Load A(i,j) to shared mem

        A_tile(threadIdx.y, threadIdx.x) <= A_gpu(i,j)

        // Load B(i,j) to shared mem

        B_tile(threadIdx.y, threadIdx.x) <= B_gpu(i,j) // No Shared Mem Bank conflict

        // Synchronize before computation

        __sync()

 

        /* Accumulate one tile of C from tiles of A and B in shared mem */

 

        for k = 0 to threadDim.x do

            // Accumulate for matrix C  // No Shared Mem Bank conflict

            accu <= accu + A_tile(threadIdx.y,k) * B_tile(k,threadIdx.x)

        end

        // Synchronize

        __sync()

 

    end

 

    // Row i of matrix C

    i <= blockIdx.y * blockDim.y + threadIdx.y

    // Column j of matrix C

    j <= blockIdx.x * blockDim.x + threadIdx.x

    // Store accumulated value to C(i,j)

    C_gpu(i,j) <= accu

 

}

计算优化

内核是计算绑定的。因此，需要在所有指令中增加有用浮点运算。由于内积消耗大部分时间，因此确保该部分高效非常重要。如果检查内积的二进制代码，会发现 CUDA 中的一行代码在二进制文件中需要两条指令。

/* CUDA code for inner product */

 

accu <= accu + A_tile(threadIdx.y,k) * B_tile(k, threadIdx.x)

/* Disassembled from cubin binary */

 

mov.b32 $r0, s[$ofs4+0x0000]

mad.rn.f32 $r9, s[$ofs1+0x002c], $r0, $r9

流多处理器 （SM） 的当前体系结构，仅允许共享内存中的一个源操作数。但是，计算内积需要来自共享内存的两个源操作数。一种解决方案是将矩阵 A 或矩阵 B 存储到寄存器文件中，但这样的话，寄存器文件中的矩阵就不能被不同的线程共享，从而降低了计算与内存比率。

更好的解决方案是执行外积，而不是内积。在这种情况下，矩阵 A 存储在共享内存中，但矩阵 B 和 C 存储在寄存器中。外积不需要共享矩阵 B 和矩阵 C，因此，每个线程仅在寄存器中存储 B 的一个元素和 C 的分块的一列。外积的计算与内存比与内积相同。

/* CUDA code for outer product */

/* accu[i] and b are stored in register file */

 

accu[i] <= accu[i] + A_tile(i) * b

/* Disassembled from cubin binary */

 

mad.rn.f32 $r9, s[$ofs2+0x0010], $r29, $r9

下面是使用外积，将分块 A0,0 和分块 B0,0 相乘，以便计算 C0,0 的示例。在此示例中，A0,0 为 16 x 16，B0,0 为 16 x 64，C0,0 为 16 x 64。一个包含 64 个线程的线程块用于执行计算 C0,0。

第 1 步：将 A0,0 加载到共享内存。

{{< image src=“/media/GPU_outerproduct_laodA.png” width=“70%”>}}

第 2 步：使用 16 次迭代更新 C0,0。

每个线程在其寄存器中，存储一个 B0,0 元素。每个线程还在其寄存器中，存储一列 C0,0。

迭代 1：A0,0 第一列和 B0,0 第一行之间的外积，并更新 C0,0。

迭代 2：A0,0 第二列和 B0,0 第二行之间的外积，并更新 C0,0。

以类似的方式，继续迭代 3、4、...、15。 迭代 16：A16,0 的第 0 列和 B16,0 的第 0 行之间的外积，并更新 C0,0。

第 3 步：每个线程将一列 C0,0，从其寄存器存储到全局内存。

 循环展开

使用编译附注pragma，使得编译器展开循环。默认情况下，nvcc 将展开内部循环。但除非有编译指示pragma，否则，它不会展开外循环。

#pragma unroll

循环展开有时会对寄存器的使用产生副作用，这可能会限制并发线程的数量。但是，循环展开不会增加 matrixMul 示例中的寄存器使用量。

预取

/* Codes running on GPU */

 

__global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

 

    __shared__ float A_tile0(blockDim.y, blockDim.x)

    __shared__ float A_tile1(blockDim.x, blockDim.y)

 

 

    float *pointer0 = A_tile0

    float *pointer1 = A_tile1

 

    fetch one tile of matrix A_gpu to pointer0

 

    __sync()

 

    /* Accumulate C tile by tile. */

 

    for tileIdx = 0 to (K/blockDim.x - 1) do

 

        prefetch one tile of matrix A_gpu to pointer1

 

        accumulate C using pointer0

 

        __sync()

 

        swap pointer0 and pointer1

 

    end

 

    store tile C to global memory

 

}