基于CUDA的矩阵相乘

这几天研究了 一下CUDA，发现其并行的思想和普通的CPU多线程思想不太一致，但还是挺不错。主要是将任务划分成一个个block，然后每个block里面再划分成细的线程。然后每个线程做自己做的

事情。这种并行思想很适用于像矩阵运算这些元素与元素之间的运算并不耦合得很厉害，但整体数据很大的情况，这只是我对CUDA的初步感觉。
矩阵相乘的CPU程序如下：

// C = A*B
void  MatrixMulCPU( float *  _C, const   float   * _A, const   float   * _B, int  _wa, int  _ha, int  _wb)
{
     float  sum  =   0 ;
     for  ( int  i  =   0 ; i  <  _ha;  ++ i)
    {
         for  ( int  j  =   0 ; j  <  _wb;  ++ j)
        {
            sum  =   0 ;
             for  ( int  k  =   0 ; k  <  _wa;  ++ k)
            {
                sum  +=  ( float )_A[i * _wa + k] * ( float )_B[k * _wb +  j];
            }
            _C[i * _wb + j]  =  ( float )sum;
        }
    }
}

从上面可以看出，C(i,j) = sum { A(i,k)*B(k,j) } 0<=k < _wa;耦合程度很小，所以我们可以通过划分区域的方法，让每个线程负责一个区域。
怎么划分呢？首先最初的想法是让每一个线程计算一个C(i,j)，那么估算一下，应该需要height_c*width_c，也就是ha*wb个线程。进一步，我们将矩阵按一个大方格Grid划分，如果一个
方格Grid大小是16*16，那么矩阵80*48的可以表示为5(*16) * 3(*16)，即16*16个大格子(block)，每一个格子内，自然就是(height_c/16) *(width_c/16)个线程了。
好了，划分完后，内核代码如下：

计算版本0：

__global__  void  matrix_kernel_0( float *  _C, const   float *  _A, const   float   * _B, int  _wa, int  _wb)
{
     float  sum  =   0 ;
     // 找出该线程所在的行列
     int  row  =  blockIdx.y * blockDim.y  +  threadIdx.y;
     int  col  =  blockIdx.x * blockDim.x  +  threadIdx.x;

     // 线程Thread(row,col)负责计算C(row,col)
     for  ( int  i  =   0 ; i  <  _wa;  ++ i)
    {
        sum  +=  _A[row * _wa  +  i] * _B[i * _wb  +  col];
    }
    _C[row * _wb  +  col]  =  sum;
}

另外一种思路，我们不让每一个线程完整计算一个C(i,j)，通过C(i,j) = sum { A(i,k)*B(k,j) }发现，我们还可以再细度划分：
Csub(i,j) = sum{A(i,ksub+offsetA)*B(ksub+offsetB,j)}  0<=ksub < blockSize
C(i,j) = sum{Csub(i,j)}
就是把矩阵分成n*n个大的子块，然后每一个block负责计算子块i 和 子块j的子乘积，计算完毕后加起来则可。这里主要使用了共享显存作优化。

计算版本1：

__global__  void  matrix_kernel_1( float *  _C, const   float *  _A, const   float   * _B, int  _wa, int  _wb)
{
     int  bx  =  blockIdx.x;
     int  by  =  blockIdx.y;
     int  tx  =  threadIdx.x;
     int  ty  =  threadIdx.y;

     // 该block要处理的A
     int  aBegin  =  _wa * (by * BLOCK_SIZE); // A(0,by)
     int  aEnd  =  aBegin  +  _wa  -   1 ;
     int  aStep  =  BLOCK_SIZE; // offsetA

     int  bBegin  =  BLOCK_SIZE * bx; // B(bx,0)
     int  bStep  =  BLOCK_SIZE * _wb; // offsetB
    
     float  cSub  =   0 ;
     for  ( int  a  =  aBegin,b  =  bBegin; a  <=  aEnd; a  +=  aStep,b  +=  bStep)
    {
        __shared__  float  As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__  float  Bs[BLOCK_SIZE][BLOCK_SIZE];
         // 每个线程负责一个元素拷贝
        As[ty][tx]  =  _A[a  +  _wa * ty  +  tx];
        Bs[ty][tx]  =  _B[b  +  _wb * ty  +  tx];

        __syncthreads();
        
         // 每个线程负责计算一个子块i 和 子块j的子乘积
         for  ( int  k  =   0 ; k  <  BLOCK_SIZE;  ++ k)
        {
            cSub  +=  As[ty][k] * Bs[k][tx];
        }

        __syncthreads();
    }

     // 全局地址，向全局寄存器写回去
     // 一个线程负责一个元素，一个block负责一个子块
     int  cIndex  =  (by * BLOCK_SIZE  +  ty) * _wb  +  (bx * BLOCK_SIZE  +  tx);
    _C[cIndex]  =  cSub;
}

最后写一个面向Host的接口函数：

void  matrixMulGPU( float *  _C, const   float   * _A, const   float   * _B, int  _wa, int  _ha, int  _wb)
{
     float *  d_a  =  myNewOnGPU < float > (_wa * _ha);
     float *  d_b  =  myNewOnGPU < float > (_wb * _wa);
     float *  d_c  =  myNewOnGPU < float > (_wb * _ha);
    copyFromCPUToGPU(_A,d_a,_wa * _ha);
    copyFromCPUToGPU(_B,d_b,_wb * _wa);
    dim3 threads(BLOCK_SIZE,BLOCK_SIZE);
    dim3 blocks(WC / BLOCK_SIZE,HC / BLOCK_SIZE);
    matrix_kernel_0 <<< blocks,threads >>> (d_c,d_a,d_b,_wa,_wb);
    cudaThreadSynchronize();
    copyFromGPUToCPU(d_c,_C,_wb * _ha);

    myDeleteOnGPU(d_a);
    myDeleteOnGPU(d_b);
    myDeleteOnGPU(d_c);
}

调用的主函数如下：

#include  < stdio.h >
#include  < cuda_runtime.h >
#include  < cutil.h >
#include  < cutil_inline.h >
#include  < stdlib.h >
#include  < time.h >
#include  < math.h >
#include  < string .h >
#include  < Windows.h >
#include  " CUDACommon.h "
#include  " MatrixMulCPU.h "
#include  " MatrixMulGPU.h "

void  randomInit( float *  _data, int  _size)
{
     for  ( int  i  =   0 ; i  <  _size;  ++ i)
    {
        _data[i]  =  rand() / ( float )RAND_MAX;
    }
}

bool  checkError( const   float *  _A, const   float *  _B, int  _size)
{
     for  ( int  i  =   0  ; i  <  _size;  ++ i)
    {
         if  (fabs(_A[i]  -  _B[i])  >   1.0e-3 )
        {
            printf( " %f \t %f\n " ,_A[i],_B[i]);
             return   false ;
        }
    }
     return   true ;
}

int  main( int  argc,  char *  argv[])
{
    srand( 13 );
     if ( ! InitCUDA()) {
         return   0 ;
    }

     float *  A  =  myNewOnCPU < float > (WA * HA);
     float *  B  =  myNewOnCPU < float > (WB * HB);
    randomInit(A,WA * HA);
    randomInit(B,WB * HB);
     float *  C  =  myNewOnCPU < float > (WC * HC);
    memset(C, 0 , sizeof ( float ) * WC * HC);
    
     float *  C2  =  myNewOnCPU < float > (WC * HC);
    memset(C2, 0 , sizeof ( float ) * WC * HC);
    
    unsigned  int  tick1  =  GetTickCount();
    MatrixMulCPU(C2,A,B,WA,HA,WB);
    printf( " CPU use Time : %dms\n " ,GetTickCount()  -  tick1);
    unsigned  int  timer  =   0 ;
    cutilCheckError(cutCreateTimer( & timer));
    cutilCheckError(cutStartTimer(timer));
    {
        matrixMulGPU(C,A,B,WA,HA,WB);
    }
    cutilCheckError(cutStopTimer(timer));
    printf( " GPU use time: %f (ms) \n " , cutGetTimerValue(timer));
    cutilCheckError(cutDeleteTimer(timer));

     if  (checkError(C,C2,WC * HC))
    {
        printf( " Accept\n " );
    }
     else
    {
        printf( " Worng Answer\n " );
    }

    myDeleteOnCPU(A);
    myDeleteOnCPU(B);
    myDeleteOnCPU(C);
    myDeleteOnCPU(C2);

     return   0 ;
}

运算结果如下：
版本0：



版本1：

可以看出，GPU并行性能比CPU好很多，而且版本1优于版本0