 
 typedef Intrepid::RealSpaceTools<double> rst;
 
+#define BLOCK_SIZE 64;
+
 //Pre-C++11 timing (thanks jeff)
 double getElapsedTime(const timespec start, const timespec end) {
 	timespec temp;

	double elapsedTime_kokkos_noCopy = 0;

	for(int i = 0; i < 5; ++i){
		contractFieldFieldTensorKokkos1D<Kokkos::OpenMP, omp_input_view_t, omp_output_view_t, omp_input_host_t, omp_output_host_t> 
			(omp_hostOut, omp_hostLeft, omp_hostRight, omp_kokkosOut, omp_kokkosLeft,
			 omp_kokkosRight, c, l, r, p, d1, d2, &elapsedTime_kokkos_noCopy);
	}


	std::cout << "cuda kokkos no-copy speedup of " << elapsedTime_serial/elapsedTime_kokkos_noCopyCuda << std::endl;

  Kokkos::finalize();
  
  //Now try the cuda version, start with warmup
  cudaDoContractFieldFieldTensor(cudaOut,cudaLeft,cudaRight, c, l, r,
   p, d1, d2);

  clock_gettime(CLOCK_MONOTONIC, &tic);
  for(int i = 0; i < 5; ++i){
    cudaDoContractFieldFieldTensor(cudaOut,cudaLeft,cudaRight, c, l, r,
     p, d1, d2);
  }

  clock_gettime(CLOCK_MONOTONIC, &toc);
  const double elapsedTime_cuda = getElapsedTime(tic, toc);

  for (int cl = 0; cl < c; ++cl) {
    for (int lbf = 0; lbf < l; ++lbf) {
      for (int rbf = 0; rbf < r; ++rbf) {
        out1_c_l_r(cl, lbf, rbf) = cudaOut[lbf*r*c + rbf*c + cl];
      }
    }
  }

  rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
  if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), Intrepid::NORM_ONE) > zero) {
    std::cout << "\n\nINCORRECT contractFieldFieldTensor (1): check cuda; "
      << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), Intrepid::NORM_ONE) << "\n\n";
  }

  std::cout << "cuda speedup of " << elapsedTime_serial/elapsedTime_cuda << std::endl;


	return 0;
}
