cuda入门：如何进行矩阵乘法优化-技术开发专区

cuda入门：如何进行矩阵乘法优化

作者：GXW 编辑：王玉圆 2011-07-06 10:47 来源：GXW的�

　　进一步的改良

　　上一节的结论显示出，矩阵乘法的程序，效率是受限于内存带宽的。那有没有办法降低内存的存取次数呢?答案当然是有的，不然就不会有这一节了 :)

　　要进一步降低内存带宽的使用，可以注意到，在上一节的方法中，虽然 A 矩阵的存取次数被减至最低，但是 B 矩阵的存取次数并没有减少。这是因为我们只将 A 矩阵的 row 加载到 shared memory 中，但是 B 矩阵的 column 也是有被重复使用的。理想上应该也可以避免重复加载才对。不过，由于 B 矩阵的 column 使用的时机，和 A 矩阵的 row 是不同的，所以并不能直接这样做。

　　解决方法是 "blocking"。也就是把整个矩阵乘法的动作，切割成很多小矩阵的乘法。例如，要计算 C 矩阵的 (0, 0) ~ (15, 15) 的值，可以把它想成：

　　A(0~15, 0~15) * B(0~15, 0~15) + A(0~15,16~31) * B(16~31, 0~15)

　　+ A(0~15, 32~47) * B(32~47, 0~15) + ...

　　这样一来，我们就可以把两个小矩阵加载到 shared memory，则小矩阵本身的乘法就不需要再存取任何外部的内存了!这样一来，假设小矩阵的大小是 k，则实际上需要的内存存取次数就会变成约 2k2(n/k)3 = 2n3/k。

　　由于目前 CUDA 每个 block 的 thread 数目最多是 512，因此 k = 16 似乎是一个相当理想的数字(共 256 个 threads)。因此，对于一个 n = 1000 的矩阵来说，我们可以把内存存取的量减少到约 500MB，也就是上一节的存取量的 1/8。理论上，这样应该可以让效率提高八倍(假设没有遇到别的瓶颈)。

　　为了方便进行区块的计算，我们让每个 block 有 16x16 个 threads，再建立 (n/16)x(n/16) 个 blocks。把呼叫 kernel 的地方改成：

    int bx = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
    dim3 blocks(bx, bx);
    dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
    matMultCUDA<<<blocks, threads>>>(ac, pitch_a / sizeof(float),
        bc, pitch_b / sizeof(float), cc, pitch_c / sizeof(float), n);

　　BLOCK_SIZE 则是定义成 16。dim3 是 CUDA 的一种数据型态，表示一个 3D 的向量。在这里，我们透过 dim3 来建立 16x16 个 threads 的 block，和 (n/16)x(n/16) 个 blocks。

　　Kernel 程序的部份，则改成：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        __shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
        const int tidc = threadIdx.x;
        const int tidr = threadIdx.y;
        const int bidc = blockIdx.x * BLOCK_SIZE;
        const int bidr = blockIdx.y * BLOCK_SIZE;
        int i, j;

        float results = 0;
        float comp = 0;

        for(j = 0; j < n; j += BLOCK_SIZE) {
          if(tidr + bidr < n && tidc + j < n) {
              matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];
            }
            else {
              matA[tidr][tidc] = 0;
            }

            if(tidr + j < n && tidc + bidc < n) {
              matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc];
            }
            else {
              matB[tidr][tidc] = 0;
            }

          __syncthreads();

            for(i = 0; i < BLOCK_SIZE; i++) {
              float t;
              comp -= matA[tidr][i] * matB[i][tidc];
                t = results - comp;
                comp = (t - results) + comp;
              results = t;
            }

          __syncthreads();
        }

        if(tidr + bidr < n && tidc + bidc < n) {
            c[(tidr + bidr) * ldc + tidc + bidc] = results;
        }
    }

　　注意到因为我们现在使用 16x16 的 threads，因此 threadIdx 变量可以取得 threadIdx.x 和 threadIdx.y，范围分别是 0 ~ 15。blockIdx.x 和 blockIdx.y 变量也是同样的情形，范围分别是 0 ~ n/16。

　　在程序中，因为矩阵的大小不一定会是 16 的倍数，因此需要使用 if 判断式检查是否超出矩阵范围。

　　这个版本在 GeForce 8800GT 上的执行结果如下：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.0780 (25.64 GFLOPS)

　　速度虽然提高了，但是似乎并没有达到预期中的八倍。当然，前面提到过，我们在计算时间时，把一些复制内存、配置内存的动作也计算在内，这些动作的时间并不会缩短。实际上 kernel 的运行时间，大约是 0.053 秒左右(约略相当于 38GFLOPS)，比上一节的版本快了将近一倍。

　　如果这一版程序已经不再限于内存带宽，那为什么没有达到预期的效率呢?这是因为这一版程序已经是限于指令周期了。除了使用 Kahan's Summation Formula 会需要更多的运算之外，程序中也有大量计算矩阵地址的乘法等等，这都会需要花费运算资源。另外，那些用来判断超出矩阵范围的 if 判断式，也会有一定的影响。

　　要把那些 if 判断式去掉，有一个方法是，在配置内存时，就配置成 16 的倍数，并在复制矩阵到显卡内存之前，先将它清为 0。如下所示：

    int newn = ((n + BLOCK_SIZE - 1) / BLOCK_SIZE) * BLOCK_SIZE;

    cudaMallocPitch((void**) &ac, &pitch_a,
        sizeof(float) * newn, newn);
   cudaMallocPitch((void**) &bc, &pitch_b,
        sizeof(float) * newn, newn);
   cudaMallocPitch((void**) &cc, &pitch_c,
        sizeof(float) * newn, newn);

   cudaMemset(ac, 0, pitch_a * newn);
   cudaMemset(bc, 0, pitch_b * newn);

　　这样一来，我们就可以把 kernel 中的 if 判断式都移除了：

    __global__ static void matMultCUDA(const float* a, size_t lda,
        const float* b, size_t ldb, float* c, size_t ldc, int n)
    {
        __shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
        const int tidc = threadIdx.x;
        const int tidr = threadIdx.y;
        const int bidc = blockIdx.x * BLOCK_SIZE;
        const int bidr = blockIdx.y * BLOCK_SIZE;
        int i, j;

        float results = 0;
        float comp = 0;

        for(j = 0; j < n; j += BLOCK_SIZE) {
          matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];
            matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc];

            __syncthreads();

          for(i = 0; i < BLOCK_SIZE; i++) {
              float t;
                comp -= matA[tidr][i] * matB[i][tidc];
                t = results - comp;
              comp = (t - results) + comp;
              results = t;
            }

          __syncthreads();
        }

        c[(tidr + bidr) * ldc + tidc + bidc] = results;
    }

　　这个版本的执行结果是：

　　Max error: 1.19209e-007 Average error: 4.22751e-008

　　Time used: 0.0780 (25.64 GFLOPS)

　　似乎没有改善。不过，实际上 kernel 的运行时间已经减少到 0.042 秒(约略相当于 48GFLOPS)。

　　结论

　　有些读者可能会想，如果把 block 再变得更大(例如 32x32)是否会有帮助呢?当然，由于最后的程序已经不再是受限于内存带宽(在 0.042 秒内存取 500MB 的数据约相当于 12GB/s 的带宽)，所以把 block 再加大并不会有帮助了。而且，由于一个 block 内的 thread 数目最多只能到 512 个，将 block 变大也会造成很多额外负担。而且 shared memory 的大小也有限制(GeForce 8800GT 的 shared memory 大小限制是 16384 bytes)，所以也不能任意增加 block 的大小。

更多内容请点击：

CUDA专区：http://cuda.it168.com/

CUDA论坛：http://cudabbs.it168.com/

第1页：cuda入门：矩阵乘法第2页：cuda入门：第一个改良第3页：cuda入门：进一步的改良和结论

关注我们