     const size_t last_y = blockIdx.y;
     const size_t last_z = blockIdx.z;
 
-    const size_t n = last_x*gridDim.y*gridDim.z + last_y*gridDim.z + last_z;
-
-    const long int min_out_x = offset_x + last_x*shift_x;
-    const long int min_out_y = offset_y + last_y*shift_y;
-
-    const long int min_x = (min_out_x < 0) ? -min_out_x : 0;
-    const long int min_y = (min_out_y < 0) ? -min_out_y : 0;
-
-    const long int max_x = (length_x + min_out_x > output_size_x) ? output_size_x - min_out_x : length_x;
-    const long int max_y = (length_y + min_out_y > output_size_y) ? output_size_y - min_out_y : length_y;
-
-    float_type inf_cst = expr_inf_cst[n];
-    float_type sup_cst = expr_sup_cst[n];
+    size_t out_z = threadIdx.x;
+
+    __shared__ size_t n;
+
+    __shared__ long int min_out_x;
+    __shared__ long int min_out_y;
+
+    __shared__ long int min_x;
+    __shared__ long int min_y;
+
+    __shared__ long int max_x;
+    __shared__ long int max_y;
+
+    if(out_z == 0)
+    {
+        n = last_x*gridDim.y*gridDim.z + last_y*gridDim.z + last_z;
+
+        min_out_x = offset_x + last_x*shift_x;
+        min_out_y = offset_y + last_y*shift_y;
+
+        min_x = (min_out_x < 0) ? -min_out_x : 0;
+        min_y = (min_out_y < 0) ? -min_out_y : 0;
+
+        max_x = (length_x + min_out_x > output_size_x) ? output_size_x - min_out_x : length_x;
+        max_y = (length_y + min_out_y > output_size_y) ? output_size_y - min_out_y : length_y;
+    }
+
+    __syncthreads();
+
+    float_type inf_cst = 0;
+    float_type sup_cst = 0;
 
     float_type tmp1, tmp2;
     float_type maxRes, maxMul;
 
-    for(long int out_x = min_x; out_x < max_x; out_x++)
-    {
-        for(long int out_y = min_y; out_y < max_y; out_y++)
-        {
-            for(size_t out_z = 0; out_z < output_size_z; out_z++)
+    while(out_z < output_size_z)
+    {
+        for(long int out_x = min_x; out_x < max_x; out_x++)
+        {
+            for(long int out_y = min_y; out_y < max_y; out_y++)
             {
                 size_t mat_out = n*length_x*length_y*output_size_z + out_x*length_y*output_size_z + out_y*output_size_z + out_z;
 
                 const float_type prev_inf_coeff = expr_inf_coeff[mat_out];
                 const float_type prev_sup_coeff = expr_sup_coeff[mat_out];
 
-                if((prev_inf_coeff != 0) || (prev_sup_coeff != 0))
-                {
-                    float_type aux_cst = aux_csts[out_z];
-
-                    elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, prev_inf_coeff, prev_sup_coeff, aux_cst);
-
-                    maxRes = max(abs(inf_cst), abs(sup_cst));
-                    maxMul = max(abs(tmp1), abs(tmp2));
-
-                    inf_cst += tmp1 - (maxRes + maxMul)*ulp - min_denormal;
-                    sup_cst += tmp2 + (maxRes + maxMul)*ulp + min_denormal;
-                }
-            }
-        }
-    }
-
-    res_inf_cst[n] = inf_cst;
-    res_sup_cst[n] = sup_cst;
+                elina_double_interval_mul_cst_coeff_const_expr(&tmp1, &tmp2, prev_inf_coeff, prev_sup_coeff, aux_csts[out_z]);
+
+                maxRes = max(abs(inf_cst), abs(sup_cst));
+                maxMul = max(abs(tmp1), abs(tmp2));
+
+                inf_cst += tmp1 - (maxRes + maxMul)*ulp - min_denormal;
+                sup_cst += tmp2 + (maxRes + maxMul)*ulp + min_denormal;
+            }
+        }
+
+        out_z += blockDim.x;
+    }
+
+    block_reduce_sum(inf_cst, blockDim.x);
+    block_reduce_sum(sup_cst, blockDim.x);
+
+    if(threadIdx.x == 0)
+    {
+        res_inf_cst[n] = inf_cst + expr_inf_cst[n];
+        res_sup_cst[n] = sup_cst + expr_sup_cst[n];
+    }
 }
 
 

const size_t maximum_backstep = 2;

const size_t num_threads = 128;

bool results[90];
bool results_calculated;

        uinf_coeff_tmp = nullptr;
        usup_coeff_tmp = nullptr;

        compute_lb_from_expr<<<num_out_neurons_last_layer, num_threads>>>(lb_array, linf_coeff, lsup_coeff, linf_cst, fp->input_inf + num_in_neurons_0_layer, fp->input_sup + num_in_neurons_0_layer, mu);
        compute_ub_from_expr<<<num_out_neurons_last_layer, num_threads>>>(ub_array, uinf_coeff, usup_coeff, usup_cst, fp->input_inf + num_in_neurons_0_layer, fp->input_sup + num_in_neurons_0_layer, mu);
    }
    else
    {
        compute_lb_from_expr<<<num_out_neurons_last_layer, num_threads>>>(lb_array, linf_coeff, lsup_coeff, linf_cst, fp->input_inf, fp->input_sup, num_in_neurons_0_layer);
        compute_ub_from_expr<<<num_out_neurons_last_layer, num_threads>>>(ub_array, uinf_coeff, usup_coeff, usup_cst, fp->input_inf, fp->input_sup, num_in_neurons_0_layer);
    }

    cudaFree(linf_coeff);
