
#ifndef _SCAN_WORKEFFICIENT_KERNEL_H_
#define _SCAN_WORKEFFICIENT_KERNEL_H_

#include "my_cutil.h"

////////////////////////////////////////////////////////////////////////////////
//   Notes for running in PUG:
//   Bitvector size: 14 bits
//   Need to turn on automatic loop refinement
////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
//! Work-efficient compute implementation of scan, one thread per 2 elements
//! Work-efficient: O(log(n)) steps, and O(n) adds.
//! Also shared storage efficient: Uses n elements in shared mem -- no ping-ponging
//! Uses a balanced tree type algorithm.  See Blelloch, 1990 "Prefix Sums 
//! and Their Applications", or Prins and Chatterjee PRAM course notes:
//! http://www.cs.unc.edu/~prins/Classes/203/Handouts/pram.pdf
//!
//! Pro: Work Efficient
//! Con: Shared memory bank conflicts due to the addressing used.
//
//! @param g_odata  output data in global memory
//! @param g_idata  input data in global memory
//! @param n        input number of elements to scan from input data
///////////////////////////////////////////////////////////////////////////////

__global__ void scan_workefficient_kernel (float *g_odata, float *g_idata, int n)
{
  assume(n >= blockDim.x);
  assume(blockDim.x <= 1 << 13);    // here the n in 1 << n is the number of bits minus 1
                                    // this assumption is for preventing bitvector overflowing
                                    // since the 2 * x is converted to x << 1

    // Dynamically allocated shared memory for scan kernels
    extern  __shared__  float temp[];

    int thid = threadIdx.x;

    int offset = 1;

    // Cache the computational window in shared memory
    temp[2*thid]   = g_idata[2*thid];
    temp[2*thid+1] = g_idata[2*thid+1];

    // build the sum in place up the tree
    for (int d = n>>1; d > 0; d >>= 1)
    {
      __syncthreads();

      offset *= 2;

      if (thid < d)
        {
	  int ai = offset/2*(2*thid+1)-1;
	  int bi = offset/2*(2*thid+2)-1;
	  
	  temp[bi] = 1;
	  temp[bi] += temp[ai];
        }
    }

    // scan back down the tree

    // clear the last element
    if (thid == 0)
    {
        temp[n - 1] = 0;
    }

    // traverse down the tree building the scan in place
    for (int d = 1; d < n; d *= 2)
    {
      __syncthreads();

      offset >>= 1;
      
      if (thid < d)
        {
	  int ai = offset*(2*thid+1)-1;
	  int bi = offset*(2*thid+2)-1;
	  
	  float t = temp[ai];
	  temp[ai]  = temp[bi];
	  temp[bi] += t;
        }
    }

#ifndef MUTATION
    __syncthreads();
#endif

    // write results to global memory
    g_odata[2*thid]   = temp[2*thid];
    g_odata[2*thid+1] = temp[2*thid+1];

}

#endif // #ifndef _SCAN_WORKEFFICIENT_KERNEL_H_
