   //  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;
-  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)
+  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)
     {
-      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::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];

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

  __syncthreads();

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

  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 = 128;

  int threads = (nr_elements < max_threads*2) ? nextPow2((nr_elements + 1)/ 2) : max_threads;
  int blocks = (nr_elements + (threads * 2 - 1)) / (threads * 2);
