矩阵按行(列)求和CUDA并行算法设计

下载请点击：矩阵按行（列）求和CUDA并行算法

 

通过矩阵按行求和与按列求和两个示例介绍CUDA并行算法设计的思路，希望对大家有所帮助。很多公司CUDA工程师面试也会考察这个题目。

1.        矩阵按行求和

矩阵A[M][N]，B[M];

B[i] =A[0]+A[1]+...+A[N-1];

串行代码：

  #define A(i)(j)  A[(i)*N+(j)]

void sum_row(T *A,T*B,int M,int N)

{

     int i,j;

     T tmp;

     for(i=0; i<M; i++)

     {

              tmp=0;

              for(j=0; j<N; j++)

                       tmp += A(i)(j);

              B[i] = tmp;

     }

}

 

CUDA程序设计：

每个block处理一行数据的求和，M设置为block数目，即

dim3 dimBlock(BLOCKSIZE, 1, 1);

dim3 dimGrid(M, 1, 1);

 

__global__ void sum_row_kernel(T*A,T *B,int M,int N)

{

     int i,j;

     T tmp;

     unsigned inttid =threadIdx.x;

     i = blockIdx.x;

     tmp=0;

     T sh_B[BLOCKSIZE];

     sh_B[tid] = 0;

     __syncthreads();

     for(j=tid; j<N; j+= BLOCKSIZE)

              tmp += A(i)(j);

     sh_B[tid] = tmp;

     __syncthreads();

    if (blockSize >= 512) {if (tid < 256) {sh_B[tid] =tmp = tmp +sh_B[tid + 256]; }__syncthreads();}

    if (blockSize >= 256) {if (tid < 128) {sh_B[tid] =tmp = tmp +sh_B[tid + 128]; }__syncthreads(); }

    if (blockSize >= 128) {if (tid < 64) {sh_B[tid] =tmp = tmp +sh_B[tid +  64]; }__syncthreads(); }

   

   

    if (tid < 32)

    {

                volatile float *smem = sh_B;

                if (BLOCKSIZE >=  64) { smem[tid] = tmp = tmp + smem[tid + 32]; }

                if (BLOCKSIZE >=  32) { smem[tid] = tmp = tmp + smem[tid + 16]; }

                if (BLOCKSIZE >=  16) { smem[tid] = tmp = tmp + smem[tid +  8]; }

                if (BLOCKSIZE >=   8) { smem[tid] = tmp = tmp + smem[tid +  4]; }

                if (BLOCKSIZE >=   4) { smem[tid] = tmp = tmp + smem[tid +  2]; }

                if (BLOCKSIZE >=   2) { smem[tid] = tmp = tmp + smem[tid +  1]; }

 

  }

   if(tid==0)

       B =sh_B[0];

 

}

备注：

1）  当M比较大，N比较小时，可以采用一个warp处理一行的方式；

2）  当M比较小，N比较大时，可以让多个block处理一行的求和

以上的设计可以保证足够多的线程数和block数，充分利用GPU资源。

 

2.        矩阵按列求和

矩阵A[M][N]，B[N];

B[j] = A[0][j]+A[1][j]+...+A[M-1][j];

串行代码：

  #define A(i)(j)  A[(i)*N+(j)]

void sum_col(T *A,T*B,int M,int N)

{

     int i,j;

     T tmp;

     for(j=0; j<N; j++)

     {

              tmp=0;

              for(i=0; i<M; i++)

                       tmp += A(i)(j);

              B[j] = tmp;

     }

}

 

CUDA程序设计：

每个线程处理一列数据的求和，N/threads设置为block数目，即

dim3 dimBlock(BLOCKSIZE, 1, 1);

dim3 dimGrid(（M+ BLOCKSIZE-1）/ BLOCKSIZE, 1, 1);

 

__global__ void sum_col_kernel(T*A,T *B,int M,int N)

{

     int i,j;

     T tmp;

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

     tmp=0;

     for(i=0; i<M; i++)

              tmp += A(i)(j);  //可以满足合并访问

     B[j] = tmp;

}

备注：

当M比较大，N比较小时，可以采用一个block处理一列的方式，然后对block内的值再求和