reduce

#include <cuda_runtime.h>

__global__ void reduction_kernel(const float* input, float* output, int N) {
    // int idx = blockIdx.x * blockDim.x + threadIdx.x;
    // if (idx >= N) return;
    // atomicAdd(output, input[idx]);
    extern __shared__ double sdata[];
    int tid = threadIdx.x;
    int idx = blockDim.x * blockIdx.x + tid;

    double local = 0;
    if (idx >= N) {
        local = 0;
    } else {
        local = input[idx];
    }

    sdata[tid] = local;

    for (int s=blockDim.x / 2; s > 0; s >>= 1) {
        __syncthreads();

        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
    }

    __syncthreads();
    if (tid == 0) {
        // output = sdata[0];
        atomicAdd(output, (float) sdata[0]);
    }
}

// input, output are device pointers
extern "C" void solve(const float* input, float* output, int N) {
    int threadsPerBlock = 1024;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    cudaMemset(output, 0, sizeof(float));

    reduction_kernel<<<blocksPerGrid, threadsPerBlock, threadsPerBlock * sizeof(double)>>>(input, output, N);
    cudaDeviceSynchronize();


}