   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 max_threads = 128;
 
   int threads = (nr_elements < max_threads*2) ? nextPow2((nr_elements + 1)/ 2) : max_threads;
   int blocks = (nr_elements + (threads * 2 - 1)) / (threads * 2);

  }

  // 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)

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

  // warm up call
  reduceJacobian(nr_elements ,threads ,blocks ,jacobian_input_device ,visibility_input_device ,residual_input_device ,gradient_output ,hessian_output);

  // ----------------------- reduce normal
  int output_cache = 1;
  int problem_size;

  //test static memory
  std::clock_t c_start_gpu_normal = std::clock();
  for(unsigned int tt = 0; tt < nr_iterations;++tt)
  {
    //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);

    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)
      {
        problem_size = blocks;
        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_2,hessian_output_2,gradient_output,hessian_output);
        output_cache = 1;
      }
    }
  }
  std::clock_t c_end_gpu_normal = std::clock();
  double time_gpu_normal = CLOCK_TO_MS(c_start_gpu_normal,c_end_gpu_normal);

  //get result
  float hessian_array_host[kHessianTriagN];
