【CUDA并行编程之四】矩阵相乘


前面介绍了基本的Cuda编程的相关知识，那么这一篇在此基础之上来看看GPU在处理数据计算上的高效能，我们拿矩阵相乘来作为例子。


1.CPU上执行矩阵相乘以及性能。

在CPU上进行矩阵相乘运算的代码：

mat_mul.cc:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. //a[i]*b[i] + c[i] = d[i]  
  2. #include<iostream>  
  3. #include<vector>  
  4. #include<map>  
  5. #include<fstream>  
  6. #include"wtime.h"   
  7.   
  8. using namespace std;  
  9.   
  10. const int N = 320;  
  11.   
  12. //矩阵有两种表达的方法用二维矩阵或者用一维矩阵表示  
  13. int a[N+1][N+1],b[N+1][N+1],c[N+1][N+1],d[N+1][N+1];  
  14. int aa[(N+1)*(N+1)],bb[(N+1)*(N+1)],cc[(N+1)*(N+1)],dd[(N+1)*(N+1)];  
  15.   
  16. void init()  
  17. {  
  18.     for(int i=0;i<N;i++)  
  19.         for(int j=0;j<N;j++)  
  20.         {  
  21.             a[i][j] = 1;  
  22.             b[i][j] = 2;  
  23.             c[i][j] = 3;  
  24.         }  
  25. }  
  26.   
  27. void init1()  
  28. {  
  29.     for(int i=0;i<N;i++)  
  30.         for(int j=0;j<N;j++)  
  31.         {  
  32.             aa[i*N+j] = 1;  
  33.             bb[i*N+j] = 2;  
  34.             cc[i*N+j] = 3;  
  35.         }  
  36. }  
  37.   
  38. void mul()  
  39. {  
  40.     for(int i=0;i<N;i++)   
  41.       for(int j=0;j<N;j++)  
  42.       {  
  43.         for(int k=0;k<N;k++)  
  44.         {  
  45.             d[i][j] += a[i][k] * b[k][j];  
  46.         }  
  47.         d[i][j] += c[i][j];  
  48.       }  
  49. }  
  50.   
  51. void mul1()  
  52. {  
  53.     for(int i=0;i<N;i++)   
  54.       for(int j=0;j<N;j++)  
  55.       {  
  56.         for(int k=0;k<N;k++)  
  57.         {  
  58.             dd[i*N+j] += aa[i*N+k] * bb[k*N+j];  
  59.         }  
  60.         dd[N*i+j] += cc[N*i+j];  
  61.       }  
  62. }  
  63.   
  64. void print()  
  65. {  
  66.     ofstream fout;  
  67.     fout.open("result.txt");  
  68.     if(!fout)  
  69.     {  
  70.         perror("can not open the file");  
  71.     }  
  72.     for(int i=0;i<N;i++)  
  73.     {  
  74.       for(int j=0;j<N;j++)  
  75.       {  
  76.           fout<<d[i][j]<<" ";  
  77.       }  
  78.       fout<<endl;  
  79.     }  
  80.     fout.close();  
  81. }  
  82.   
  83. int main()  
  84. {  
  85.     init1();      
  86.       
  87.     double t = wtime();  
  88.     mul1();  
  89.     t = wtime()-t;  
  90.     printf("computation timing = %10.10f sec ",t);  
  91.       
  92.     //print();  
  93.   
  94.     return 0;  
  95. }  
wtime.h:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. #ifndef _WTIME_  
  2. #define _WTIME_  
  3.   
  4. double wtime();  
  5.   
  6. #endif  


wtime.cc:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. #include <stdio.h>  
  2. #include <sys/time.h>  
  3. #include <iostream>  
  4. #include <cstdlib>  
  5.   
  6. double wtime(void)  
  7. {  
  8.     double now_time;  
  9.     struct timeval etstart;  
  10.     struct timezone tzp;  
  11.   
  12.     if(gettimeofday(&etstart,&tzp)==-1)  
  13.     {  
  14.         perror("Error:calling gettimeofday() not successfully. ");  
  15.     }  
  16.   
  17.     now_time = ( (double)etstart.tv_sec ) + ((double)etstart.tv_usec) / 1000000.0;  
  18.   
  19.     return now_time;  
  20. }  
  21.   
  22. #if 0  
  23. int main()  
  24. {  
  25.     double time;  
  26.     time = wtime();  
  27.   
  28.     printf("time of day = %10.4f ",time);  
  29.   
  30.     return 0;  
  31. }  
  32. #endif  

makefile:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. target:  
  2.     g++ mat_mul.cc wtime.cc  
  3.     ./a.out  

结果：



2.GPU上执行矩阵相乘以及性能。

代码：

cuda_mat_mul_v1.cu:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. //matrix multiplication with global memory   
  2. #include<iostream>  
  3. #include<fstream>  
  4. #include "wtime.h"  
  5.   
  6. using namespace std;  
  7.   
  8.   
  9. const int BLOCK_SIZE = 16;  
  10. const int GRID_SIZE = 20;  
  11.   
  12. //D = A * B + C;  
  13. __global__ void mat_mul(int *da,int *db,int *dc,int *dd,int N)  
  14. {  
  15.     int row = blockIdx.y * blockDim.y + threadIdx.y;  
  16.     int col = blockIdx.x * blockDim.x + threadIdx.x;  
  17.   
  18.     int sum = 0;  
  19.     for(int i=0;i<N;i++)  
  20.     {  
  21.         sum += da[row*N + i] * db[row*i+col];  
  22.     }  
  23.     dd[row*N + col] = sum + dc[row*N + col];  
  24. }  
  25.   
  26. int main()  
  27. {  
  28.     int N = BLOCK_SIZE * GRID_SIZE;  
  29.     int *ha,*hb,*hc,*hd;  
  30.     int *da,*db,*dc,*dd;  
  31.     double time;  
  32.     ha = new int[N*N];  
  33.     hb = new int[N*N];  
  34.     hc = new int[N*N];  
  35.     hd = new int[N*N];  
  36.     cudaError_t err;  
  37.   
  38.     //initialize  
  39.     for(int i=0;i<N;i++)  
  40.         for(int j=0;j<N;j++)  
  41.         {  
  42.             ha[i*N+j] = 1;  
  43.             hb[i*N+j] = 2;  
  44.             hc[i*N+j] = 3;  
  45.         }  
  46.   
  47.     //malloc</strong>  
  48.     cudaMalloc(&da,N*N*sizeof(int));  
  49.     cudaMalloc(&db,N*N*sizeof(int));  
  50.     cudaMalloc(&dc,N*N*sizeof(int));  
  51.     err = cudaMalloc(&dd,N*N*sizeof(int));  
  52.     printf("Cuda Malloc C : %s ",cudaGetErrorString(err));  
  53.   
  54.     //host to device  
  55.     cudaMemcpy(da,ha,N*N*sizeof(int),cudaMemcpyHostToDevice);  
  56.     cudaMemcpy(db,hb,N*N*sizeof(int),cudaMemcpyHostToDevice);  
  57.     cudaMemcpy(dc,hc,N*N*sizeof(int),cudaMemcpyHostToDevice);  
  58.     cudaMemcpy(dd,hd,N*N*sizeof(int),cudaMemcpyHostToDevice);  
  59.   
  60.     dim3 threadBlock(BLOCK_SIZE,BLOCK_SIZE);  
  61.     dim3 grid(GRID_SIZE,GRID_SIZE);  
  62.     //kernel  
  63.     time = wtime();  
  64.     mat_mul<<<grid,threadBlock>>>(da,db,dc,dd,N);  
  65.     printf("Computation time is %10.10f ",wtime()-time);  
  66.   
  67.     //device to host  
  68.     cudaMemcpy(hd,dd,N*N*sizeof(int),cudaMemcpyDeviceToHost);  
  69.   
  70.     //print result to file  
  71.     ofstream fout;  
  72.     fout.open("result_v1.txt");  
  73.     if(!fout)    
  74.     {  
  75.         cerr<<"open the file error"<<endl;  
  76.         exit(-1);  
  77.     }  
  78.     for(int i=0;i<N;i++)   
  79.     {  
  80.         for(int j=0;j<N;j++)  
  81.         {  
  82.             fout<<hd[i*N+j]<<" ";  
  83.         }  
  84.         fout<<endl;  
  85.     }  
  86.       
  87.     delete []ha;delete []hb;delete []hc;delete []hd;  
  88.     cudaFree(da);cudaFree(db);cudaFree(dc);cudaFree(dd);  
  89.   
  90.     return 0;  
  91. }  
cuda_wtime.cu:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. #include <stdio.h>  
  2. #include <sys/time.h>  
  3. #include <iostream>  
  4. #include <cstdlib>  
  5.   
  6. double wtime(void)  
  7. {  
  8.     double now_time;  
  9.     struct timeval etstart;  
  10.     struct timezone tzp;  
  11.   
  12.     if(gettimeofday(&etstart,&tzp)==-1)  
  13.     {  
  14.         perror("Error:calling gettimeofday() not successfully. ");  
  15.     }  
  16.   
  17.     now_time = ( (double)etstart.tv_sec ) + ((double)etstart.tv_usec) / 1000000.0;  
  18.   
  19.     return now_time;  
  20. }  
  21.   
  22. #if 0  
  23. int main()  
  24. {  
  25.     double time;  
  26.     time = wtime();  
  27.   
  28.     printf("time of day = %10.4f ",time);  
  29.   
  30.     return 0;  
  31. }  
  32. #endif  
wtime.h:
[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. #ifndef _WTIME_  
  2. #define _WTIME_  
  3.   
  4. double wtime();  
  5.   
  6. #endif  


cuda_wtime.cu:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. #include <stdio.h>  
  2. #include <sys/time.h>  
  3. #include <iostream>  
  4. #include <cstdlib>  
  5.   
  6. double wtime(void)  
  7. {  
  8.     double now_time;  
  9.     struct timeval etstart;  
  10.     struct timezone tzp;  
  11.   
  12.     if(gettimeofday(&etstart,&tzp)==-1)  
  13.     {  
  14.         perror("Error:calling gettimeofday() not successfully. ");  
  15.     }  
  16.   
  17.     now_time = ( (double)etstart.tv_sec ) + ((double)etstart.tv_usec) / 1000000.0;  
  18.   
  19.     return now_time;  
  20. }  
  21.   
  22. #if 0  
  23. int main()  
  24. {  
  25.     double time;  
  26.     time = wtime();  
  27.   
  28.     printf("time of day = %10.4f ",time);  
  29.   
  30.     return 0;  
  31. }  
  32. #endif  

makefile:

[cpp] view plaincopy在CODE上查看代码片派生到我的代码片
  1. cu:  
  2.     nvcc cuda_mat_mul_v1.cu cuda_wtime.cu  
  3.     ./a.out  

结果：



3.计算性能对比：

矩阵大小
 1600*1600
1200*1200
800*800
320*320
串行时间/s
30.9
11.49865
2.597987
0.162311
并行时间
 grid=100/block=16
grid=75/block=16
grid=50/block=16
 grid=20/block=16
kernel执行时间/s
0.0000319
0.0000309944
0.0000309944
0.0000231266
并行计算总时间（分配内存加+数据拷贝+计算）/s
0.70796
0.439213
0.310214
0.237676

可见，在矩阵规模大的时候非常明显的体现出了GPU强大的计算能力。


注明出处：http://blog.csdn.net/lavorange/article/details/41896591


【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/walccott/p/4957569.html