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

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

1.基于CPU的矢量求和：

代码非常简单：

#include<iostream>

using namespace std;

const int N =10;

void add( int *a ,int *b , int *c)
{
    int tid = 0;
    while(tid < N)
    {
        c[tid] = a[tid] + b[tid];
        tid += 1;
    }
}

int main()
{
    int a[N],b[N],c[N];
    //在CPU上对数组'a'和'b'赋值
    for(int i=0;i<N;i++)
    {
        a[i] = -1;
        b[i] = i * i;
    }
    add(a,b,c);
    //打印结果
    for(int i=0;i<N;i++)
    {
        cout<<a[i]<<" + "<<b[i]<<" = "<<c[i]<<endl;
    }
    return 0;
}

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

<strong>一个CPU核： </strong>
void add( int *a ,int *b , int *c)
{
   <strong> int tid = 0;</strong>
    while(tid < N)
    {
        c[tid] = a[tid] + b[tid];
        <strong>tid += 2;  </strong>
    }
}

<strong>第2个CPU核：</strong>
void add( int *a ,int *b , int *c)
{
    <strong>int tid = 1;</strong>
    while(tid < N)
    {
        c[tid] = a[tid] + b[tid];
        <strong>tid += 2;</strong>
    }
}

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

2.基于GPU的矢量求和：

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

#include<iostream>

using namespace std;

#define N 10

<strong>__global__</strong> void add(int *a , int *b , int *c)
{
    int tid = blockIdx.x;
    if(tid<N)
    {
        c[tid] = a[tid]+b[tid];
    }
}

int main( void )
{
    int a[N],b[N],c[N];
    int *dev_a , *dev_b, *dev_c;

    //allocate memory on GPU
    cudaMalloc( (void**)&dev_a, N*sizeof(int) ) ;
    cudaMalloc( (void**)&dev_b, N*sizeof(int) ) ;
    cudaMalloc( (void**)&dev_c, N*sizeof(int) ) ;

    for(int i=0;i<N;i++)
    {
        a[i] = -1;
        b[i] = i * i ;
    }

    cudaMemcpy( dev_a , a, N*sizeof(int), cudaMemcpyHostToDevice ) ;
    cudaMemcpy( dev_b , b, N*sizeof(int), cudaMemcpyHostToDevice ) ;

    <strong>add<<<N,1>>>(dev_a,dev_b,dev_c);</strong>

    cudaMemcpy( c , dev_c , N*sizeof(int), cudaMemcpyDeviceToHost) ;

    for(int i=0;i<N;i++)
    {
        cout<<a[i]<<"+"<<b[i]<<"="<<c[i]<<endl;
    }

    //release the memory on GPU
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}

运行结果：

解释一下代码：

+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动态分配数组。

代码：

#include<iostream>
#include<vector>

using namespace std;

const int N = 10;

__global__ void add(int* a , int* b ,int* c)
{
    int tid = blockIdx.x;
    if(tid<N)
    {
        c[tid]  = a[tid] + b[tid];
    }
}

int main()
{
    vector<int> vec_a,vec_b;
    int *va,*vb,*vc;
    int *dev_a,*dev_b,*dev_c;

    cudaMalloc( (void**)&dev_a,N*sizeof(int) ) ;
    cudaMalloc( (void**)&dev_b,N*sizeof(int) ) ;
    cudaMalloc( (void**)&dev_c,N*sizeof(int) ) ;

    for(int i=0;i<N;i++)
    {
        vec_a.push_back(-1);//vec_a[i] = -1;
        vec_b.push_back(i*i);//vec_b[i] = i * i;
    }

    <strong>/*
     * 第一种方式
     */
    va = new int[N];
    vb = new int[N];
    copy(vec_a.begin(),vec_a.end(),va);
    copy(vec_b.begin(),vec_b.end(),vb);
    /*
     * 第二种方式
    va = (int *)&vec_a[0];//vector to array
    vb = (int *)&vec_b[0];
    */</strong>

    cudaMemcpy(dev_a,va,N*sizeof(int),cudaMemcpyHostToDevice) ;
    cudaMemcpy(dev_b,vb,N*sizeof(int),cudaMemcpyHostToDevice) ;

    add<<<N,1>>>(dev_a,dev_b,dev_c);

    vc = new int[N];
    cudaMemcpy(vc,dev_c,N*sizeof(int),cudaMemcpyDeviceToHost) ;

#if 1
    for(int i=0;i<N;i++)
    {
        cout<<va[i]<<"+"<<vb[i]<<"="<<vc[i]<<endl;
    }
#endif
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}

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

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

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

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

double arr[100];
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