#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[];

/* This version uses contiguous threads, but its interleaved 
   addressing results in many shared memory bank conflicts. */

__global__ void reduce1kernel(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) 
    {
        int index = 2 * s * tid;
        if (index < blockDim.x) 
        {
            sdata[index] += sdata[index + s];
        }
	      __syncthreads();
    }

    // write result for this block to global mem
     /* BUGINJECT: ADD_BARRIER, DOWN */
    if (tid == 0) {
#ifdef MUTATION
      __syncthreads();
#endif
      g_odata[blockIdx.x] = sdata[0];
    }
}
