 			data_t temp_result;
 			cudaMemcpy(&temp_result, &output[0], sizeof(data_t), cudaMemcpyDeviceToHost);
 			result = op(result, temp_result);
-			//std::cout << grid << ". result: " << result << std::endl;
 		} // for
 
 		output = base_output;

	} // reduce_block()


	// this performs better than thrust::reduce
	template <typename iterator_t, typename data_t, typename functor_t>
	data_t reduce_test(iterator_t start, iterator_t end, data_t init, functor_t op) {

		const unsigned int num_subtiles = NUM_SUBTILES_;
		const unsigned int block_size = BLOCK_DIM_;
		const unsigned int tile_size = block_size * num_subtiles;
		const unsigned int n_grid_max = GRID_SIZE_ * tile_size;
		const unsigned int n = end - start;
		const unsigned int num_grids = ceil((float) n / n_grid_max);

		unsigned int d_shmem_size = block_size * sizeof(data_t);
		unsigned int output_size = GRID_SIZE_;

		data_t *output[num_grids];
		data_t *base_output[num_grids];

namespace cuda {

	const unsigned int MAX_CUDA_THREADS_ = 1024;
	const unsigned int MAX_GRID_SIZE_ = 65535;	// max number of blocks
	const unsigned int BLOCK_DIM_REDUCE_TEST_ = 256;	// number of threads - assuming power of 2
	const unsigned int BLOCK_DIM_REDUCE_ = 256;	// number of threads - assuming power of 2
	const unsigned int BLOCK_DIM_REDUCE_ONE_ = 256;	// number of threads - assuming power of 2
	const unsigned int NUM_SUBTILES_ = 8;	// number of subtiles in a block (processed by one thread block)

	extern __shared__ unsigned char d_data_sh[];

	/* a test imlementation - performs OK */

	template <typename data_t, typename functor_t>
	__global__
