 
 void reductionExperiment(int _nr_ele)
 {
+  unsigned int nr_iterations = 10;
   // Generate Test data
   int nr_patches = _nr_ele;
   int nr_elements = nr_patches*kPatchArea;

  }

  // print result
  std::cout << "Hessian Normal" << std::endl;
  float hessian_output_array[kJacobianSize*kJacobianSize];
  for(unsigned int row = 0, index = 0; row < kJacobianSize; ++row)
  {
    for(unsigned int col = row; col < kJacobianSize; ++col,++index)
    {
      hessian_output_array[row*kJacobianSize + col] = hessian_output_array[col*kJacobianSize + row] =  hessian_array_host[index];
    }
  }

  for(unsigned int row = 0; row < kJacobianSize; ++row)
  {
    for(unsigned int col = 0; col < kJacobianSize; ++col)
    {
      std::cout << std::setw( 12 ) << hessian_output_array[row*kJacobianSize + col] << " " ;
    }
    std::cout << std::endl;
  }
  std::cout << std::endl;

  std::cout << "Gradient Normal" << std::endl;
  for(unsigned int ii = 0; ii < kJacobianSize; ++ii)
  {
    std::cout << gradient_array_host[ii] << " " ;
  }
  std::cout << std::endl;


  // ----------------------- reduce brent
  output_cache = 1;
  std::cout << "reducing brent blocks/problem_size: " << blocks_brent << "," << nr_elements << std::endl;
  reduceJacobian(nr_elements,threads,blocks_brent,jacobian_input_device,visibility_input_device,residual_input_device,gradient_output,hessian_output);

  problem_size = blocks_brent;
  while(blocks_brent > 1)
  {
    threads = (problem_size < max_threads*2) ? nextPow2((problem_size + 1)/ 2) : max_threads;
    blocks_brent = (problem_size + (threads * 2 - 1)) / (threads * 2);
    std::cout << "reducing brent blocks/problem_size: " << blocks_brent << "," << problem_size << std::endl;
    reduceHessianGradient(problem_size,threads,blocks_brent,gradient_output,hessian_output,gradient_output_2,hessian_output_2);
    output_cache = 2;
    if(blocks_brent > 1)
    {
      problem_size = blocks_brent;
      threads = (problem_size < max_threads*2) ? nextPow2((problem_size + 1)/ 2) : max_threads;
      blocks_brent = (problem_size + (threads * 2 - 1)) / (threads * 2);
      std::cout << "reducing brent blocks/problem_size: " << blocks_brent << "," << problem_size << std::endl;
      reduceHessianGradient(problem_size,threads,blocks_brent,gradient_output_2,hessian_output_2,gradient_output,hessian_output);
      output_cache = 1;
    }
  }

  //get result
  if(output_cache == 1)

  }

  // print result
  //  std::cout << "Hessian Brent" << std::endl;
  //  for(unsigned int row = 0, index = 0; row < kJacobianSize; ++row)
  //  {
  //    for(unsigned int col = row; col < kJacobianSize; ++col,++index)
  //    {
  //      hessian_output_array[row*kJacobianSize + col] = hessian_output_array[col*kJacobianSize + row] =  hessian_array_host[index];
  //    }
  //  }

  //  for(unsigned int row = 0; row < kJacobianSize; ++row)
  //  {
  //    for(unsigned int col = 0; col < kJacobianSize; ++col)
  //    {
  //      std::cout << std::setw( 12 ) << hessian_output_array[row*kJacobianSize + col] << " " ;
  //    }
  //    std::cout << std::endl;
  //  }
  //  std::cout << std::endl;

  //  std::cout << "Gradient Brent" << std::endl;
  //  for(unsigned int ii = 0; ii < kJacobianSize; ++ii)
  //  {
  //    std::cout << gradient_array_host[ii] << " " ;
  //  }
  //  std::cout << std::endl;

  std::clock_t c_start_cpu = std::clock();
  for(unsigned int tt = 0; tt < nr_iterations;++tt)
  {
    Eigen::Matrix<float,kJacobianSize,kJacobianSize> H;
    Eigen::Matrix<float,kJacobianSize,1> g;

    reduceJacobianCPU(nr_elements,
                      jacobian_input_host,
                      visibility_input_host,
                      residual_input_host,
                      H,
                      g);
  }
  std::clock_t c_end_cpu = std::clock();
  double time_cpu= CLOCK_TO_MS(c_start_cpu,c_end_cpu);

  //  std::cout << "Hessian CPU" << std::endl;
  //  std::cout << H << std::endl;
  //  std::cout << "Gradient CPU" << std::endl;
  //  std::cout << g.transpose() << std::endl;

  std::cout << "Time normal " << time_gpu_normal << std::endl;
  std::cout << "Time brent " << time_gpu_brent << std::endl;
  std::cout << "Time CPU "  << time_cpu << std::endl;

  free(jacobian_input_host);
  free(visibility_input_host);
