     return a + b;
   }
 };
-
 
 template <typename T, typename Op, unsigned int block_size>
 __global__ void reduce(const T *a, unsigned int n, T *b) {

    // Reduce in shared memory
    sdata[tid] = reduced;
    __syncthreads();
    if (block_size >= 512) {
      if (tid < 256) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 256]);
      }
      __syncthreads();
    }
    if (block_size >= 256) {
      if (tid < 128) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 128]);
      }
      __syncthreads();
    }
    if (block_size >= 128) {
      if (tid <  64) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 64]);
      }
      __syncthreads();
    }
    // No need to sync threads in the same warp
    if (tid < 32) {
      if (block_size >= 64) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 32]);
      }
      if (block_size >= 32) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 16]);
      }
      if (block_size >= 16) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 8]);
      }
      if (block_size >= 8) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 4]);
      }
      if (block_size >= 4) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 2]);
      }
      if (block_size >= 2) {
        sdata[tid] = reduced = op(reduced, sdata[tid + 1]);
      }
    }


  x |= x >> 16;
  return ++x;
}

unsigned int n_reduce_blocks(unsigned int n) {
  return min(max_blocks, (n + (2*reduce_cta_size - 1)) / (2*reduce_cta_size));
