     {
       hessian_output_array[row*kJacobianSize + col] = hessian_output_array[col*kJacobianSize + row] =  hessian_array_host[index];
     }
-    std::cout << std::endl;
   }
 
   for(unsigned int row = 0; row < kJacobianSize; ++row)

      break;

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

      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)
  {
    cudaMemcpy(hessian_array_host,hessian_output,kHessianTriagN*sizeof(float),cudaMemcpyDeviceToHost);
    cudaMemcpy(gradient_array_host,gradient_output,kJacobianSize*sizeof(float),cudaMemcpyDeviceToHost);
