diff --git a/Project.toml b/Project.toml
index 659ba2cba..9c0ed4519 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DynamicSparseArrays = "8086fd22-9a0c-46a5-a6c8-6e24676501fe"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
diff --git a/src/Algorithm/Algorithm.jl b/src/Algorithm/Algorithm.jl
index b6a1a7535..d63be923f 100644
--- a/src/Algorithm/Algorithm.jl
+++ b/src/Algorithm/Algorithm.jl
@@ -6,7 +6,7 @@ import TimerOutputs
 
 using ..Coluna, ..ColunaBase, ..MathProg
 
-using Crayons, DynamicSparseArrays, Logging, Parameters, Printf, Random, Statistics
+using Crayons, DynamicSparseArrays, Logging, Parameters, Printf, Random, Statistics, SparseArrays, LinearAlgebra
 
 const TO = TimerOutputs
 const DS = DataStructures
diff --git a/src/Algorithm/colgen.jl b/src/Algorithm/colgen.jl
index ffdc0fe12..18af3f9dc 100644
--- a/src/Algorithm/colgen.jl
+++ b/src/Algorithm/colgen.jl
@@ -146,40 +146,105 @@ function get_units_usage(algo::ColumnGeneration, reform::Reformulation)
     return units_usage
 end
 
+############################################################################################
+# Information extracted to speed-up some computations.
+############################################################################################
+function _submatrix(
+    form::Formulation, 
+    keep_constr::Function, 
+    keep_var::Function,
+)
+    matrix = getcoefmatrix(form)
+    constr_ids = ConstrId[]
+    var_ids = VarId[]
+    nz = Float64[]
+    for constr_id in Iterators.filter(keep_constr, Iterators.keys(getconstrs(form)))
+        for (var_id, coeff) in @view matrix[constr_id, :]
+            if keep_var(var_id)
+                push!(constr_ids, constr_id)
+                push!(var_ids, var_id)
+                push!(nz, coeff)
+            end
+        end
+    end
+    return dynamicsparse(
+        constr_ids, var_ids, nz, ConstrId(Coluna.MAX_NB_ELEMS), VarId(Coluna.MAX_NB_ELEMS)
+    )
+end
+
+"""
+Extracted information to speed-up calculation of reduced costs of master variables.
+We extract from the master only the information we need to compute the reduced cost of DW 
+subproblem variables:
+- `c` contains the perenial cost of DW subproblem representative variables
+- `A` is a submatrix of the master coefficient matrix that involves only DW subproblem
+  representative variables.
+
+Calculation is `c - transpose(A) * master_lp_dual_solution`.
+"""
 struct ReducedCostsCalculationHelper
-    length::Int
-    dwspvarids::Vector{VarId}
-    perencosts::Vector{Float64}
-    dwsprep_coefmatrix::DynamicSparseArrays.Transposed{DynamicSparseArrays.DynamicSparseMatrix{Coluna.MathProg.ConstrId, Coluna.MathProg.VarId, Float64}}
+    c::SparseVector{Float64,VarId}
+    A::DynamicSparseMatrix{ConstrId,VarId,Float64}
 end
 
-# Precompute information to speed-up calculation of reduced costs of original variables
-# in the master of a given reformulation.
-function ReducedCostsCalculationHelper(reform::Reformulation)
-    dwspvarids = VarId[]
-    perencosts = Float64[]
+function ReducedCostsCalculationHelper(master)
+    dwspvar_ids = VarId[]
+    peren_costs = Float64[]
 
-    master = getmaster(reform)
-    for (varid, _) in getvars(master)
-        if iscuractive(master, varid) && getduty(varid) <= AbstractMasterRepDwSpVar
-            push!(dwspvarids, varid)
-            push!(perencosts, getcurcost(master, varid))
+    for var_id in Iterators.keys(getvars(master))
+        if iscuractive(master, var_id) && getduty(var_id) <= AbstractMasterRepDwSpVar
+            push!(dwspvar_ids, var_id)
+            push!(peren_costs, getcurcost(master, var_id))
         end
     end
 
-    master_coefmatrix = getcoefmatrix(master)
-    dwsprep_coefmatrix = dynamicsparse(ConstrId, VarId, Float64)
-    for (constrid, _) in getconstrs(master)
-        for (varid, coeff) in @view master_coefmatrix[constrid, :]
-            if getduty(varid) <= AbstractMasterRepDwSpVar
-                dwsprep_coefmatrix[constrid, varid] = coeff
-            end
-        end
-    end
-    closefillmode!(dwsprep_coefmatrix)
-    return ReducedCostsCalculationHelper(
-        length(dwspvarids), dwspvarids, perencosts, transpose(dwsprep_coefmatrix)
+    dwsprep_costs = sparsevec(dwspvar_ids, peren_costs, Coluna.MAX_NB_ELEMS)
+    dwsprep_coefmatrix = _submatrix(
+        master, 
+        constr_id -> true,
+        var_id -> getduty(var_id) <= AbstractMasterRepDwSpVar
     )
+    return ReducedCostsCalculationHelper(dwsprep_costs, dwsprep_coefmatrix)
+end
+
+"""
+Precompute information to speed-up calculation of subgradient of master variables.
+We extract from the master follwowing information:
+- `a` contains the perenial rhs of all master constraints except convexity constraints;
+- `A` is a submatrix of the master coefficient matrix that involves only representative of
+  original variables (pure master vars + DW subproblem represtative vars) 
+
+Calculation is `a - A * (m .* z)`
+where :
+ - `m` contains a multiplicity factor for each variable involved in the calculation
+       (lower or upper sp multiplicity depending on variable reduced cost);
+ - `z` is the concatenation of the solution to the master (for pure master vars) and pricing
+       subproblems (for DW subproblem represtative vars).
+
+Operation `m .* z` "mimics" a solution in the original space.
+"""
+struct SubgradientCalculationHelper
+    a::SparseVector{Float64,ConstrId}
+    A::DynamicSparseMatrix{ConstrId,VarId,Float64}
+end
+
+function SubgradientCalculationHelper(master)
+    constr_ids = ConstrId[]
+    constr_rhs = Float64[]
+    for (constr_id, constr) in getconstrs(master)
+        if !(getduty(constr_id) <= MasterConvexityConstr) && 
+           iscuractive(master, constr) && isexplicit(master, constr)
+            push!(constr_ids, constr_id)
+            push!(constr_rhs, getcurrhs(master, constr_id))
+        end 
+    end
+    a = sparsevec(constr_ids, constr_rhs, Coluna.MAX_NB_ELEMS)
+    A = _submatrix(
+        master, 
+        constr_id -> !(getduty(constr_id) <= MasterConvexityConstr),
+        var_id -> getduty(var_id) <= MasterPureVar || getduty(var_id) <= MasterRepPricingVar
+    )
+    return SubgradientCalculationHelper(a, A)
 end
 
 function run!(algo::ColumnGeneration, env::Env, reform::Reformulation, input::OptimizationState)
@@ -288,6 +353,31 @@ end
 
 set_bestcol_id!(spinfo::SubprobInfo, varid::VarId) = spinfo.bestcol_id = varid
 
+function compute_subgradient_contribution(
+    algo::ColumnGeneration, stabunit::ColGenStabilizationUnit, master::Formulation,
+    puremastervars::Vector{Pair{VarId,Float64}}, spinfos::Dict{FormId,SubprobInfo}
+)
+    sense = getobjsense(master)
+    var_ids = VarId[]
+    var_vals = Float64[]
+
+    for (var_id, mult) in puremastervars
+        push!(var_ids, var_id)
+        push!(var_vals, mult) 
+    end
+
+    for (_, spinfo) in spinfos
+        iszero(spinfo.ub) && continue
+        mult = improving_red_cost(getbound(spinfo.bestsol), algo, sense) ? spinfo.ub : spinfo.lb
+        for (sp_var_id, sp_var_val) in spinfo.bestsol
+            push!(var_ids, sp_var_id)
+            push!(var_vals, sp_var_val * mult)
+        end
+    end
+    return sparsevec(var_ids, var_vals)
+end
+
+
 function compute_db_contributions!(
     spinfo::SubprobInfo, dualbound::DualBound{MaxSense}, primalbound::PrimalBound{MaxSense}
 )
@@ -425,22 +515,18 @@ end
 
 # this method must be redefined if subproblem is a custom model
 function updatemodel!(
-    form::Formulation, repr_vars_red_costs::Dict{VarId, Float64}, ::DualSolution
+    form::Formulation, redcosts, ::DualSolution
 )
-    for (varid, _) in getvars(form)
-        setcurcost!(form, varid, get(repr_vars_red_costs, varid, 0.0))
+    for (var_id, _) in getvars(form)
+        setcurcost!(form, var_id, redcosts[var_id])
     end
     return
 end
 
 function updatereducedcosts!(
-    reform::Reformulation, redcostshelper::ReducedCostsCalculationHelper, masterdualsol::DualSolution
+    reform::Reformulation, helper::ReducedCostsCalculationHelper, masterdualsol::DualSolution
 )
-    redcosts = Dict{VarId,Float64}()
-    result = redcostshelper.dwsprep_coefmatrix * masterdualsol.solution.sol
-    for (i, varid) in enumerate(redcostshelper.dwspvarids)
-        redcosts[varid] = redcostshelper.perencosts[i] - get(result, varid, 0.0)
-    end
+    redcosts = helper.c - transpose(helper.A) * masterdualsol
     for (_, spform) in get_dw_pricing_sps(reform)
         updatemodel!(spform, redcosts, masterdualsol)
     end
@@ -609,43 +695,6 @@ function update_lagrangian_dual_bound!(
     return
 end
 
-function compute_subgradient_contribution(
-    algo::ColumnGeneration, stabunit::ColGenStabilizationUnit, master::Formulation,
-    puremastervars::Vector{Pair{VarId,Float64}}, spinfos::Dict{FormId,SubprobInfo}
-)
-    sense = getobjsense(master)
-    constrids = ConstrId[]
-    constrvals = Float64[]
-
-    if subgradient_is_needed(stabunit, algo.smoothing_stabilization)
-        master_coef_matrix = getcoefmatrix(master)
-
-        for (varid, mult) in puremastervars
-            for (constrid, var_coeff) in @view master_coef_matrix[:,varid]
-                push!(constrids, constrid)
-                push!(constrvals, var_coeff * mult)
-            end
-        end
-
-        for (_, spinfo) in spinfos
-            iszero(spinfo.ub) && continue
-            mult = improving_red_cost(getbound(spinfo.bestsol), algo, sense) ? spinfo.ub : spinfo.lb
-            for (sp_var_id, sp_var_val) in spinfo.bestsol
-                for (master_constrid, sp_var_coef) in @view master_coef_matrix[:,sp_var_id]
-                    if !(getduty(master_constrid) <= MasterConvexityConstr)
-                        push!(constrids, master_constrid)
-                        push!(constrvals, sp_var_coef * sp_var_val * mult)
-                    end
-                end
-            end
-        end
-    end
-
-    return DualSolution(
-        master, constrids, constrvals, VarId[], Float64[], ActiveBound[], 0.0, 
-        UNKNOWN_SOLUTION_STATUS
-    )
-end
 
 function move_convexity_constrs_dual_values!(
     spinfos::Dict{FormId,SubprobInfo}, dualsol::DualSolution
@@ -710,7 +759,8 @@ function cg_main_loop!(
         spinfos[spid] = SubprobInfo(reform, spid)
     end
 
-    redcostshelper = ReducedCostsCalculationHelper(reform)
+    redcostshelper = ReducedCostsCalculationHelper(getmaster(reform))
+    sg_helper = SubgradientCalculationHelper(getmaster(reform))
     iteration = 0
     essential_cuts_separated = false
 
@@ -793,7 +843,8 @@ function cg_main_loop!(
                         if phase == 2 # because the new cuts may make the master infeasible
                             return false, true
                         end
-                        redcostshelper = ReducedCostsCalculationHelper(reform)
+                        redcostshelper = ReducedCostsCalculationHelper(getmaster(reform))
+                        sg_helper = SubgradientCalculationHelper(getmaster(reform))
                     end
                 end
             end
@@ -843,6 +894,7 @@ function cg_main_loop!(
                 TO.@timeit Coluna._to "Smoothing update" begin
                     smooth_dual_sol = update_stab_after_gencols!(
                         stabunit, algo.smoothing_stabilization, nb_new_col, lp_dual_sol, smooth_dual_sol,
+                        sg_helper,
                         compute_subgradient_contribution(algo, stabunit, masterform, pure_master_vars, spinfos)
                     )
                 end
diff --git a/src/Algorithm/colgenstabilization.jl b/src/Algorithm/colgenstabilization.jl
index 14b9a917e..7500cb51a 100644
--- a/src/Algorithm/colgenstabilization.jl
+++ b/src/Algorithm/colgenstabilization.jl
@@ -70,66 +70,6 @@ function init_stab_before_colgen_loop!(unit::ColGenStabilizationUnit)
     return
 end
 
-function componentwisefunction(in_dual_sol::DualSolution, out_dual_sol::DualSolution, f::Function)
-    form = getmodel(out_dual_sol)
-    constrids = Vector{ConstrId}()
-    constrvals = Vector{Float64}()
-    value::Float64 = 0.0
-    out_next = iterate(out_dual_sol)
-    in_next = iterate(in_dual_sol)
-    while out_next !== nothing || in_next !== nothing
-        if out_next !== nothing 
-            ((out_constrid, out_val), out_state) = out_next
-            if in_next !== nothing
-                ((in_constrid, in_val), in_state) = in_next
-                if in_constrid < out_constrid
-                    value = f(in_val, 0.0)
-                    push!(constrids, in_constrid)
-                    push!(constrvals, value)
-                    in_next = iterate(in_dual_sol, in_state)
-                elseif out_constrid < in_constrid    
-                    value = f(0.0, out_val)
-                    push!(constrids, out_constrid)
-                    push!(constrvals, value)
-                    out_next = iterate(out_dual_sol, out_state)    
-                else
-                    push!(constrids, out_constrid)
-                    value = f(in_val, out_val)
-                    push!(constrvals, value)
-                    in_next = iterate(in_dual_sol, in_state)
-                    out_next = iterate(out_dual_sol, out_state)    
-                end
-            else
-                value = f(0.0, out_val)
-                push!(constrids, out_constrid)
-                push!(constrvals, value)
-                out_next = iterate(out_dual_sol, out_state)    
-            end
-        else    
-            ((in_constrid, in_val), in_state) = in_next
-            value = f(in_val, 0.0)
-            push!(constrids, in_constrid)
-            push!(constrvals, value)
-            in_next = iterate(in_dual_sol, in_state)
-        end
-    end
-    return (constrids, constrvals)
-end
-
-function linear_combination(in_dual_sol::DualSolution, out_dual_sol::DualSolution, coeff::Float64)
-    return coeff * in_dual_sol + (1.0 - coeff) * out_dual_sol
-
-    # form = getmodel(in_dual_sol)
-    # bound = 0.0
-    # for (i, constrid) in enumerate(constrids)
-    #     bound += constrvals[i] * getcurrhs(form, constrid) 
-    # end
-    # return DualSolution(
-    #     form, constrids, constrvals, VarId[], Float64[], ActiveBound[], bound, 
-    #     UNKNOWN_FEASIBILITY
-    # )
-end
-
 function update_stab_after_rm_solve!(
     unit::ColGenStabilizationUnit, smoothparam::Float64, lp_dual_sol::DualSolution
 )
@@ -150,63 +90,31 @@ function update_stab_after_rm_solve!(
         # in this case Lagrangian bound calculation is simplified in col.gen.
         # (we use the fact that the contribution of pure master variables 
         #  is included in the value of the LP dual solution)
-        # thus, LP dual solution should be retured, as linear_combination() 
+        # thus, LP dual solution should be retured, as linear combination
         # does not include pure master variables contribution to the bound
         return lp_dual_sol
     end
-    return linear_combination(unit.stabcenter, lp_dual_sol, unit.curalpha)
-end
-
-function norm(dualsol::DualSolution)
-    product_sum = 0.0
-    for (constrid, val) in dualsol
-        product_sum += val * val
-    end
-    return sqrt(product_sum)
+    return unit.curalpha * unit.stabcenter + (1.0 - unit.curalpha) * lp_dual_sol
 end
 
 function update_alpha_automatically!(
     unit::ColGenStabilizationUnit, nb_new_col::Int64, lp_dual_sol::DualSolution{M},  
-    smooth_dual_sol::DualSolution{M}, subgradient_contribution::DualSolution{M}
+    smooth_dual_sol::DualSolution{M}, h, subgradient_contribution
 ) where {M}    
 
     master = getmodel(lp_dual_sol)
 
     # first we calculate the in-sep direction
-    constrids, constrvals = componentwisefunction(smooth_dual_sol, unit.stabcenter, -)
-    in_sep_direction = DualSolution(
-        master, constrids, constrvals, VarId[], Float64[], ActiveBound[], 0.0, 
-        UNKNOWN_FEASIBILITY
-    )
+    in_sep_direction = smooth_dual_sol - unit.stabcenter
     in_sep_dir_norm = norm(in_sep_direction)
 
-    # we initialize the subgradient with the right-hand-side of all master constraints
-    # except the convexity constraints
-    constrids = Vector{ConstrId}()
-    constrrhs = Vector{Float64}() 
-    for (constrid, constr) in getconstrs(master)
-        if !(getduty(constrid) <= MasterConvexityConstr) && 
-           iscuractive(master, constr) && isexplicit(master, constr)
-            push!(constrids, constrid)
-            push!(constrrhs, getcurrhs(master, constrid))
-        end 
-    end
-    subgradient = DualSolution(
-        master, constrids, constrrhs, VarId[], Float64[], ActiveBound[], 0.0, 
-        UNKNOWN_FEASIBILITY
-    )
-    
-    # we calculate the subgradient at the sep point
-    for (constrid, value) in subgradient_contribution
-        subgradient[constrid] = subgradient[constrid] - value
-    end
+    subgradient = h.a - h.A * subgradient_contribution
     subgradient_norm = norm(subgradient)
 
     # we now calculate the angle between the in-sep direction and the subgradient 
-    constrids, constrvals = componentwisefunction(in_sep_direction, subgradient, *)
-    angle = sum(constrvals) / (in_sep_dir_norm * subgradient_norm)
+    angle = (transpose(in_sep_direction) * subgradient) / (in_sep_dir_norm * subgradient_norm)
     if getobjsense(master) == MaxSense
-        angle = -angle
+        angle *= -1
     end
 
     # we modify the alpha parameter based on the calculated angle
@@ -220,15 +128,15 @@ end
 
 function update_stab_after_gencols!(
     unit::ColGenStabilizationUnit, smoothparam::Float64, nb_new_col::Int64, 
-    lp_dual_sol::DualSolution{M}, smooth_dual_sol::DualSolution{M}, 
-    subgradient_contribution::DualSolution{M}
+    lp_dual_sol::DualSolution{M}, smooth_dual_sol::DualSolution{M}, h,
+    subgradient_contribution
 ) where {M}
 
     iszero(smoothparam) && return nothing
 
     if smoothparam == 1.0 && unit.nb_misprices == 0
         update_alpha_automatically!(
-            unit, nb_new_col, lp_dual_sol, smooth_dual_sol, subgradient_contribution
+            unit, nb_new_col, lp_dual_sol, smooth_dual_sol, h, subgradient_contribution
         )
     end
 
@@ -249,8 +157,7 @@ function update_stab_after_gencols!(
         # stabilization is deactivated, thus we need to return the original LP dual solution
         return lp_dual_sol
     end
-
-    return linear_combination(unit.stabcenter, lp_dual_sol, unit.curalpha)
+    return unit.curalpha * unit.stabcenter + (1.0 - unit.curalpha) * lp_dual_sol
 end
 
 function update_stability_center!(
diff --git a/src/Coluna.jl b/src/Coluna.jl
index aba37cbfd..824282257 100644
--- a/src/Coluna.jl
+++ b/src/Coluna.jl
@@ -19,6 +19,8 @@ const DEF_OPTIMALITY_RTOL = 1e-9
 
 const TOL = 1e-8 # if - ϵ_tol < val < ϵ_tol, we consider val = 0
 const TOL_DIGITS = 8 # because round(val, digits = n) where n is from 1e-n
+
+const MAX_NB_ELEMS = typemax(Int32) # max number of variables or constraints.
 ###
 
 # submodules
diff --git a/src/ColunaBase/solsandbounds.jl b/src/ColunaBase/solsandbounds.jl
index 551c64db2..75b735c55 100644
--- a/src/ColunaBase/solsandbounds.jl
+++ b/src/ColunaBase/solsandbounds.jl
@@ -285,7 +285,7 @@ Create a solution to the `model`. Other arguments are:
 function Solution{Mo,De,Va}(
     model::Mo, decisions::Vector{De}, values::Vector{Va}, solution_value::Float64, status::SolutionStatus
 ) where {Mo<:AbstractModel,De,Va}
-    sol = sparsevec(decisions, values, typemax(De))
+    sol = sparsevec(decisions, values, Coluna.MAX_NB_ELEMS)
     return Solution(model, solution_value, status, sol)
 end
 
diff --git a/src/MathProg/MathProg.jl b/src/MathProg/MathProg.jl
index 921e4f3bd..b6434d3be 100644
--- a/src/MathProg/MathProg.jl
+++ b/src/MathProg/MathProg.jl
@@ -9,7 +9,7 @@ using ..ColunaBase
 
 import Base: haskey, length, iterate, diff, delete!, contains, setindex!, getindex, view
 
-using DynamicSparseArrays, SparseArrays, Logging, Printf
+using DynamicSparseArrays, SparseArrays, Logging, Printf, LinearAlgebra
 
 const BD = BlockDecomposition
 const ClB = ColunaBase
diff --git a/src/MathProg/solutions.jl b/src/MathProg/solutions.jl
index 542289166..5eec94ae2 100644
--- a/src/MathProg/solutions.jl
+++ b/src/MathProg/solutions.jl
@@ -1,9 +1,54 @@
+############################################################################################
 # MathProg > Solutions
 # Representations of the primal & dual solutions to a MILP formulation
+############################################################################################
 
+"Supertype for solutions operated by Coluna."
 abstract type AbstractSolution end
 
-# Primal Solution
+# The API for `AbstractSolution` is not very clear yet.
+
+# Redefine methods from ColunaBase to access the formulation, the value, the
+# status of a Solution, and other specific information
+ColunaBase.getmodel(s::AbstractSolution) = getmodel(s.solution)
+ColunaBase.getvalue(s::AbstractSolution) = getvalue(s.solution)
+ColunaBase.getbound(s::AbstractSolution) = getbound(s.solution)
+ColunaBase.getstatus(s::AbstractSolution) = getstatus(s.solution)
+
+Base.length(s::AbstractSolution) = length(s.solution)
+Base.get(s::AbstractSolution, id, default) = get(s.solution, id, default)
+Base.getindex(s::AbstractSolution, id) = getindex(s.solution, id)
+Base.setindex!(s::AbstractSolution, val, id) = setindex!(s.solution, val, id)
+
+# Iterating over a PrimalSolution or a DualSolution is similar to iterating over
+# ColunaBase.Solution
+Base.iterate(s::AbstractSolution) = iterate(s.solution)
+Base.iterate(s::AbstractSolution, state) = iterate(s.solution, state)
+
+function contains(sol::AbstractSolution, f::Function)
+    for (elemid, _) in sol
+        f(elemid) && return true
+    end
+    return false
+end
+
+function _sols_from_same_model(sols::NTuple{N, S}) where {N,S<:AbstractSolution}
+    for i in 2:length(sols)
+        getmodel(sols[i-1]) != getmodel(sols[i]) && return false
+    end
+    return true
+end
+
+# To check if a solution is part of solutions from the pool.
+Base.:(==)(v1::DynamicMatrixColView, v2::AbstractSolution) = v1 == v2.solution
+
+# To allocate an array with size equals to the number of non-zero elements when using
+# "generation" syntax.
+Base.length(gen::Base.Generator{<:AbstractSolution}) = nnz(gen.iter.solution)
+
+############################################################################################
+# Primal Solution                                                                          
+############################################################################################
 struct PrimalSolution{M} <: AbstractSolution
     solution::Solution{M,VarId,Float64}
     custom_data::Union{Nothing, BlockDecomposition.AbstractCustomData}
@@ -40,7 +85,42 @@ end
 
 Base.copy(s::P) where {P<:PrimalSolution}= P(copy(s.solution), copy(s.custom_data))
 
+function Base.isinteger(sol::PrimalSolution)
+    for (vc_id, val) in sol
+        if getperenkind(getmodel(sol), vc_id) !== Continuous && abs(round(val) - val) > 1e-5
+            return false
+        end
+    end
+    return true
+end
+
+function Base.isless(s1::PrimalSolution, s2::PrimalSolution)
+    getobjsense(getmodel(s1)) == MinSense && return s1.solution.bound > s2.solution.bound
+    return s1.solution.bound < s2.solution.bound
+end
+
+# Method `cat` is not implemented for a set of DualSolutions because @guimarqu don't know 
+# how to concatenate var red cost of a variable if both bounds are active in different 
+# solutions and because we don't need it for now.
+function Base.cat(sols::PrimalSolution...)
+    if !_sols_from_same_model(sols)
+        error("Cannot concatenate solutions not attached to the same model.")
+    end
+
+    ids = VarId[]
+    vals = Float64[]
+    for sol in sols, (id, value) in sol
+        push!(ids, id)
+        push!(vals, value)
+    end
+    return PrimalSolution(
+        getmodel(sols[1]), ids, vals, sum(getvalue.(sols)), getstatus(sols[1])
+    )
+end
+
+############################################################################################
 # Dual Solution
+############################################################################################
 
 # Indicate whether the active bound of a variable is the lower or the upper one.
 @enum ActiveBound LOWER UPPER
@@ -96,83 +176,11 @@ Base.copy(s::D) where {D<:DualSolution} = D(copy(s.solution), copy(s.var_redcost
 
 get_var_redcosts(s::DualSolution) = s.var_redcosts
 
-# Redefine methods from ColunaBase to access the formulation, the value, the
-# status of a Solution, and other specific information
-ColunaBase.getmodel(s::AbstractSolution) = getmodel(s.solution)
-ColunaBase.getvalue(s::AbstractSolution) = getvalue(s.solution)
-ColunaBase.getbound(s::AbstractSolution) = getbound(s.solution)
-ColunaBase.getstatus(s::AbstractSolution) = getstatus(s.solution)
-
-Base.length(s::AbstractSolution) = length(s.solution)
-Base.get(s::AbstractSolution, id, default) = get(s.solution, id, default)
-Base.getindex(s::AbstractSolution, id) = getindex(s.solution, id)
-Base.setindex!(s::AbstractSolution, val, id) = setindex!(s.solution, val, id)
-
-# Iterating over a PrimalSolution or a DualSolution is similar to iterating over
-# ColunaBase.Solution
-Base.iterate(s::AbstractSolution) = iterate(s.solution)
-Base.iterate(s::AbstractSolution, state) = iterate(s.solution, state)
-
-function Base.isinteger(sol::PrimalSolution)
-    for (vc_id, val) in sol
-        if getperenkind(getmodel(sol), vc_id) !== Continuous && abs(round(val) - val) > 1e-5
-            return false
-        end
-    end
-    return true
-end
-
-function Base.isless(s1::PrimalSolution, s2::PrimalSolution)
-    getobjsense(getmodel(s1)) == MinSense && return s1.solution.bound > s2.solution.bound
-    return s1.solution.bound < s2.solution.bound
-end
-
 function Base.isless(s1::DualSolution, s2::DualSolution)
     getobjsense(getmodel(s1)) == MinSense && return s1.solution.bound < s2.solution.bound
     return s1.solution.bound > s2.solution.bound
 end
 
-function contains(sol::AbstractSolution, f::Function)
-    for (elemid, _) in sol
-        f(elemid) && return true
-    end
-    return false
-end
-
-function _sols_from_same_model(sols::NTuple{N, S}) where {N,S<:AbstractSolution}
-    for i in 2:length(sols)
-        getmodel(sols[i-1]) != getmodel(sols[i]) && return false
-    end
-    return true
-end
-
-# Method `cat` is not implemented for a set of DualSolutions because @guimarqu don't know 
-# how to concatenate var red cost of a variable if both bounds are active in different 
-# solutions and because we don't need it for now.
-function Base.cat(sols::PrimalSolution...)
-    if !_sols_from_same_model(sols)
-        error("Cannot concatenate solutions not attached to the same model.")
-    end
-
-    ids = VarId[]
-    vals = Float64[]
-    for sol in sols, (id, value) in sol
-        push!(ids, id)
-        push!(vals, value)
-    end
-    return PrimalSolution(
-        getmodel(sols[1]), ids, vals, sum(getvalue.(sols)), getstatus(sols[1])
-    )
-end
-
-function Base.print(io::IO, form::AbstractFormulation, sol::Solution)
-    println(io, "Solution")
-    for (id, val) in sol
-        println(io, getname(form, id), " = ", val)
-    end
-    return
-end
-
 function Base.show(io::IO, solution::DualSolution{M}) where {M}
     println(io, "Dual solution")
     for (constrid, value) in solution
@@ -192,28 +200,23 @@ function Base.show(io::IO, solution::PrimalSolution{M}) where {M}
     Printf.@printf(io, "└ value = %.2f \n", getvalue(solution))
 end
 
+############################################################################################
+# Linear Algebra
+############################################################################################
 
-# To check if a solution is part of solutions from the pool.
-Base.:(==)(v1::DynamicMatrixColView, v2::AbstractSolution) = v1 == v2.solution
+# op(::S, ::S) has return type `S` for op ∈ (:+, :-) and S <: AbstractSolution 
 
-# To allocate an array with size equals to the number of non-zero elements when using
-# "generation" syntax.
-Base.length(gen::Base.Generator{<:AbstractSolution}) = nnz(gen.iter.solution)
-
-### Math operation
-
-## op(::S, ::S) has return type `S` for op ∈ (:+, :-) and S <: AbstractSolution 
-_math_op_constructor(::Type{S}, form, varids, varvals, cost) where {S<:PrimalSolution} =
+_math_op_constructor(::Type{S}, form::F, varids, varvals, cost) where {S<:PrimalSolution,F} =
     PrimalSolution(form, varids, varvals, cost, ClB.UNKNOWN_SOLUTION_STATUS)
 
-_math_op_constructor(::Type{<:S}, form, constrids, constrvals, cost) where {S<:DualSolution} = 
+_math_op_constructor(::Type{<:S}, form::F, constrids, constrvals, cost) where {S<:DualSolution,F} = 
     DualSolution(form, constrids, constrvals, [], [], [], cost, ClB.UNKNOWN_SOLUTION_STATUS)
 
 _math_op_cost(::Type{<:S}, form, varids, varvals) where {S<:PrimalSolution} = 
-    mapreduce(((id,val),) -> getcurcost(form, id) * val, +, Iterators.zip(varids, varvals))
+    mapreduce(((id,val),) -> getcurcost(form, id) * val, +, Iterators.zip(varids, varvals); init = 0.0)
 
 _math_op_cost(::Type{<:S}, form, constrids, constrvals) where {S<:DualSolution} =
-    mapreduce(((id, val),) -> getcurrhs(form, id) * val, +, Iterators.zip(constrids, constrvals))
+    mapreduce(((id, val),) -> getcurrhs(form, id) * val, +, Iterators.zip(constrids, constrvals); init = 0.0)
 
 function Base.:(*)(a::Real, s::S) where {S<:AbstractSolution}
     ids, vals = findnz(a * s.solution.sol)
@@ -232,14 +235,35 @@ for op in (:+, :-)
     end
 end
 
-## *(::M, ::S) has return type `SparseVector` for:
-##  - M <: DynamicSparseMatrix
-##  - S <: AbstractSolution
+# transpose
+struct Transposed{S<:AbstractSolution}
+    sol::S
+end
 
-## We don't support operation with classic sparse matrix because row and col ids
-## must be of the same type. 
-## In Coluna, we use VarId to index the cols and 
-## ConstrId to index the rows.
+Base.transpose(s::AbstractSolution) = Transposed(s)
+
+Base.:(*)(s1::Transposed{S}, s2::S) where {S<:AbstractSolution} =
+    transpose(s1.sol.solution.sol) * s2.solution.sol
+
+function Base.:(*)(s::Transposed{<:AbstractSolution}, vec::SparseVector)
+    # We multiply two sparse vectors that may have different sizes.
+    sol_vec = s.sol.solution.sol
+    len = Coluna.MAX_NB_ELEMS
+    vec1 = sparsevec(findnz(sol_vec)..., len)
+    vec2 = sparsevec(findnz(vec)..., len)
+    return transpose(vec1) * vec2
+end
+
+# *(::M, ::S) has return type `SparseVector` for:
+#  - M <: DynamicSparseMatrix
+#  - S <: AbstractSolution
+
+# We don't support operation with classic sparse matrix because row and col ids
+# must be of the same type. 
+# In Coluna, we use VarId to index the cols and 
+# ConstrId to index the rows.
 
 Base.:(*)(m::DynamicSparseMatrix, s::AbstractSolution) = m * s.solution.sol
-Base.:(*)(m::DynamicSparseArrays.Transposed{<:DynamicSparseMatrix}, s::AbstractSolution) = m * s.solution.sol
\ No newline at end of file
+Base.:(*)(m::DynamicSparseArrays.Transposed{<:DynamicSparseMatrix}, s::AbstractSolution) = m * s.solution.sol
+
+LinearAlgebra.norm(s::AbstractSolution) = norm(s.solution.sol)
\ No newline at end of file
diff --git a/src/MathProg/vcids.jl b/src/MathProg/vcids.jl
index 07c08447b..ac65fa5b0 100644
--- a/src/MathProg/vcids.jl
+++ b/src/MathProg/vcids.jl
@@ -58,7 +58,7 @@ Base.hash(a::Id, h::UInt) = hash(a.uid, h)
 Base.zero(I::Type{Id{VC}}) where {VC} = I(0)
 Base.zero(::Id{VC}) where {VC} = Id{VC}(0)
 Base.one(I::Type{Id{VC}}) where {VC} = I(1)
-Base.typemax(I::Type{Id{VC}}) where {VC} = I(typemax(Int))
+Base.typemax(I::Type{Id{VC}}) where {VC} = I(Coluna.MAX_NB_ELEMS)
 Base.isequal(a::Id{VC}, b::Id{VC}) where {VC} = isequal(a.uid, b.uid)
 
 Base.promote_rule(::Type{T}, ::Type{<:Id}) where {T<:Integer} = T
diff --git a/test/runtests.jl b/test/runtests.jl
index b54c0c361..44f31988b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,6 @@ include("ColunaTests.jl")
 retest(Coluna, ColunaTests)
 
 # Run a specific test:
-# retest(ColunaTests, "MathProg - variable and constraint ids")
+#retest(ColunaTests, "strong branching")
 
 
diff --git a/test/unit/ColunaBase/solsandbounds.jl b/test/unit/ColunaBase/solsandbounds.jl
index 43e78526f..4557f7717 100644
--- a/test/unit/ColunaBase/solsandbounds.jl
+++ b/test/unit/ColunaBase/solsandbounds.jl
@@ -307,7 +307,7 @@ end
         dict_sol = Dict(1 => 2.0, 2 => 3.0, 3 => 4.0)
         primal_sol = Solution(model, collect(keys(dict_sol)), collect(values(dict_sol)), 0.0, ClB.FEASIBLE_SOL)
         
-        @test length(primal_sol) == typemax(Int)
+        @test length(primal_sol) == typemax(Coluna.MAX_NB_ELEMS)
         @test nnz(primal_sol) == 3
         @test primal_sol[1] == 2.0
         primal_sol[1] = 5.0 # change the value
diff --git a/test/unit/MathProg/solutions.jl b/test/unit/MathProg/solutions.jl
index 847913f2f..082bbb5c0 100644
--- a/test/unit/MathProg/solutions.jl
+++ b/test/unit/MathProg/solutions.jl
@@ -100,7 +100,30 @@
         @test findnz(c.solution.sol) == ([ClMP.getid(var3)], [-4.0])
     end
 
-    @testset "lin alg 2 - spMv" begin
+    @testset "lin alg 2 - transpose" begin
+        form = ClMP.create_formulation!(
+            Env{ClMP.VarId}(Coluna.Params()), ClMP.Original(), obj_sense = Coluna.MathProg.MaxSense
+        )
+        var1 = ClMP.setvar!(form, "var1", ClMP.OriginalVar; cost = 1.0)
+        var2 = ClMP.setvar!(form, "var2", ClMP.OriginalVar; cost = 0.5)
+        var3 = ClMP.setvar!(form, "var3", ClMP.OriginalVar; cost = -0.25)
+
+        nzinds1 = [ClMP.getid(var1), ClMP.getid(var2)]
+        nzvals1 = [12.0, 4.0]
+        vec1 = sparsevec(ClMP.getuid.(nzinds1), nzvals1, 4)
+        primalsol1 = ClMP.PrimalSolution(form, nzinds1, nzvals1, 14.0, ClB.FEASIBLE_SOL)
+
+        nzinds2 = [ClMP.getid(var1), ClMP.getid(var2), ClMP.getid(var3)]
+        nzvals2 = [1.0, 2.0, 4.0]
+        vec2 = sparsevec(ClMP.getuid.(nzinds2), nzvals2, 4)
+        primalsol2 = ClMP.PrimalSolution(form, nzinds2, nzvals2, 1.0, ClB.FEASIBLE_SOL)
+
+        a = transpose(vec1) * vec2
+        b = transpose(primalsol1) * primalsol2
+        @test a == b
+    end
+
+    @testset "lin alg 3 - spMv" begin
         form = ClMP.create_formulation!(
             Env{ClMP.VarId}(Coluna.Params()), ClMP.Original(), obj_sense = Coluna.MathProg.MaxSense
         )