Cuda Sum-Reduce in Shared Memory

#include "../common/book.h"
#define N   (1024*2048)
#define NumSMPs 14    #Machine dependent
#define ThreadsPerBlock   1024
#define BlocksPerGrid     N/ThreadsPerBlock



__global__ void add_reduce( float *a, float *b, float *c ) {
// Shared memory for results of multiplication

__shared__ float cache[ThreadsPerBlock];
int gindex = threadIdx.x + blockIdx.x * ThreadsPerBlock;
int lindex = threadIdx.x;
cache[lindex] = a[gindex]*b[gindex];

int offset = ThreadsPerBlock/2;
while (offset != 0) {
     if (lindex < offset){
         cache[lindex] += cache[lindex + offset];}
    __syncthreads();
    offset /= 2; }
if (lindex == 0)
c[blockIdx.x] = cache[0];
}


int main( void ) {
    float  *a, *b, *c;
    float  *dev_a, *dev_b, *dev_c;

    // allocate the memory on the CPU
    a = (float*)malloc( N * sizeof(float) );
    b = (float*)malloc( N * sizeof(float) );
    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_c, BlocksPerGrid * sizeof(float) ) );

    // fill the arrays 'a' and 'b' on the CPU
    for (int i=0; i<N; i++) {
        a[i] = 1.0/(i+1);
        b[i] = 2.0/(i+1);
    }

    // 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 ) );

    add_reduce<<<BlocksPerGrid,ThreadsPerBlock>>>( dev_a, dev_b, dev_c );

    // copy the array 'c' back from the GPU to the CPU
    HANDLE_ERROR( cudaMemcpy( c, dev_c, BlocksPerGrid * sizeof(float),
                              cudaMemcpyDeviceToHost ) );

    // verify that the GPU did the work we requested

    printf("\nResults of GPU\n");
    float sum = 0.0;
     for (int i=0 ; i < BlocksPerGrid; ++i) {
         sum += c[i];}
      printf("Sum of %d Block results: %f\n", BlocksPerGrid, sum);


    // free the memory we allocated on the GPU
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_b ) );
    HANDLE_ERROR( cudaFree( dev_c ) );

    // free the memory we allocated on the CPU
    free( a );
    free( b );
    free( c );

    return 0;
}