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

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

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

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

mat_mul.cc:

[cpp]  view plain copy

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\n",t);  
91.       
92.     //print();  
93.   
94.     return 0;  
95. }  

wtime.h:

[cpp]  view plain copy

1. #ifndef _WTIME_  
2. #define _WTIME_  
3.   
4. double wtime();  
5.   
6. #endif  

wtime.cc:

[cpp]  view plain copy

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.\n");  
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\n",time);  
29.   
30.     return 0;  
31. }  
32. #endif  

makefile:

[cpp]  view plain copy

1. target:  
2.     g++ mat_mul.cc wtime.cc  
3.     ./a.out  

结果：

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

代码：

cuda_mat_mul_v1.cu:

[cpp]  view plain copy

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\n",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\n",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 plain copy

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.\n");  
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\n",time);  
29.   
30.     return 0;  
31. }  
32. #endif  

wtime.h:

[cpp]  view plain copy

1. #ifndef _WTIME_  
2. #define _WTIME_  
3.   
4. double wtime();  
5.   
6. #endif  

cuda_wtime.cu:

[cpp]  view plain copy

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.\n");  
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\n",time);  
29.   
30.     return 0;  
31. }  
32. #endif  

makefile:

[cpp]  view plain copy

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