CUDA/GPU下矩阵乘法的几种实现的C++源码——转载

 

[cpp] view plaincopy
 
  1. #include "cutil_inline.h"  
  2. #include "cublas.h"  
  3.   
  4.   
  5. void MatrixMul(const float *A, const float *B, float *C, int Width) {  
  6.     int i, j, k;  
  7.     for(i=0; i<Width; i++)  
  8.         for(j=0; j<Width; j++){  
  9.             float s=0;  
  10.             for(k=0; k<Width; k++)   
  11.                 s+=A[i*Width+k]*B[k*Width+j];  
  12.             C[i*Width+j]=s;  
  13.         }  
  14. }  
  15.   
  16. #define TILE_WIDTH 16  
  17. //简单方法  
  18. __global__ void MatrixMulKernel1(float* Md, float* Nd, float* Pd, int Width)  
  19. {  
  20.     int x = threadIdx.x+blockIdx.x*blockDim.x;  
  21.     int y = threadIdx.y+blockIdx.y*blockDim.y;  
  22.   
  23.     float Pvalue = 0;  
  24.     for (int k = 0; k < Width; ++k)  
  25.         Pvalue+=Md[y * Width + k]*Nd[k * Width + x];  
  26.     Pd[y*Width + x] = Pvalue;  
  27. }  
  28.   
  29. //块方法  
  30. __global__ void MatrixMulKernel2(float* Md, float* Nd, float* Pd, int Width)  
  31. {  
  32.     int bx = blockIdx.x;  
  33.     int by = blockIdx.y;  
  34.     int tx = threadIdx.x;  
  35.     int ty = threadIdx.y;  
  36.     float Pvalue = 0;  
  37.       
  38.     for (int m = 0; m < gridDim.x; ++m) {  
  39.         __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];  
  40.         __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];  
  41.         Mds[ty][tx] = *(Md + (by*blockDim.y + ty) * Width + m*blockDim.x + tx);  
  42.         Nds[ty][tx] = *(Nd + (m*blockDim.y + ty) * Width + bx*blockDim.x + tx);  
  43.         __syncthreads();  
  44.         for (int k = 0; k < blockDim.x; ++k)  
  45.             Pvalue += Mds[ty][k] * Nds[k][tx];  
  46.         __syncthreads();  
  47.     }  
  48.     Pd[(by*blockDim.y+ty)*Width+bx*blockDim.x+tx] = Pvalue;  
  49. }  
  50. //块方法, 循环展开  
  51. __global__ void MatrixMulKernel3(float* Md, float* Nd, float* Pd, int Width)  
  52. {  
  53.     int bx = blockIdx.x;  
  54.     int by = blockIdx.y;  
  55.     int tx = threadIdx.x;  
  56.     int ty = threadIdx.y;  
  57.     float Pvalue = 0;  
  58.       
  59.     for (int m = 0; m < gridDim.x; ++m) {  
  60.         __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];  
  61.         __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];  
  62.         Mds[ty][tx] = *(Md + (by*blockDim.y + ty) * Width + m*blockDim.x + tx);  
  63.         Nds[ty][tx] = *(Nd + (m*blockDim.y + ty) * Width + bx*blockDim.x + tx);  
  64.         __syncthreads();  
  65.             Pvalue += Mds[ty][0] * Nds[0][tx] + Mds[ty][1] * Nds[1][tx] +  
  66.                 Mds[ty][2] * Nds[2][tx] + Mds[ty][3] * Nds[3][tx] +  
  67.                 Mds[ty][4] * Nds[4][tx] + Mds[ty][5] * Nds[5][tx] +  
  68.                 Mds[ty][6] * Nds[6][tx] + Mds[ty][7] * Nds[7][tx] +  
  69.                 Mds[ty][8] * Nds[8][tx] + Mds[ty][9] * Nds[9][tx] +  
  70.                 Mds[ty][10] * Nds[10][tx] + Mds[ty][11] * Nds[11][tx] +  
  71.                 Mds[ty][12] * Nds[12][tx] + Mds[ty][13] * Nds[13][tx] +  
  72.                 Mds[ty][14] * Nds[14][tx] + Mds[ty][15] * Nds[15][tx];  
  73.         __syncthreads();  
  74.     }  
  75.     Pd[(by*blockDim.y+ty)*Width+bx*blockDim.x+tx] = Pvalue;  
  76. }  
  77.   
  78. #define ThreadColumn 2  
  79. //块方法, 线程粒度  
  80. __global__ void MatrixMulKernel4(float* Md, float* Nd, float* Pd, int Width)  
  81. {  
  82.     int bx = blockIdx.x;  
  83.     int by = blockIdx.y;  
  84.     int tx = threadIdx.x;  
  85.     int ty = threadIdx.y;  
  86.     float Pvalue1 = 0, Pvalue2=0;  
  87.       
  88.     for (int m = 0; m < gridDim.x; ++m) {  
  89.         __shared__ float Mds[TILE_WIDTH][ThreadColumn*TILE_WIDTH];  
  90.         __shared__ float Nds[ThreadColumn*TILE_WIDTH][ThreadColumn*TILE_WIDTH];  
  91.         for(int k=0;k<ThreadColumn; ++k)  
  92.             Mds[ty][tx+k*TILE_WIDTH] = *(Md + (by*blockDim.y + ty) * Width +  ThreadColumn*m*blockDim.x + tx + k*TILE_WIDTH );  
  93.         for(int h=0;h<ThreadColumn; ++h)  
  94.             for(int k=0; k<ThreadColumn; ++k)  
  95.                 Nds[ty+h*TILE_WIDTH][tx+k*TILE_WIDTH] = *(Nd + (m*blockDim.y*ThreadColumn + ty + h*TILE_WIDTH) * Width  +   
  96.                   bx*blockDim.x*ThreadColumn + tx + k*TILE_WIDTH);  
  97.         __syncthreads();  
  98.         for(int h=0;h<ThreadColumn*TILE_WIDTH;++h)  
  99.             Pvalue1 +=Mds[ty][h] * Nds[h][tx];  
  100.         for(int k=0;k<ThreadColumn*TILE_WIDTH;++k)  
  101.             Pvalue2 +=Mds[ty][k] * Nds[k][tx+TILE_WIDTH];  
  102.         __syncthreads();  
  103.     }  
  104.       
  105.     Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx] = Pvalue1;  
  106.     Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx + TILE_WIDTH] = Pvalue2;  
  107. }  
  108. //块方法, 循环展开, 线程粒度  
  109. __global__ void MatrixMulKernel5(float* Md, float* Nd, float* Pd, int Width)  
  110. {  
  111.     int bx = blockIdx.x;  
  112.     int by = blockIdx.y;  
  113.     int tx = threadIdx.x;  
  114.     int ty = threadIdx.y;  
  115.     float Pvalue1 = 0, Pvalue2=0;  
  116.       
  117.     for (int m = 0; m < gridDim.x; ++m) {  
  118.         __shared__ float Mds[TILE_WIDTH][ThreadColumn*TILE_WIDTH];  
  119.         __shared__ float Nds[ThreadColumn*TILE_WIDTH][ThreadColumn*TILE_WIDTH];  
  120.         for(int k=0;k<ThreadColumn; ++k)  
  121.             Mds[ty][tx+k*TILE_WIDTH] = *(Md + (by*blockDim.y + ty) * Width +  ThreadColumn*m*blockDim.x + tx + k*TILE_WIDTH );  
  122.         for(int h=0;h<ThreadColumn; ++h)  
  123.             for(int k=0; k<ThreadColumn; ++k)  
  124.                 Nds[ty+h*TILE_WIDTH][tx+k*TILE_WIDTH] = *(Nd + (m*blockDim.y*ThreadColumn + ty + h*TILE_WIDTH) * Width  +   
  125.                   bx*blockDim.x*ThreadColumn + tx + k*TILE_WIDTH);  
  126.         __syncthreads();  
  127.             Pvalue1 += Mds[ty][0] * Nds[0][tx] + Mds[ty][1] * Nds[1][tx] +  
  128.                 Mds[ty][2] * Nds[2][tx] + Mds[ty][3] * Nds[3][tx] +  
  129.                 Mds[ty][4] * Nds[4][tx] + Mds[ty][5] * Nds[5][tx] +  
  130.                 Mds[ty][6] * Nds[6][tx] + Mds[ty][7] * Nds[7][tx] +  
  131.                 Mds[ty][8] * Nds[8][tx] + Mds[ty][9] * Nds[9][tx] +  
  132.                 Mds[ty][10] * Nds[10][tx] + Mds[ty][11] * Nds[11][tx] +  
  133.                 Mds[ty][12] * Nds[12][tx] + Mds[ty][13] * Nds[13][tx] +  
  134.                 Mds[ty][14] * Nds[14][tx] + Mds[ty][15] * Nds[15][tx] +  
  135.                 Mds[ty][16] * Nds[16][tx] + Mds[ty][17] * Nds[17][tx] +  
  136.                 Mds[ty][18] * Nds[18][tx] + Mds[ty][19] * Nds[19][tx] +  
  137.                 Mds[ty][20] * Nds[20][tx] + Mds[ty][21] * Nds[21][tx] +  
  138.                 Mds[ty][22] * Nds[22][tx] + Mds[ty][23] * Nds[23][tx] +  
  139.                 Mds[ty][24] * Nds[24][tx] + Mds[ty][25] * Nds[25][tx] +  
  140.                 Mds[ty][26] * Nds[26][tx] + Mds[ty][27] * Nds[27][tx] +  
  141.                 Mds[ty][28] * Nds[28][tx] + Mds[ty][29] * Nds[29][tx] +  
  142.                 Mds[ty][30] * Nds[30][tx] + Mds[ty][31] * Nds[31][tx];  
  143.   
  144.             Pvalue2 += Mds[ty][0] * Nds[0][tx+TILE_WIDTH] + Mds[ty][1] * Nds[1][tx+TILE_WIDTH] +  
  145.                 Mds[ty][2] * Nds[2][tx+TILE_WIDTH] + Mds[ty][3] * Nds[3][tx+TILE_WIDTH] +  
  146.                 Mds[ty][4] * Nds[4][tx+TILE_WIDTH] + Mds[ty][5] * Nds[5][tx+TILE_WIDTH] +  
  147.                 Mds[ty][6] * Nds[6][tx+TILE_WIDTH] + Mds[ty][7] * Nds[7][tx+TILE_WIDTH] +  
  148.                 Mds[ty][8] * Nds[8][tx+TILE_WIDTH] + Mds[ty][9] * Nds[9][tx+TILE_WIDTH] +  
  149.                 Mds[ty][10] * Nds[10][tx+TILE_WIDTH] + Mds[ty][11] * Nds[11][tx+TILE_WIDTH] +  
  150.                 Mds[ty][12] * Nds[12][tx+TILE_WIDTH] + Mds[ty][13] * Nds[13][tx+TILE_WIDTH] +  
  151.                 Mds[ty][14] * Nds[14][tx+TILE_WIDTH] + Mds[ty][15] * Nds[15][tx+TILE_WIDTH] +  
  152.                 Mds[ty][16] * Nds[16][tx+TILE_WIDTH] + Mds[ty][17] * Nds[17][tx+TILE_WIDTH] +  
  153.                 Mds[ty][18] * Nds[18][tx+TILE_WIDTH] + Mds[ty][19] * Nds[19][tx+TILE_WIDTH] +  
  154.                 Mds[ty][20] * Nds[20][tx+TILE_WIDTH] + Mds[ty][21] * Nds[21][tx+TILE_WIDTH] +  
  155.                 Mds[ty][22] * Nds[22][tx+TILE_WIDTH] + Mds[ty][23] * Nds[23][tx+TILE_WIDTH] +  
  156.                 Mds[ty][24] * Nds[24][tx+TILE_WIDTH] + Mds[ty][25] * Nds[25][tx+TILE_WIDTH] +  
  157.                 Mds[ty][26] * Nds[26][tx+TILE_WIDTH] + Mds[ty][27] * Nds[27][tx+TILE_WIDTH] +  
  158.                 Mds[ty][28] * Nds[28][tx+TILE_WIDTH] + Mds[ty][29] * Nds[29][tx+TILE_WIDTH] +  
  159.                 Mds[ty][30] * Nds[30][tx+TILE_WIDTH] + Mds[ty][31] * Nds[31][tx+TILE_WIDTH];  
  160.         __syncthreads();  
  161.     }  
  162.       
  163.     Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx] = Pvalue1;  
  164.     Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx + TILE_WIDTH] = Pvalue2;  
  165. }  
  166. int main(int argc, char **argv) {  
  167.     int Width(0);  
  168.     bool is_test=false;  
  169.     if(argc==1){  
  170.         Width=512;  
  171.         is_test=false;  
  172.     }  
  173.     else if(argc==2){  
  174.         Width=atoi(argv[1]);  
  175.         is_test=false;  
  176.     }  
  177.     else{  
  178.         Width=atoi(argv[1]);  
  179.         is_test=(bool)atoi(argv[2]);  
  180.     }  
  181.     float *h_A=(float*)malloc(Width*Width*sizeof(float));  
  182.     float *h_B=(float*)malloc(Width*Width*sizeof(float));  
  183.     float *h_C=(float*)malloc(Width*Width*sizeof(float));  
  184.     float *h_C_ref=(float*)malloc(Width*Width*sizeof(float));  
  185.     float *d_A, *d_B, *d_C;  
  186.     float t0, error_norm=0, ref_norm=0;  
  187.     unsigned int timer1=0;  
  188.     cutCreateTimer(&timer1);  
  189.     cutStartTimer(timer1);  
  190.     for(int i=0; i<Width*Width; i++) {  
  191.         h_A[i]=rand()/(float)RAND_MAX;  
  192.         h_B[i]=rand()/(float)RAND_MAX;  
  193.     }  
  194.     cublasInit();  
  195.     cublasAlloc(Width*Width, sizeof(float), (void**)&d_A);  
  196.     cublasAlloc(Width*Width, sizeof(float), (void**)&d_B);  
  197.     cublasAlloc(Width*Width, sizeof(float), (void**)&d_C);  
  198.     cublasSetVector(Width*Width, sizeof(float), h_A, 1, d_A, 1);  
  199.     cublasSetVector(Width*Width, sizeof(float), h_B, 1, d_B, 1);  
  200.   
  201.     if(is_test)  
  202.         MatrixMul(h_A, h_B, h_C_ref,Width);  
  203.   
  204.     dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);  
  205.     dim3 dimGrid(Width/TILE_WIDTH, Width/TILE_WIDTH);  
  206.       
  207.     t0=cutGetTimerValue(timer1);  
  208.     cudaThreadSynchronize();  
  209.     MatrixMulKernel1<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);  
  210.     cudaThreadSynchronize();  
  211.     float gpu_t1=(cutGetTimerValue(timer1)-t0)/1000.0f;  
  212.     cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);  
  213.   
  214.     if(is_test){  
  215.         error_norm=0;  
  216.         ref_norm=0;  
  217.         for(int i=0; i<Width*Width; i++) {  
  218.             float diff=h_C_ref[i]-h_C[i];  
  219.             error_norm+=diff*diff;  
  220.             ref_norm+=h_C_ref[i]*h_C_ref[i];  
  221.         }  
  222.         printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");  
  223.     }  
  224.     t0=cutGetTimerValue(timer1);  
  225.     cudaThreadSynchronize();  
  226.     MatrixMulKernel2<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);  
  227.     cudaThreadSynchronize();  
  228.     float gpu_t2=(cutGetTimerValue(timer1)-t0)/1000.0f;  
  229.     cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);  
  230.     if(is_test){  
  231.         error_norm=0;  
  232.         ref_norm=0;  
  233.         for(int i=0; i<Width*Width; i++) {  
  234.             float diff=h_C_ref[i]-h_C[i];  
  235.             error_norm+=diff*diff;  
  236.             ref_norm+=h_C_ref[i]*h_C_ref[i];  
  237.         }  
  238.         printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");  
  239.     }  
  240.   
  241.     t0=cutGetTimerValue(timer1);  
  242.     cudaThreadSynchronize();  
  243.     MatrixMulKernel3<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);  
  244.     cudaThreadSynchronize();  
  245.     float gpu_t3=(cutGetTimerValue(timer1)-t0)/1000.0f;  
  246.     cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);  
  247.     if(is_test){  
  248.         error_norm=0;  
  249.         ref_norm=0;  
  250.         for(int i=0; i<Width*Width; i++) {  
  251.             float diff=h_C_ref[i]-h_C[i];  
  252.             error_norm+=diff*diff;  
  253.             ref_norm+=h_C_ref[i]*h_C_ref[i];  
  254.         }  
  255.         printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");  
  256.     }  
  257.   
  258.     dimGrid.x/=ThreadColumn;  
  259.   
  260.     t0=cutGetTimerValue(timer1);  
  261.     cudaThreadSynchronize();  
  262.     MatrixMulKernel4<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);  
  263.     cudaThreadSynchronize();  
  264.     float gpu_t4=(cutGetTimerValue(timer1)-t0)/1000.0f;  
  265.     cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);  
  266.     if(is_test){  
  267.         error_norm=0;  
  268.         ref_norm=0;  
  269.         for(int i=0; i<Width*Width; i++) {  
  270.             float diff=h_C_ref[i]-h_C[i];  
  271.             error_norm+=diff*diff;  
  272.             ref_norm+=h_C_ref[i]*h_C_ref[i];  
  273.         }  
  274.         printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");  
  275.     }  
  276.   
  277.     t0=cutGetTimerValue(timer1);  
  278.     cudaThreadSynchronize();  
  279.     MatrixMulKernel5<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);  
  280.     cudaThreadSynchronize();  
  281.     float gpu_t5=(cutGetTimerValue(timer1)-t0)/1000.0f;  
  282.     cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);  
  283.     if(is_test){  
  284.         error_norm=0;  
  285.         ref_norm=0;  
  286.         for(int i=0; i<Width*Width; i++) {  
  287.             float diff=h_C_ref[i]-h_C[i];  
  288.             error_norm+=diff*diff;  
  289.             ref_norm+=h_C_ref[i]*h_C_ref[i];  
  290.         }  
  291.         printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");  
  292.     }  
  293.   
  294.     printf("矩阵阶数为%4d,", Width);  
  295.     printf("简单方法: %.6fs(%.3fGflops),", gpu_t1, 1e-9*Width*Width*Width*2/gpu_t1);  
  296.     printf("块方法: %.6fs(%.3fGflops),", gpu_t2, 1e-9*Width*Width*Width*2/gpu_t2);  
  297.     printf("块+循环展开方法: %.6fs(%.3fGflops),", gpu_t3, 1e-9*Width*Width*Width*2/gpu_t3);  
  298.     printf("块+线程粒度2: %.6fs(%.3fGflops),", gpu_t4, 1e-9*Width*Width*Width*2/gpu_t4);  
  299.     printf("块+循环展开方法+线程粒度2: %.6fs(%.3fGflops)/n", gpu_t5, 1e-9*Width*Width*Width*2/gpu_t5);  
  300. }  

 

 

环境：CUDA toolkit3.2+Windows XP+CUDA SDK中的vs2008模板release编译通过，显卡是GeForce GT240。大家可根据自己的情况进行测试，报告一下结果吧.


不同规模的运行结果:
矩阵阶数为 512,简单方法: 0.033887s(7.922Gflops),块方法: 0.008424s(31.865Gflops),块+循环展开方法: 0.003995s(67.191Gflops),块+线程粒度2: 0.008091s(33.179Gflops),块+循环展开方法+线程粒度2: 0.003399s(78.984Gflops)

矩阵阶数为1024,简单方法: 0.264424s(8.121Gflops),块方法: 0.066201s(32.439Gflops),块+循环展开方法: 0.030827s(69.662Gflops),块+线程粒度2: 0.063264s(33.945Gflops),块+循环展开方法+线程粒度2: 0.026282s(81.710Gflops)

矩阵阶数为2048,简单方法: 2.112513s(8.132Gflops),块方法: 0.527002s(32.599Gflops),块+循环展开方法: 0.244058s(70.392Gflops),块+线程粒度2: 0.505495s(33.986Gflops),块+循环展开方法+线程粒度2: 0.208845s(82.261Gflops)

矩阵阶数为4096,简单方法: 17.705070s(7.763Gflops),块方法: 4.205308s(32.682Gflops),块+循环展开方法: 1.966043s(69.906Gflops),块+线程粒度2: 4.059064s(33.860Gflops),块+循环展开方法+线程粒度2: 1.677771s(81.918Gflops) 

 

测试结果表明在块大小16x16时，块方法和循环展开对速度有很显著的影响，而线程粒度的使用对速度只有很小的提高。而块大小是8x8时的情形没有测试，线程粒度对速度应该有较大影响，家可以做一做。

posted @ 2015-11-25 15:38  uestc_summer  阅读(646)  评论(0)    收藏  举报