         {
             if(fp->layers[k]->activation == RELU)
             {
-                lexpr_replace_relu_bounds<<<1, 1>>>(lexpr_ptr, aux_lb_array, aux_ub_array, 1);
-            }
-
-            expr_from_previous_layer<<<1, 1>>>(lexpr_ptr, aux_expr_array, 1);
+                lexpr_replace_relu_bounds_specialized<<<num_blocks_relu, num_threads>>>(lexpr_array, aux_lb_array, aux_ub_array, 1, num_out_neurons_current_layer);
+            }
+
+            layer_allocate_exprs<<<1, 1>>>(lexpr_array_tmp, 1, num_in_neurons_current_layer);
+
+            expr_from_previous_layer_specialized<<<num_blocks_linear, num_threads>>>(lexpr_array, lexpr_array_tmp, aux_expr_array, 1, num_out_neurons_current_layer, num_in_neurons_current_layer);
+
+            std::swap(lexpr_array, lexpr_array_tmp);
+
+            free_expr_array<<<1, 1>>>(lexpr_array_tmp, num_out_neurons_current_layer);
         }
         else
         {
-            lexpr_replace_maxpool_bounds<<<1, 1>>>(lexpr_ptr, aux_maxpool_lexpr_array, aux_maxpool_uexpr_array, aux_lb_array, aux_ub_array, 1);
+            lexpr_replace_maxpool_bounds<<<1, 1>>>(lexpr_array, aux_maxpool_lexpr_array, aux_maxpool_uexpr_array, aux_lb_array, aux_ub_array, 1);
         }
     }
 
     if((fp->input_lexpr != nullptr) && (fp->input_uexpr != nullptr))
     {
-        replace_input_poly_cons_in_lexpr<<<1, 1>>>(lexpr_ptr, fp->input_inf, fp->input_sup, fp->input_lexpr, fp->input_uexpr, 1);
-    }
-
-    compute_lb_from_expr<<<1, 1>>>(lb_dev, lexpr_ptr, fp->input_inf, fp->input_sup, 1);
-}
-
-
-__global__
-void create_sub_expr(expr_t** sub, const elina_dim_t y, const elina_dim_t x)
-{
-
-    expr_t* expr = (expr_t*) malloc(sizeof(expr_t));
-
-    expr->inf_cst = 0;
-    expr->sup_cst = 0;
-
-    expr->inf_coeff = (double*) malloc(2*sizeof(double));
-    expr->sup_coeff = (double*) malloc(2*sizeof(double));
-    expr->dim = (size_t*) malloc(2*sizeof(size_t));
-    expr->size = 2;
-    expr->type = SPARSE;
-
-    expr->inf_coeff[0] = -1;
-    expr->sup_coeff[0] = 1;
-    expr->dim[0] = y;
-
-    expr->inf_coeff[1] = 1;
-    expr->sup_coeff[1] = -1;
-    expr->dim[1] = x;
-
-    printf("y, %lu\n", expr->dim[0]);
-    printf("x, %lu\n", expr->dim[1]);
-
-    *sub = expr;
+        replace_input_poly_cons_in_lexpr<<<1, 1>>>(lexpr_array, fp->input_inf, fp->input_sup, fp->input_lexpr, fp->input_uexpr, 1);
+    }
+
+    compute_lb_from_expr<<<1, 1>>>(lb_dev, lexpr_array, fp->input_inf, fp->input_sup, 1);
+
+    double lb;
+    cudaMemcpy(&lb, lb_dev, sizeof(double), cudaMemcpyDeviceToHost);
+
+    free_expr_array<<<1, 1>>>(lexpr_array, 1);
+    cudaFree(lexpr_array);
+    cudaFree(lexpr_array_tmp);
+    cudaFree(lb_dev);
+
+    return lb;
 }
 
 

}


void get_lb_using_previous_layers(double* lb_dev, elina_manager_t* man, const fppoly_t* const fp, expr_t** expr)
{
    const size_t numlayers = fp->numlayers;
    expr_t** lexpr_ptr = expr;
    fppoly_internal_t* pr = fppoly_init_from_manager(man, ELINA_FUNID_ASSIGN_LINEXPR_ARRAY);

    for(int k = numlayers - 1; k >= 0; k--)
    {
        expr_t** aux_expr_array = fp->layers[k]->expr_array;
        expr_t** aux_maxpool_lexpr_array = fp->layers[k]->maxpool_lexpr_array;
        expr_t** aux_maxpool_uexpr_array = fp->layers[k]->maxpool_uexpr_array;

    {
        const size_t num_out_neurons_current_layer = fp->layers[k]->num_out_neurons;
        const size_t num_in_neurons_current_layer  = fp->layers[k]->num_in_neurons;
        std::cout << "num_out_neurons_current " << num_out_neurons_current_layer << " num_in_neurons_current " << num_in_neurons_current_layer << std::endl;

        const size_t num_threads = 512;

