矩阵乘法:在CUDA中合并全局内存访问后性能下降

Question

我最近开始使用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_coalesced并且void MatrixMult_tiling分别是有和没有colalesced memroy访问B元素的函数.

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

编辑1:在罗伯特的回答(见下文)之后,我理解问题是在elemntwise乘法期间的加载操作.

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 ms.因此,合并的例程仍然较慢.我不明白现在的问题在哪里.

编辑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的代码将是最佳的.转置被乘数矩阵之一可能不是一个好主意.:-(

谢谢大家的帮助.

Answer 1

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

我不确定你为什么这么想:

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

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

为了确保我们不会超越术语,"合并"或"未合并"是我们应用于访问全局内存模式(不是共享内存)的术语.您的全局内存访问模式位于平铺内核的这些行中:

    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_c,如果执行替换,你会发现threadIdx.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路序列化(或者在16x16线程块的情况下可能是16路序列化),并且性能肯定更差.为什么？请记住,对于给定的载荷,j在经纱上是恒定的,并且tx在经纱上线性增加.因此,假设j在特定的循环迭代中为1.warp 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.