CUDA编程（1）

/** Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.** NVIDIA Corporation and its licensors retain all intellectual property and* proprietary rights in and to this software and related documentation.* Any use, reproduction, disclosure, or distribution of this software* and related documentation without an express license agreement from* NVIDIA Corporation is strictly prohibited.** Please refer to the applicable NVIDIA end user license agreement (EULA)* associated with this source code for terms and conditions that govern* your use of this NVIDIA software.**/#include "book.h"#define imin(a,b) (a<b?a:b)const int N = 33 * 1024;const int threadsPerBlock = 256;const int blocksPerGrid =imin(32, (N + threadsPerBlock - 1) / threadsPerBlock);__global__ void dot(float *a, float *b, float *c) {    //对于每个block cache都产生一个副本    __shared__ float cache[threadsPerBlock];    int tid = threadIdx.x + blockIdx.x * blockDim.x;    int cacheIndex = threadIdx.x;    float   temp = 0;    while (tid < N) {        temp += a[tid] * b[tid];        tid += blockDim.x * gridDim.x;    }    // set the cache values    cache[cacheIndex] = temp;    // 如果一个线程快的线程没有执行完syncthreads前的语句的话，不执行后面的语句    __syncthreads();    //块内线程释放    // for reductions, threadsPerBlock must be a power of 2    // because of the following code    int i = blockDim.x / 2;    while (i != 0) {        if (cacheIndex < i)            cache[cacheIndex] += cache[cacheIndex + i];        //确保前面线程已经OK        __syncthreads();        i /= 2;    }    if (cacheIndex == 0)        c[blockIdx.x] = cache[0];}int main(void) {    float   *a, *b, c, *partial_c;    float   *dev_a, *dev_b, *dev_partial_c;    // allocate memory on the cpu side    a = (float*)malloc(N*sizeof(float));    b = (float*)malloc(N*sizeof(float));    partial_c = (float*)malloc(blocksPerGrid*sizeof(float));    // allocate the memory on the GPU    HANDLE_ERROR(cudaMalloc((void**)&dev_a,        N*sizeof(float)));    HANDLE_ERROR(cudaMalloc((void**)&dev_b,        N*sizeof(float)));    HANDLE_ERROR(cudaMalloc((void**)&dev_partial_c,        blocksPerGrid*sizeof(float)));    // fill in the host memory with data    for (int i = 0; i<N; i++) {        a[i] = i;        b[i] = i * 2;    }    // copy the arrays 'a' and 'b' to the GPU    HANDLE_ERROR(cudaMemcpy(dev_a, a, N*sizeof(float),        cudaMemcpyHostToDevice));    HANDLE_ERROR(cudaMemcpy(dev_b, b, N*sizeof(float),        cudaMemcpyHostToDevice));    dot << <blocksPerGrid, threadsPerBlock >> >(dev_a, dev_b,        dev_partial_c);    // copy the array 'c' back from the GPU to the CPU    HANDLE_ERROR(cudaMemcpy(partial_c, dev_partial_c,        blocksPerGrid*sizeof(float),        cudaMemcpyDeviceToHost));    // finish up on the CPU side    c = 0;    for (int i = 0; i<blocksPerGrid; i++) {        c += partial_c[i];    }#define sum_squares(x)  (x*(x+1)*(2*x+1)/6)    printf("Does GPU value %.6g = %.6g?\n", c,        2 * sum_squares((float)(N - 1)));    // free memory on the gpu side    HANDLE_ERROR(cudaFree(dev_a));    HANDLE_ERROR(cudaFree(dev_b));    HANDLE_ERROR(cudaFree(dev_partial_c));    // free memory on the cpu side    free(a);    free(b);    free(partial_c);    getchar();}