 			cudaStreamDestroy(stream[grid]);
 			cudaFree(output[grid]);
 		} // for
-		mytimer2.stop();
-		//std::cout << "Housekeeping time: " << mytimer2.elapsed_msec() << " ms." << std::endl;
 
 		return result;
 	} // reduce()
 
 
+	/* the main reduce - performs better than thrust on double, worse on single *
+	 * it performs much better than thrust in RMC code */
 
 	template <typename data_t, typename functor_t>
 	__global__

	} // 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];


			while(1) {
				unsigned int num_blocks = ceil((float) n_grid / (block_size * num_subtiles));
				reduce_block <<< num_blocks, block_size, d_shmem_size >>>
								(input, n_grid, init, op, output, num_subtiles);
				cudaError_t err = cudaGetLastError();
