     */
 }
 
-//What we should test:
-//-scan
-//-reduce
-//-histogram
+
+struct ConvertRV32ToGrey{
+    static bool init;
+    static pop::UI8 _look_up_table[256][256][256];
+    static pop::UI8 lumi(const pop::VecN<4,pop::UI8> &rgb){
+        if(init==false){
+            init= true;
+            for(unsigned int i=0;i<256;i++){
+                for(unsigned int j=0;j<256;j++){
+                    for(unsigned int k=0;k<256;k++){
+                        _look_up_table[i][j][k]=pop::ArithmeticsSaturation<pop::UI8,pop::F64>::Range(0.299*i + 0.587*j + 0.114*k+0.000001);
+                    }
+                }
+            }
+        }
+        return _look_up_table[rgb(2)][rgb(1)][rgb(0)];
+    }
+};
+bool ConvertRV32ToGrey::init =false;
+pop::UI8 ConvertRV32ToGrey::_look_up_table[256][256][256];
+
+template<typename T>
+__device__ T arithmeticSaturationRangeGPU(pop::F64 p, T limitMin, T limitMax) {
+    if(p>=limitMax) return limitMax;
+    else if(p<limitMin) return limitMin;
+    else return static_cast<T>(p);
+}
+
+template<typename T>
+__global__ void lumi(T *initImage, T *greyImage, const int sizeGreyImage, const T limitMin, const T limitMax) {
+	// thread tx of block bx: works on values from blockSize*bx + tx*4 to blockSize*bx + tx*4 + 3
+	int idx = blockDim.x*blockIdx.x + threadIdx.x;
+	int start = idx<<2;
+
+	if (idx < sizeGreyImage) {
+		pop::F64 p = 0.000001;
+		p += 0.299*initImage[start+2];
+		p += 0.587*initImage[start+1];
+		p += 0.114*initImage[start];
+		greyImage[idx] = arithmeticSaturationRangeGPU(p, limitMin, limitMax);
+	}
+}
+
+template<typename T>
+void lumi_test() {
+    uint64_t start_time, stop_time, diff_time;
+    cudaEvent_t cuda_start_time, cuda_stop_time;
+    cudaError_t error;
+    float msecTotal;
+
+    std::cout << "---\tWorking with " << type_name<T>() << "\t---" << std::endl;
+
+    pop::MatN<2,pop::VecN<4,T> > h_image_RV32(MATRIX_DIM,MATRIX_DIM); // the initial image
+    for (int i=0; i<MATRIX_DIM; i++) {
+        for (int j=0; j<MATRIX_DIM; j++) {
+        	pop::VecN<4,T> v;
+        	v(0) = i%256;
+        	v(1) = (255-(j%256));
+        	v(2) = (i+j*10)%256;
+        	v(3) = 0;
+        	h_image_RV32(i, j) = v;
+        }
+    }
+
+	pop::MatN<2, T> imggrey_pop(MATRIX_DIM, MATRIX_DIM); // the image after computation on the CPU
+	pop::MatN<2, T> h_imggrey(MATRIX_DIM, MATRIX_DIM); // the image on host side after computation on the GPU
+
+
+    /********************** CPU lumi **********************/
+	rdtsc(start_time);
+	for (int i=0; i<N_ROUNDS; i++) {
+		std::transform(h_image_RV32.begin(),h_image_RV32.end(),imggrey_pop.begin(),ConvertRV32ToGrey::lumi);
+	}
+	rdtsc(stop_time);
+
+	diff_time = diffTime(stop_time, start_time)/(N_ROUNDS*1000.0);
+	std::cout << "CPU lumi: " << diff_time << "ms" << std::endl;
+
+
+    /********************** CUDA initializations **********************/
+
+    error = cudaEventCreate(&cuda_start_time);
+    if (error != cudaSuccess)
+    {
+        std::cerr << "Failed to create start event (error code " << cudaGetErrorString(error) << ")!" << std::endl;
+    }
+    error = cudaEventCreate(&cuda_stop_time);
+    if (error != cudaSuccess)
+    {
+        std::cerr << "Failed to create stop event (error code " << cudaGetErrorString(error) << ")!" << std::endl;
+    }
+
+    T* d_image_RV32;
+    const int sizeInitialImage = h_image_RV32.size() * 4 * sizeof(T);
+    cudaMalloc((void**)&d_image_RV32, sizeInitialImage);
+
+    T* d_imggrey;
+    const int sizeGreyImage = h_imggrey.size() * sizeof(T);
+    cudaMalloc((void**)&d_imggrey, sizeGreyImage);
+
+    /********************** GPU lumi **********************/
+
+    error = cudaEventRecord(cuda_start_time, NULL);
+    if (error != cudaSuccess)
+    {
+    	std::cerr << "Failed to record start event (error code " << cudaGetErrorString(error) << ")!" << std::endl;
+    }
+
+    for (int i=0; i<N_ROUNDS; i++) {
+    	cudaMemcpy(d_image_RV32, h_image_RV32.data(), sizeInitialImage, cudaMemcpyHostToDevice);
+
+    	lumi<<<h_imggrey.size()/MAX_NB_THREADS, MAX_NB_THREADS>>>(d_image_RV32, d_imggrey, h_imggrey.size(), pop::NumericLimits<T>::minimumRange(), pop::NumericLimits<T>::maximumRange());
+
+    	error = cudaPeekAtLastError();
+    	if (error != cudaSuccess) {
+    		std::cerr << "There was an error executing the kernel: " << cudaGetErrorString(error) << std::endl;
+    	}
+
+        cudaMemcpy(h_imggrey.data(), d_imggrey, sizeGreyImage, cudaMemcpyDeviceToHost);
+    }
+
+    error = cudaEventRecord(cuda_stop_time, NULL);
+    if (error != cudaSuccess)
+    {
+    	std::cerr << "Failed to record stop event (error code " << cudaGetErrorString(error) << ")!" << std::endl;
+    }
+
+    error = cudaEventSynchronize(cuda_stop_time);
+    if (error != cudaSuccess)
+    {
+    	std::cerr << "Failed to synchronize on the stop event (error code " << cudaGetErrorString(error) << ")!" << std::endl;
+    }
+
+    msecTotal = 0.0f;
+    cudaEventElapsedTime(&msecTotal, cuda_start_time, cuda_stop_time);
+    std::cout << "GPU lumi: " << msecTotal / N_ROUNDS << "ms" << std::endl;
+
+
+    /********************** Ok, done **********************/
+
+    cudaMemcpy(h_imggrey.data(), d_imggrey, sizeGreyImage, cudaMemcpyDeviceToHost);
+    cudaFree(d_image_RV32);
+    cudaFree(d_imggrey);
+
+    std::cout << "Is the result correct? " << (h_imggrey == imggrey_pop ? "true" : "false") << std::endl << std::endl;
+    /*
+    for (int i=0; i<MATRIX_DIM; i++) {
+    	for (int j=0; j<MATRIX_DIM; j++) {
+    		if (h_cudaResult(i, j) != h_imggrey(i, j)) {
+    			std::cout << "Diff: c(" << i << ", " << j << ") = " << h_cudaResult(i, j) << "\t!=\tc_pop(" << i << ", " << j << ") = " << h_imggrey(i, j) << std::endl;
+    		}
+    	}
+    }
+    */
+}
 
 void cuda_test()
 {
 	init_clock_mhz();
 
-	multiply_matrix_test<pop::UI32>();
-	multiply_matrix_test<pop::F32>();
-	multiply_matrix_test<pop::F64>();
+	// Testing matrix multiply
+	if (false) {
+		multiply_matrix_test<pop::UI32>();
+		multiply_matrix_test<pop::F32>();
+		multiply_matrix_test<pop::F64>();
+	}
+
+	// cost of lumi
+	lumi_test<pop::UI8>();
 }
 #endif

#include <unistd.h>
#include <cstdlib>
#include <cxxabi.h>

const int MATRIX_DIM = 1024;
const int BLOCK_SIZE = 32;


template<typename T>
std::string type_name()

    std::cout << "GPU parallel matrix multiplication (shared): " << msecTotal << "ms" << std::endl;
#endif


    /********************** Ok, done **********************/

    pop::MatN<2, T> h_C(MATRIX_DIM, MATRIX_DIM);
