Cuda计算点积-蒲公英云

Cuda计算点积

调试代码的过程中，发现了这样一个问题

#include <iostream>
using namespace std;
int main()
{ 
    int test[10];
    cout << sizeof(test) << endl;
    float *a;
    //a = new float [10];
    cout << sizeof(a) << endl;
    int *b;
    b = new int (10);
    cout << sizeof(b) << endl;
}

本以为直接指定数组大小和先生成指针再指定空间大小的sizeof()会是一样大的，结果不是一样大。

因为这个错误，在cuda c中调试了很久，最后终于发现是这里错了。

cuda c，使用归约法求点积

有下面几个关键点：

1、指定block和thread的数量，不宜过少也不宜过多。

2、使用共享内存，一个线程块中每个线程都可以访问此共享内存

3、使用归约法，一个线程块中对线程求和的部分复杂度节约到了logn，不过线程数量应该是2的k次方，k为整数。

4、使用同步线程__syncthreads()

cuda c求一维数组点积源代码

#include <iostream>
#define imin(a, b) (a<b?a:b)
#define sum_squares(x) (x*(x+1)*(2*x+1)/6) // 平方和
using namespace std;
const int N = 33 * 1024;
const int threadsPerBlocks = 256;
const int blocksPerGrid = imin(32, (N + threadsPerBlocks - 1) / threadsPerBlocks); // 向上取整获取线程块的数量
__global__ void dot(float *a, float *b, float *c) {  // 计算两个一维数组的点积
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    __shared__ float cache[threadsPerBlocks]; // 共享意味着在一个线程块里线程都可以访问这个数组
    int cacheIndex = threadIdx.x;
    float temp = 0;
    while (tid < N) { 
        temp += a[tid] * b[tid];
        tid += gridDim.x * blockDim.x;
    }
    cache[cacheIndex] = temp; // 将所有的乘法结果聚集到一个grid（网格）里面
    __syncthreads(); // 直到上面的所有语句在所有线程里都执行结束才可以进行下面的语句，即为同步
    int i = threadsPerBlocks / 2;
    while (i != 0) {  // 归约法，logn
        if (i > cacheIndex)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }
    if (cacheIndex == 0)
        c[blockIdx.x] = cache[0];
}
int main() { 
    float *a, *b, *partial_c, c;
    float *dev_a, *dev_b, *dev_partial_c;
    a = new float[N];
    b = new float[N];
    partial_c = new float[blocksPerGrid];
    cudaMalloc((void**)&dev_a, N * sizeof(float));
    cudaMalloc((void**)&dev_b, N * sizeof(float));
    cudaMalloc((void**)&dev_partial_c, blocksPerGrid * sizeof(float));
    for (int i = 0; i < N; ++i) { 
        a[i] = i;
        b[i] = 2 * i;
    }
    cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice);
    dot << <blocksPerGrid, threadsPerBlocks >> > (dev_a, dev_b, dev_partial_c);
    cudaMemcpy(partial_c, dev_partial_c, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);
    c = 0.0;
    for (int i = 0; i < blocksPerGrid; ++i) { 
        c += partial_c[i];
    }
    cout << c << " =? " << 2*(float)sum_squares(float(N - 1)) << endl;
}