#define __1D_GRID
#define __1D_THREAD_BLOCK
#include <cuda.h>

#define N 1024
__axiom(blockDim.x == N);

/*
    Parallel reduction kernels
*/

/*
    Parallel sum reduction using shared memory
    - takes log(n) steps for n input elements
    - uses n threads
    - only works for power-of-2 arrays
*/

//REVISIT: cannot refer to this array in invariants
//extern __shared__ int sdata[N];

/* This reduction interleaves which threads are active by using the modulo
   operator.  This operator is very expensive on GPUs, and the interleaved 
   inactivity means that no whole warps are active, which is also very 
   inefficient */

__global__ void reduce0kernel(int *g_idata, int *g_odata, unsigned int n)
{
    __requires(__is_pow2(blockDim.x));

    // load shared mem
    __shared__ int sdata[N];
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;

    sdata[tid] = (i < n) ? g_idata[i] : 0;

    __syncthreads();

    // do reduction in shared mem
    for(unsigned int s=1;
#ifndef INFERENCE
        __invariant(s >= 0),
        __invariant(s ==  1  | s ==   2 | s ==    4 | s ==    8 |
                    s == 16  | s ==  32 | s ==   64 | s ==  128 |
                    s == 256 | s == 512 | s == 1024 | s == 2048),
        __invariant(__no_read(g_odata)),
        __invariant(__no_write(g_odata)),
        __invariant(__no_read(sdata)),
        __invariant(__no_write(sdata)),
#endif
                s < blockDim.x; s *= 2) {
        // modulo arithmetic is slow!
          /* BUGINJECT: ADD_BARRIER, DOWN */
        if ((tid % (2*s)) == 0) 
        {
#ifdef MUTATION
            __syncthreads();
#endif
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
