 //
 // For more details on the reduction algorithm (notably the multi-pass approach), see
 // the "reduction" sample in the CUDA SDK.
-template <unsigned int blockSize, bool nIsPow2>
+template <
+    unsigned int blockSize
+  , bool nIsPow2
+  , typename typ
+  , typename styp
+  , class task
+  >
 __device__  __inline__
-void reduceSinglePass_dev(const cuDoubleComplex *g_idata, double *g_odata, unsigned int n)
+void reduceSinglePass_devGen(const styp *g_idata, typ *g_odata, unsigned int n)
 {
 
     //
     // PHASE 1: Process all inputs assigned to this block
     //
 
-    reduceBlocks<blockSize, nIsPow2>(g_idata, g_odata, n);
+    reduceBlocksGen<blockSize, nIsPow2, typ, styp, task>(g_idata, g_odata, n);
 
     //
     // PHASE 2: Last block finished will process all partial sums

    }
}

#endif // #ifndef _REDUCE_KERNEL_H_

    See the CUDA SDK "reduction" sample for more information.
*/

template <
    unsigned int blockSize
  , typename typ
  , class task
  >
__device__ __inline__ void
reduceBlockGen(volatile typ *sdata, typ myAcc, const unsigned int tid)
{
    sdata[tid] = myAcc;
    __syncthreads();

    // do reduction in shared mem
    if (blockSize >= 512)
    {
        if (tid < 256) sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid + 256]);
        __syncthreads();
    }

    if (blockSize >= 256)
    {
        if (tid < 128) sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid + 128]);
        __syncthreads();
    }

    if (blockSize >= 128)
    {
        if (tid <  64) sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid +  64]);
        __syncthreads();
    }

    if (tid < 32)
    {
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid + 32]);
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid + 16]);
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid +  8]);
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid +  4]);
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid +  2]);
        sdata[tid] = myAcc = task::reduce(myAcc, sdata[tid +  1]);
    }
}

template <
    unsigned int blockSize
  , bool nIsPow2
  , typename typ
  , typename styp
  , class task
  >
__device__ __inline__ void
reduceBlocksGen(const styp *g_idata, typ *g_odata, unsigned int n)
{
    extern __shared__ typ sdata[];

    // perform first level of reduction,
    // reading from global memory, writing to shared memory
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x*(blockSize*2) + threadIdx.x;
    unsigned int gridSize = blockSize*2*gridDim.x;
    typ myAcc = task::init();

    // we reduce multiple elements per thread.  The number is determined by the
    // number of active thread blocks (via gridDim).  More blocks will result
    // in a larger gridSize and therefore fewer elements per thread
    while (i < n)
    {
        myAcc = task::reduce(myAcc, task::f(g_idata[i]));

        // ensure we don't read out of bounds -- this is optimized away for powerOf2 sized arrays
        if (nIsPow2 || i + blockSize < n)
            myAcc = task::reduce(myAcc, task::f(g_idata[i+blockSize]));

        i += gridSize;
    }

    // do reduction in shared mem
    reduceBlockGen<blockSize, typ, task>(sdata, myAcc, tid);

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