   {
     char rand_init = (static_cast<float>(rand())/RAND_MAX - visibility_ratio) < 0 ? 1 : 0;
     visibility_input_host[ii] = rand_init;
+    //visibility_input_host[ii] = 1;
     //std::cout << static_cast<int>(visibility_input_host[ii]) << " ";
   }
   //std::cout << std::endl;

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

  // ----------------------- reduce normal
  int output_cache = 1;
  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);
    reduceHessianGradient(problem_size,threads,blocks,gradient_output,hessian_output,gradient_output_2,hessian_output_2);
    output_cache = 2;
    if(blocks > 1)

      break;

    case 64:
      std::cout << "reducing hessian gradient" << std::endl;
      cudaFuncSetCacheConfig (k_reduceHessianGradient<64,false>, cudaFuncCachePreferL1);
      k_reduceHessianGradient<64, false><<< dimGrid, dimBlock >>>(gradient_cache,
                                                                  hessian_cache,
