             float_type tmp1, tmp2;
             elina_double_interval_mul2(&tmp1, &tmp2, old_inf_coeff, old_sup_coeff, 0, ub);
 
-            atomicAdd(&inf_cst[n], tmp1);
-            atomicAdd(&sup_cst[n], tmp1);
-        }
-
-        i += blockDim.x;
+            inf_cst[n] += tmp1;
+            sup_cst[n] += tmp1;
+        }
     }
 }
 

        {
            if(fp->layers[k]->activation==RELU)
            {
                lexpr_replace_relu_bounds<<<num_out_neurons_last_layer, num_threads>>>(linf_coeff, lsup_coeff, linf_cst, lsup_cst, fp->layers[k]->lb_array, fp->layers[k]->ub_array, fp->layers[k]->num_out_neurons, use_area_heuristic);
            }

            std::cout << "DENSE RESIDUAL" << std::endl;

void uexpr_replace_relu_bounds(float_type* __restrict__ inf_coeff, float_type* __restrict__ sup_coeff, float_type* __restrict__ inf_cst, float_type* __restrict__ sup_cst, const float_type* __restrict__ lb_array, const float_type* __restrict__ ub_array, const size_t num_out_neurons_current_layer, const bool use_area_heuristic)
{
    const size_t n = blockIdx.x;

    for(size_t i = 0; i < num_out_neurons_current_layer; i++)
    {
        const size_t a = n*num_out_neurons_current_layer + i;

