首页 > 解决方案 > 矩阵乘法:在 CUDA 中合并全局内存访问后性能下降

问题描述

我最近开始使用 CUDA 处理 GPU。作为一个入门程序,我尝试有效地实现一个简单的矩阵乘法

C = AB

,从朴素的矩阵乘法开始(每个线程为 C 中的一个元素加载 A 和 B 的所有元素),平铺实现(线程在共享内存中的一个平铺中协同加载来自 A 和 B 的元素平铺以减少全局内存交通)提供良好的加速。然而,在平铺实现中,对全局内存的访问也不是按合并顺序进行的。因此,为了提高性能,最好转置矩阵 B 然后相乘。下面是我的代码,

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

#include <time.h>

#include <sys/time.h>


void querydeviceprop();
void allocate_matrix(float *h_a, float *h_b, int matDim);
void verify(float *h_c, float *h_c_check, int matDim);
void print_matrix(float *ha, int m,int n);
void transpose_matrix(float *ha, int matDim);

void mat_mul();

#define TILE_WIDTH 16 //should be equal to numThread for tiling implementation

__global__ void MatrixMult_tiling(float *d_a,float *d_b,float *d_c, int dim){

    __shared__ float ta[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
    __shared__ float tb[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
    int bx,by,tx,ty,i,j;
    float res;
    int row, col;

    bx=blockIdx.x;  by=blockIdx.y;
    tx=threadIdx.x; ty=threadIdx.y;

    row=by*TILE_WIDTH+ty;
    col=bx*TILE_WIDTH+tx;

    res=0;
    for(i=0;i<dim/TILE_WIDTH;i++){
        //collaboratively load the elements. Each thread loads a single element.
        ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
        tb[ty][tx]=d_b[(ty+i*TILE_WIDTH)*dim+col];

        __syncthreads();
        for(j=0;j<TILE_WIDTH;j++){

            res=res+ta[ty][j]*tb[j][tx];
        }
        __syncthreads();
    }
    d_c[row*dim+col]=res;
}

__global__ void MatrixMult_tiling_coalesced(float *d_a,float *d_b,float *d_c, int dim){

    __shared__ float ta[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
    __shared__ float tb[TILE_WIDTH][TILE_WIDTH]; //to load one tile of A
    int bx,by,tx,ty,i,j;
    float res;
    int row, col;

    bx=blockIdx.x;  by=blockIdx.y;
    tx=threadIdx.x; ty=threadIdx.y;

    row=by*TILE_WIDTH+ty;
    col=bx*TILE_WIDTH+tx;

    res=0;
    for(i=0;i<dim/TILE_WIDTH;i++){
        //collaboratively load the elements. Each thread loads a single element.
        ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
        tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
        __syncthreads();

        for(j=0;j<TILE_WIDTH;j++){
            res=res+ta[ty][j]*tb[tx][j];
        }
        __syncthreads();
    }
    d_c[row*dim+col]=res;
}

__global__ void MatrixMult_naive(float *d_a,float *d_b,float *d_c, int dim){

    int row,col,i;

    col=blockIdx.y*blockDim.y+threadIdx.y;
    row=blockIdx.x*blockDim.x+threadIdx.x;

    float res;

    if(row<dim && col<dim){
        res=0;
        for(i=0;i<dim;i++){
            res=res+(d_a[row*dim+i]*d_b[i*dim+col]);
        }
        d_c[row*dim+col]=res;
    }
}



int main(){
    mat_mul();
return 0;
}

void mat_mul(){

    cudaSetDevice(0);

    time_t t;
    cudaError_t err = cudaSuccess;
    srand((unsigned) time(&t));

    cudaEvent_t start, stop;
    float milliseconds=0;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);


    int matDim = 8192;

    float *h_a, *h_b, *h_c, *h_c_check;

    /*declare the host memories*/
    h_a=(float *)malloc(matDim*matDim*sizeof(float));
    h_b=(float *)malloc(matDim*matDim*sizeof(float));
    h_c=(float *)malloc(matDim*matDim*sizeof(float));
    h_c_check=(float *)malloc(matDim*matDim*sizeof(float));

    // Verify that allocations succeeded
    if (h_a == NULL || h_b == NULL || h_c == NULL || h_c_check ==NULL)
    {
    fprintf(stderr, "Failed to allocate host vectors!\n");
    exit(EXIT_FAILURE);
    }

    allocate_matrix(h_a,h_b,matDim); // allocate memory to hold the matrix

    //allocate cuda memory
        float *d_a=NULL;
        float *d_b=NULL;
        float *d_c=NULL;

        err=cudaMalloc((void **)&d_a, matDim*matDim*sizeof(float));
    if (err != cudaSuccess)
    {
        fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }
        err=cudaMalloc((void **)&d_b, matDim*matDim*sizeof(float));
    if (err != cudaSuccess)
    {
        fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }
        err=cudaMalloc((void **)&d_c, matDim*matDim*sizeof(float));
    if (err != cudaSuccess)
    {
        fprintf(stderr, "Failed to allocate device matrix A (error code %s)!\n", cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }

    printf("Matrix dimension is : %d\n",matDim);

    // Copy the host input matrix A and B in host memory to the device matrix in device memory
    //printf("Copy input data from the host memory to the CUDA device\n");

    cudaEventRecord(start);
    err = cudaMemcpy(d_a, h_a, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&milliseconds, start, stop);
    //printf("GPU memcpy matrix A %10.10f ms\n",milliseconds);


    if (err != cudaSuccess)
    {
        fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }

    cudaEventRecord(start);
    err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&milliseconds, start, stop);
    //printf("GPU memcpy matrix B %10.10f ms\n",milliseconds);

    if (err != cudaSuccess)
    {
        fprintf(stderr, "Failed to copy vector B from host to device (error code %s)!\n", cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }


    /*constants for kernel launch*/
    int numThread=16; //number of threads per Block axis
    int numBlocks=matDim/numThread;
    if(matDim%numThread)
        numBlocks++;
    dim3 dimGrid(numBlocks,numBlocks);
    dim3 dimBlock(numThread,numThread);

    //-------------------------------------------------------------
           //-------transpose and copy to GPU-------
           transpose_matrix(h_b, matDim);//transpose first the b matrix
           err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
           cudaEventSynchronize(stop);
           if (err != cudaSuccess){
                   fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
                   exit(EXIT_FAILURE);
               }
           //--------transpose and copy ends-------------

           cudaEventRecord(start);
           MatrixMult_tiling_coalesced<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
           cudaEventRecord(stop);
           err = cudaGetLastError();

           if (err != cudaSuccess){
               fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
               exit(EXIT_FAILURE);
           }

           cudaEventSynchronize(stop);
           cudaEventElapsedTime(&milliseconds, start, stop);
           printf("GPU time tiled & coalesced %10.10f ms\n",milliseconds);

           //printf("Copy output data from the CUDA device to the host memory\n");
           cudaEventRecord(start);
           err = cudaMemcpy(h_c_check, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
           cudaEventRecord(stop);
           cudaEventSynchronize(stop);
           cudaEventElapsedTime(&milliseconds, start, stop);
           //printf("GPU memcpy time tiled & coalesced %10.10f ms\n",milliseconds);

           //------------transpose back the original B matrix----------------
           transpose_matrix(h_b, matDim);//transpose first the b matrix
           err = cudaMemcpy(d_b, h_b, matDim*matDim*sizeof(float), cudaMemcpyHostToDevice);
           cudaEventSynchronize(stop);
           if (err != cudaSuccess){
                   fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
                   exit(EXIT_FAILURE);
               }

           //------------transpose back the original matrix ends-------------
//-------------------------------------------------------------

    cudaEventRecord(start);
    MatrixMult_tiling<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
    cudaEventRecord(stop);
    err = cudaGetLastError();

     if (err != cudaSuccess)
     {
         fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
         exit(EXIT_FAILURE);
     }

     cudaEventSynchronize(stop);
     cudaEventElapsedTime(&milliseconds, start, stop);
     printf("GPU time tiled %10.10f ms\n",milliseconds);

     //printf("Copy output data from the CUDA device to the host memory\n");
     cudaEventRecord(start);
     err = cudaMemcpy(h_c, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
     cudaEventRecord(stop);
     cudaEventSynchronize(stop);
     cudaEventElapsedTime(&milliseconds, start, stop);
     //printf("GPU memcpy time tiled %10.10f ms\n",milliseconds);


//-------------------------------------------------------------

    /*
    cudaEventRecord(start);
    MatrixMult_naive<<<dimGrid,dimBlock>>>(d_a, d_b, d_c, matDim);
    cudaEventRecord(stop);
    err = cudaGetLastError();

     if (err != cudaSuccess)
     {
         fprintf(stderr, "Failed to launch vector matrix kernel (error code %s)!\n", cudaGetErrorString(err));
         exit(EXIT_FAILURE);
     }

     cudaEventSynchronize(stop);
     cudaEventElapsedTime(&milliseconds, start, stop);
     printf("GPU time naive %10.10f ms\n",milliseconds);

     printf("Copy output data from the CUDA device to the host memory\n");
     cudaEventRecord(start);
     err = cudaMemcpy(h_c, d_c, matDim*matDim*sizeof(float), cudaMemcpyDeviceToHost);
     cudaEventRecord(stop);
     cudaEventSynchronize(stop);
     cudaEventElapsedTime(&milliseconds, start, stop);
     printf("GPU memcpy time naive %10.10f ms\n",milliseconds);
    */
//-------------------------------------------------------------

     verify(h_c, h_c_check, matDim);

     // Free device global memory
     err = cudaFree(d_a);

     if (err != cudaSuccess)
     {
         fprintf(stderr, "Failed to free device vector A (error code %s)!\n", cudaGetErrorString(err));
         exit(EXIT_FAILURE);
     }

     err = cudaFree(d_b);

     if (err != cudaSuccess)
     {
         fprintf(stderr, "Failed to free device vector B (error code %s)!\n", cudaGetErrorString(err));
         exit(EXIT_FAILURE);
     }

     err = cudaFree(d_c);

     if (err != cudaSuccess)
     {
         fprintf(stderr, "Failed to free device vector C (error code %s)!\n", cudaGetErrorString(err));
         exit(EXIT_FAILURE);
     }

     // Free host memory
     free(h_a);
     free(h_b);
     free(h_c);

     printf("Done\n");

}

void allocate_matrix(float *h_a, float *h_b, int matDim){

    int i,j;
    // Initialize the host input vectors
    for (i = 0; i < matDim; i++)
    {
        for(j=0;j< matDim;j++){
            h_a[i*matDim+j] = rand()%10;
            h_b[i*matDim+j] = rand()%10;
        }
    }

}

void print_matrix(float *ha, int m,int n){

    int i, j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            printf("  %.1f ",ha[i*m+j]);
        }
    printf("\n");
    }
}

void transpose_matrix(float *h_a, int matDim){

    int i, j;
    int temp;
    for(i=0;i<matDim;i++)
    {
    for(j=0;j<i;j++)
    {
        temp=h_a[i*matDim+j];
        h_a[i*matDim+j]=h_a[j*matDim+i];
        h_a[j*matDim+i]=temp;
    }
    }
}


void verify(float *h_c, float *h_c_check, int matDim){

    int i,j;
    //check the code
     for (i = 0; i < matDim; i++)
     {
         for(j=0;j<matDim;j++){
         if (fabs(h_c[i*matDim+j] - h_c_check[i*matDim+j]) > 1e-5)
         {
             printf("cpu : %f , gpu : %f\t",h_c[i*matDim+j],h_c_check[i*matDim+j]);
             fprintf(stderr, "Result verification failed at element %d,%d !\n\n", i,j);
             exit(EXIT_FAILURE);
         }

         }
     }

     printf("Test PASSED\n");

}

MatrixMult_tiling_coalescedvoid MatrixMult_tiling分别是对 B 的元素有和没有合并内存访问的函数。

现在,问题是所花费的时间MatrixMult_tiling_coalesced几乎是所花费时间的两倍MatrixMult_tiling。我知道在MatrixMult_tiling元素中以合并的方式(即按行主要顺序)为每个图块加载到图块中,但是图块沿列排列,而图块中的图块MatrixMult_tiling_coalesced沿行排列,因此该功能MatrixMult_tiling_coalesced应该比另一个更快。但在实践中,我可以看到相反的情况。如果有人能指出原因,我将不胜感激。提前谢谢..

编辑1:在罗伯特回答之后(见下文),我明白问题出在元素乘法期间的加载操作中。

tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];]

tb[tx][ty]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];

res=res+ta[ty][j]*tb[tx][j];

res=res+ta[ty][j]*tb[j][tx];

这将MatrixMult_tiling_coalesced函数的性能从 1500 毫秒提高到 1000 毫秒。但是,该函数MatrixMult_tiling只需要 879 毫秒。因此,合并后的例程仍然较慢。我不明白现在问题出在哪里。

编辑 2:我意识到在编辑 1 中,我刚刚将银行冲突问题从元素乘法移到了元素加载部分。代码中的以下更改没有银行冲突,

    tb[tx][ty]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];

    tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];

res=res+ta[ty][j]*tb[j][tx];

res=res+ta[ty][j]*tb[ty][j];

但它仍然比功能稍慢MatrixMult_tiling。该MatrixMult_tiling_coalesced函数需要 982 ms 与 870 ms 的MatrixMult_tiling函数。它应该至少类似于MatrixMult_tiling如果不是更快的话。

最后编辑:

编辑 2 不会产生正确的结果。所以编辑 1 的代码将是最佳的。转置被乘数矩阵之一可能不是一个好主意。:-(

感谢大家的帮助。

标签: parallel-processingcudagpumatrix-multiplication

解决方案


B当然不是我要转置的矩阵C=AB。但这既不是这里也不是那里。

我不确定你为什么认为:

在平铺实现中,对全局内存的访问也不是合并顺序

我没有在您的代码中看到任何MatrixMult_tiling导致未合并访问的代码行。

只是为了确保我们不会被术语绊倒,“coalesced”或“uncoalesced”是我们应用于访问全局内存(不是共享内存)的模式的术语。您的全局内存访问模式位于平铺内核中的这些行中:

    ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
    tb[ty][tx]=d_b[(ty+i*TILE_WIDTH)*dim+col];
    ...
    d_c[row*dim+col]=res;

并且这些模式都没有被合并到全局内存中。在每个生成的索引d_a中,如果执行替换,您会发现变量存在于所有索引中,d_b并且不乘以任何值,常量或其他值。因此,这些模式都会(很好地)合并​​。d_cthreadIdx.x

如果有人能指出原因,我将不胜感激。

当涉及到共享内存时,您做了一些坏事。

在您的平铺内核中,您的乘法运算如下所示:

        res=res+ta[ty][j]*tb[j][tx];

对于这种情况:

   ta[ty][j]

我们有一种情况,warp 中的所有线程(具有线性增加tx的值但相同的ty值)正在读取共享内存中的相同位置。这是一种“最佳”访问模式——它不会出现任何银行冲突,并且会在尽可能短的时间内得到服务。

对于这种情况:

  tb[j][tx]

我们有一种情况,warp 中的相邻线程正在读取共享内存中的相邻位置。这也是一种“最佳”的、非银行冲突的模式,并将在尽可能短的时间内得到服务。

但是在您的MatrixMult_tiling_coalesced内核中,相应的乘法运算是:

  res=res+ta[ty][j]*tb[tx][j];

同样,在这种情况下:

  ta[ty][j]

我们有一个共享内存“广播”模式(warp 中的所有线程从同一位置读取),这是最佳且快速的。但在这种情况下:

  tb[tx][j]

您实际上已经创建了对共享内存的列访问。这是共享内存最糟糕的访问模式,它将导致加载过程的 32 路序列化(或者可能是 16 路序列化,在您的 16x16 线程块的情况下),并且性能肯定更差。为什么?请记住,对于给定的负载,j在整个经线中是恒定的,并且tx在整个经线中线性增加。因此,假设j在特定循环迭代中为 1。经线 0 中的线程将读取:

tb[0][1], tb[1][1], tb[2][1], tb[3][1], ...

并且这些位置都属于共享内存的特定“列”,即它们都属于同一个共享内存库。这是共享内存的最坏情况模式。

为了完整起见,我声称MatrixMult_tiling_coalesced内核中的所有全局内存访问模式也被合并:

    ta[ty][tx]=d_a[row*dim+TILE_WIDTH*i+tx];
    tb[ty][tx]=d_b[bx*TILE_WIDTH*dim + TILE_WIDTH*i+ty*dim+tx];
    ...
    d_c[row*dim+col]=res;

因此,您的两个内核实现之间的全局内存访问模式/活动/效率应该没有重大差异。

作为旁注,我认为这都是一个学习练习。如果您对 GPU 上的高性能矩阵乘法感兴趣,我建议您考虑使用CUBLAS


推荐阅读