【CUDA并行编程之三】Cuda矢量求和运算

本文将通过矢量求和运算来说明基本的Cuda并行编程的基本概念。所谓矢量求和运算，就是两个数组数据中对应的元素两两相加，并将结果保存在第三个数组中。如下图所示：

1.基于CPU的矢量求和：

代码非常简单：

[cpp] view plaincopy

1. #include<iostream>  
2.   
3. using namespace std;  
4.   
5. const int N =10;  
6.   
7. void add( int *a ,int *b , int *c)  
8. {  
9.     int tid = 0;  
10.     while(tid < N)       
11.     {  
12.         c[tid] = a[tid] + b[tid];  
13.         tid += 1;            
14.     }  
15. }  
16.   
17. int main()  
18. {  
19.     int a[N],b[N],c[N];  
20.     //在CPU上对数组'a'和'b'赋值  
21.     for(int i=0;i<N;i++)     
22.     {  
23.         a[i] = -1;  
24.         b[i] = i * i;          
25.     }  
26.     add(a,b,c);  
27.     //打印结果  
28.     for(int i=0;i<N;i++)   
29.     {  
30.         cout<<a[i]<<" + "<<b[i]<<" = "<<c[i]<<endl;          
31.     }  
32.     return 0;  
33. }  

上面采用while循环虽然有些复杂，但这是为了使得代码能够在拥有多个CPU或者CPU核的系统上并行运行。例如，在双核处理器上可以将每次递增的大小改为2，这样其中一个核从tid=0开始执行循环，而另一个核从tid=1开始执行循环。第一个核将偶数索引的元素相加，而第二个核则将奇数索引的元素相加。这相当于在每个CPU核上执行以下代码：

[cpp] view plaincopy

1. <strong>一个CPU核： </strong>  
2. void add( int *a ,int *b , int *c)  
3. {  
4.    <strong> int tid = 0;</strong>  
5.     while(tid < N)       
6.     {  
7.         c[tid] = a[tid] + b[tid];  
8.         <strong>tid += 2;  </strong>          
9.     }  
10. }  
11.   
12. <strong>第2个CPU核：</strong>  
13. void add( int *a ,int *b , int *c)  
14. {  
15.     <strong>int tid = 1;</strong>  
16.     while(tid < N)       
17.     {  
18.         c[tid] = a[tid] + b[tid];  
19.         <strong>tid += 2;</strong>            
20.     }  
21. }  

当然，要在CPU实际执行这个运算，还需要增加更多的代码。例如，需要编写一定数量的代码来创建工作线程，每个线程都执行函数add()，并假设每个线程都将并行执行。然而，这种假设是一种理想但不实际的想法，线程调度机制的实际运行情况往往并非如此。

2.基于GPU的矢量求和：

我们可以在GPU上实现相同的加法运算，这需要将add()编写为一个设备函数。先把代码呈上：

[cpp] view plaincopy

1. #include<iostream>  
2.   
3. using namespace std;  
4.   
5. #define N 10  
6.   
7. <strong>__global__</strong> void add(int *a , int *b , int *c)  
8. {  
9.     int tid = blockIdx.x;  
10.     if(tid<N)  
11.     {  
12.         c[tid] = a[tid]+b[tid];  
13.     }  
14. }  
15.   
16. int main( void )  
17. {  
18.     int a[N],b[N],c[N];  
19.     int *dev_a , *dev_b, *dev_c;  
20.   
21.     //allocate memory on GPU  
22.     cudaMalloc( (void**)&dev_a, N*sizeof(int) ) ;  
23.     cudaMalloc( (void**)&dev_b, N*sizeof(int) ) ;  
24.     cudaMalloc( (void**)&dev_c, N*sizeof(int) ) ;  
25.   
26.     for(int i=0;i<N;i++)  
27.     {  
28.         a[i] = -1;  
29.         b[i] = i * i ;  
30.     }  
31.   
32.     cudaMemcpy( dev_a , a, N*sizeof(int), cudaMemcpyHostToDevice ) ;  
33.     cudaMemcpy( dev_b , b, N*sizeof(int), cudaMemcpyHostToDevice ) ;  
34.   
35.     <strong>add<<<N,1>>>(dev_a,dev_b,dev_c);</strong>  
36.   
37.     cudaMemcpy( c , dev_c , N*sizeof(int), cudaMemcpyDeviceToHost) ;  
38.   
39.     for(int i=0;i<N;i++)  
40.     {  
41.         cout<<a[i]<<"+"<<b[i]<<"="<<c[i]<<endl;  
42.     }  
43.   
44.     //release the memory on GPU  
45.     cudaFree(dev_a);  
46.     cudaFree(dev_b);  
47.     cudaFree(dev_c);  
48.   
49.     return 0;  
50. }  

运行结果：

解释一下代码：

+cudaMalloc()：在设备上三个数组分配内存，其中dev_a，dev_b中包含了输入值，而在数组dev_c中包含了计算结果。

+cudaFree()：避免内存泄露，在使用完GPU内存后通过cudaFree()释放它们。

+cudaMemcpy()：将输入数据复制到设备中，同时制定参数cudaMemcpyHostToDevice，在计算完成后，将计算结果通过参数cudaMemcpyDeviceToHost复制回主机。

+通过尖括号语法，在主机代码main()中执行add()中的设备代码。

__global__：为了函数add()能够在设备上执行，在函数名前面添加了修饰符__global__。

核函数：kernel<<<1,1>>>(param1,param2,...);

但是在这个示例中，尖括号中的数值并不是1：add<<<N,1>>>( dev_a, dev_b ,dev_c );

核函数中的第一个参数：number of blocks.即块的个数。

核函数中的第二个参数：thread per block.即每个线程块中线程的个数。

例如，如果指定了kernel<<<2,1>>>，那么可以认为运行时将创建核函数的两个副本，并以并行的方式来运行它们。我们将每个执行环境都成为一个线程块（block）。如果指定的kernel<<<256,1>>>()，那么将有256个线程块在GPU上运行。

3.用vector动态分配数组。

代码：

[cpp] view plaincopy

1. #include<iostream>  
2. #include<vector>  
3.   
4. using namespace std;  
5.   
6. const int N = 10;  
7.   
8. __global__ void add(int* a , int* b ,int* c)  
9. {  
10.     int tid = blockIdx.x;  
11.     if(tid<N)  
12.     {  
13.         c[tid]  = a[tid] + b[tid];  
14.     }  
15. }  
16.   
17. int main()  
18. {  
19.     vector<int> vec_a,vec_b;  
20.     int *va,*vb,*vc;  
21.     int *dev_a,*dev_b,*dev_c;  
22.   
23.     cudaMalloc( (void**)&dev_a,N*sizeof(int) ) ;  
24.     cudaMalloc( (void**)&dev_b,N*sizeof(int) ) ;  
25.     cudaMalloc( (void**)&dev_c,N*sizeof(int) ) ;  
26.   
27.     for(int i=0;i<N;i++)  
28.     {  
29.         vec_a.push_back(-1);//vec_a[i] = -1;  
30.         vec_b.push_back(i*i);//vec_b[i] = i * i;  
31.     }  
32.       
33.     <strong>/* 
34.      * 第一种方式 
35.      */  
36.     va = new int[N];  
37.     vb = new int[N];  
38.     copy(vec_a.begin(),vec_a.end(),va);  
39.     copy(vec_b.begin(),vec_b.end(),vb);  
40.     /* 
41.      * 第二种方式 
42.     va = (int *)&vec_a[0];//vector to array 
43.     vb = (int *)&vec_b[0]; 
44.     */</strong>  
45.   
46.     cudaMemcpy(dev_a,va,N*sizeof(int),cudaMemcpyHostToDevice) ;  
47.     cudaMemcpy(dev_b,vb,N*sizeof(int),cudaMemcpyHostToDevice) ;  
48.   
49.     add<<<N,1>>>(dev_a,dev_b,dev_c);  
50.   
51.     vc = new int[N];  
52.     cudaMemcpy(vc,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost) ;  
53.   
54. #if 1  
55.     for(int i=0;i<N;i++)  
56.     {  
57.         cout<<va[i]<<"+"<<vb[i]<<"="<<vc[i]<<endl;  
58.     }  
59. #endif  
60.     cudaFree(dev_a);  
61.     cudaFree(dev_b);  
62.     cudaFree(dev_c);  
63.   
64.     return 0;  
65. }  

在这里主要还是讨论一下从vector转换成为数组array的问题。

因为对于vector来将，它在内存中的存储一定是连续的，那么按照如下方式写就非常简单而且没有问题：

[cpp] view plaincopy

1. std::vector<double> v;  
2. double* a = &v[0];  

而如果对于在内存中存储不连续的话，那么就要用令一种方法，copy：

[cpp] view plaincopy

1. double arr[100];  
2. std::copy(v.begin(), v.end(), arr);  

有三个链接讨论这个问题：

1.http://stackoverflow.com/questions/2923272/how-to-convert-vector-to-array-c?answertab=active#tab-top

2.http://www.cplusplus.com/forum/beginner/7477/

3.http://www.cplusplus.com/reference/algorithm/copy/

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