   cudaMalloc((void **)& jacobian_input_device,jacobian_input_size*sizeof(float));
   cudaMalloc((void **)& visibility_input_device,visibility_input_size*sizeof(char));
   cudaMalloc((void **)& residual_input_device,residual_input_size*sizeof(float));
+  cudaMemcpy(jacobian_input_device,jacobian_input_host,jacobian_input_size*sizeof(float),cudaMemcpyHostToDevice);
+  cudaMemcpy(visibility_input_device,visibility_input_host,visibility_input_size*sizeof(char),cudaMemcpyHostToDevice);
+  cudaMemcpy(residual_input_device,residual_input_host,residual_input_size*sizeof(float),cudaMemcpyHostToDevice);
+
+  int max_threads = 64;
+
+  int threads = (nr_elements < max_threads*2) ? nextPow2((nr_elements + 1)/ 2) : max_threads;
+  int blocks = (nr_elements + (threads * 2 - 1)) / (threads * 2);
+
+  //Blocks such that each thread sums log(n) elements
+  int nr_ele_per_thread = std::floor(log2 (static_cast<double>(nr_elements)));
+  int blocks_brent = (nr_elements + (threads*nr_ele_per_thread - 1)) / (threads*nr_ele_per_thread);
+
   cudaMalloc((void **)& gradient_output,std::max(blocks,blocks_brent)*kJacobianSize*sizeof(float));
   cudaMalloc((void **)& hessian_output,std::max(blocks,blocks_brent)*kHessianTriagN*sizeof(float));
   cudaMalloc((void **)& gradient_output_2,std::max(blocks,blocks_brent)*kJacobianSize*sizeof(float));
   cudaMalloc((void **)& hessian_output_2,std::max(blocks,blocks_brent)*kHessianTriagN*sizeof(float));
 
-  cudaMemcpy(jacobian_input_device,jacobian_input_host,jacobian_input_size*sizeof(float),cudaMemcpyHostToDevice);
-  cudaMemcpy(visibility_input_device,visibility_input_host,visibility_input_size*sizeof(char),cudaMemcpyHostToDevice);
-  cudaMemcpy(residual_input_device,residual_input_host,residual_input_size*sizeof(float),cudaMemcpyHostToDevice);
-
-  //  cudaFuncSetCacheConfig (k_jacobianReduceHessianGradient<16,true>, cudaFuncCachePreferShared);
-  //  cudaFuncSetCacheConfig (k_reduceHessianGradient<16,true>, cudaFuncCachePreferShared);
+
+  //  std::cout << "nr_elements = " << nr_elements << std::endl;
+  //  std::cout << "threads = " << threads << std::endl;
+  //  std::cout << "nr_ele_per_thread = " << nr_ele_per_thread << std::endl;
+  //  std::cout << "blocks = " << blocks << std::endl;
+  //  std::cout << "blocks_brent = " << blocks_brent << std::endl;
+  //  std::cout << "total nr threads " << blocks*threads << std::endl;
+  //  std::cout << "brent nr threads " << blocks_brent*threads << std::endl;
+  //  std::cout << "effective elements per thread " << static_cast<double>(nr_elements)/static_cast<double>(blocks*threads) << std::endl;
+  //  std::cout << "brent effective elements per thread " << static_cast<double>(nr_elements)/static_cast<double>(blocks_brent*threads) << std::endl;
 
   // ----------------------- reduce normal
   int output_cache = 1;
+  std::cout << "reducing normal blocks/threads/problem_size: " << blocks << "," <<threads<<","<<nr_elements << std::endl;
   reduceJacobian(nr_elements,threads,blocks,jacobian_input_device,visibility_input_device,residual_input_device,gradient_output,hessian_output);
+
   int problem_size = blocks;
   while(blocks > 1)
   {
     threads = (problem_size < max_threads*2) ? nextPow2((problem_size + 1)/ 2) : max_threads;
     blocks = (problem_size + (threads * 2 - 1)) / (threads * 2);
+    std::cout << "reducing normal blocks/threads/problem_size: " << blocks << "," <<threads<<","<<problem_size << std::endl;
     reduceHessianGradient(problem_size,threads,blocks,gradient_output,hessian_output,gradient_output_2,hessian_output_2);
     output_cache = 2;
     if(blocks > 1)

      break;

    case 64:
      cudaFuncSetCacheConfig (k_reduceHessianGradient<64,false>, cudaFuncCachePreferL1);
      k_reduceHessianGradient<64, false><<< dimGrid, dimBlock >>>(gradient_cache,
                                                                  hessian_cache,

    copyVector<kHessianTriagN>(&s_hessian_data[hessian_index],hessian);
  }

  //__syncthreads();

  if ((_block_size >=   4) && (tid <  2))
  {
