diff --git a/src/CoolPDLP.jl b/src/CoolPDLP.jl
index cc577e9..6ad8804 100644
--- a/src/CoolPDLP.jl
+++ b/src/CoolPDLP.jl
@@ -16,7 +16,7 @@ using StableRNGs: StableRNG
 using LinearAlgebra: LinearAlgebra, Diagonal, axpby!, diag, dot, mul!, norm
 using Printf: @sprintf
 using Random: randn!
-using SparseArrays: SparseArrays, SparseMatrixCSC, AbstractSparseMatrix, findnz, nnz, nonzeros, nzrange, sparse, sprandn
+using SparseArrays: SparseArrays, SparseMatrixCSC, AbstractSparseMatrix, findnz, nnz, nonzeros, nzrange, sparse, sprandn, spzeros
 
 include("public.jl")
 
@@ -28,7 +28,7 @@ include("public.jl")
     include("utils/linalg.jl")
     include("utils/test.jl")
 
-    include("problems/milp.jl")
+    include("problems/program.jl")
     include("problems/solution.jl")
     include("problems/modify.jl")
 
@@ -52,7 +52,8 @@ include("MOI_wrapper.jl")
 
 export GPUSparseMatrixCOO, GPUSparseMatrixCSR, GPUSparseMatrixELL
 
-export MILP, nbvar, nbvar_int, nbvar_cont, nbcons, nbcons_eq, nbcons_ineq
+export AbstractProgram, LinearProgram, QuadraticProgram
+export nbvar, nbvar_int, nbvar_cont, nbcons, nbcons_eq, nbcons_ineq
 export PrimalDualSolution
 
 export preprocess, initialize, solve, solve!
diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl
index dd3006c..561a66e 100644
--- a/src/MOI_wrapper.jl
+++ b/src/MOI_wrapper.jl
@@ -71,6 +71,7 @@ const SUPPORTED_SET_TYPE{T} = Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.Less
 MOI.supports_constraint(::Optimizer{T}, ::Type{MOI.VariableIndex}, ::Type{<:SUPPORTED_SET_TYPE{T}}) where {T} = true
 MOI.supports_constraint(::Optimizer{T}, ::Type{MOI.ScalarAffineFunction{T}}, ::Type{<:SUPPORTED_SET_TYPE{T}}) where {T} = true
 MOI.supports(::Optimizer{T}, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}) where {T} = true
+MOI.supports(::Optimizer{T}, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}}) where {T} = true
 
 MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
 MOI.supports(::Optimizer, ::MOI.Silent) = true
@@ -200,14 +201,34 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
 
     c = zeros(T, n)
     obj_constant = zero(T)
-    if cache.objective.scalar_affine !== nothing
+    if !isnothing(cache.objective.scalar_affine)
         for term in cache.objective.scalar_affine.terms
             c[term.variable.value] += term.coefficient
         end
         obj_constant = cache.objective.scalar_affine.constant
     end
+
+    Q = spzeros(T, n, n)
+    if !isnothing(cache.objective.scalar_quadratic)
+        for term in cache.objective.scalar_quadratic.quadratic_terms
+            i = term.variable_1.value
+            j = term.variable_2.value
+            # MOI stores upper-triangular entries of the symmetric Q in ½xᵀQx.
+            # Each ScalarQuadraticTerm(c, xi, xj) with i ≤ j means Q[i,j] = c.
+            Q[i, j] += term.coefficient
+            if i != j
+                Q[j, i] += term.coefficient
+            end
+        end
+        for term in cache.objective.scalar_quadratic.affine_terms
+            c[term.variable.value] += term.coefficient
+        end
+        obj_constant += cache.objective.scalar_quadratic.constant
+    end
+
     if max_sense
         c .*= -one(T)
+        Q .*= -one(T)
     end
 
     dest.sets = cache.constraints.sets
@@ -216,7 +237,11 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
     lc = cache.constraints.constants.lower
     uc = cache.constraints.constants.upper
 
-    milp = MILP(; c, lv, uv, A, lc, uc)
+    milp = if nnz(Q) == 0
+        LinearProgram(; c, lv, uv, A, lc, uc)
+    else
+        QuadraticProgram(; c, lv, uv, A, Q, lc, uc)
+    end
 
     algorithm = pop!(dest.options, :algorithm, PDLP)
 
@@ -237,15 +262,26 @@ function MOI.optimize!(dest::Optimizer{T}, fcache::MOI.Utilities.UniversalFallba
 
     dest.x = Array(sol.x)
     dest.y = Array(sol.y)
-    dest.z = proj_multiplier.(c .- milp.At * dest.y, lv, uv)
+
+    grad = if milp isa QuadraticProgram
+        c .+ milp.Q * dest.x
+    else
+        c
+    end
+    dest.z = proj_multiplier.(grad .- milp.At * dest.y, lv, uv)
 
     raw_obj = objective_value(dest.x, milp)
-    raw_dual_obj = (  # lᵀ|y|⁺ - uᵀ|y|⁻ + lᵥᵀ|z|⁺ - uᵥᵀ|z|⁻
+    raw_dual_obj_lp = (  # lᵀ|y|⁺ - uᵀ|y|⁻ + lᵥᵀ|z|⁺ - uᵥᵀ|z|⁻
         sum(safeprod_left.(lc, positive_part.(dest.y)))
             - sum(safeprod_left.(uc, negative_part.(dest.y)))
             + sum(safeprod_left.(lv, positive_part.(dest.z)))
             - sum(safeprod_left.(uv, negative_part.(dest.z)))
     )
+    raw_dual_obj = if milp isa QuadraticProgram
+        raw_dual_obj_lp - dot(dest.x, milp.Q * dest.x) / 2
+    else
+        raw_dual_obj_lp
+    end
     dest.obj_value = (max_sense ? -raw_obj : raw_obj) + obj_constant
     dest.dual_obj_value = (max_sense ? -raw_dual_obj : raw_dual_obj) + obj_constant
     dest.solve_time = stats.time_elapsed
diff --git a/src/algorithms/common.jl b/src/algorithms/common.jl
index ae864bb..351ad69 100644
--- a/src/algorithms/common.jl
+++ b/src/algorithms/common.jl
@@ -152,7 +152,7 @@ Apply preconditioning, type conversion and device transfer to `milp_init` and `s
 Return a tuple `(milp, sol)`.
 """
 function preprocess(
-        milp_init_cpu::MILP,
+        milp_init_cpu::AbstractProgram,
         sol_init_cpu::PrimalDualSolution,
         algo::Algorithm,
     )
@@ -178,20 +178,20 @@ function initialize end
 """
     solve(milp, sol, algo)
     solve(milp, algo)
-    
+
 Solve the continuous relaxation of `milp` starting from solution `sol` using the algorithm defined by `algo`.
 
 Return a couple `(sol, stats)` where `sol` is the last solution and `stats` contains convergence information.
 """
 function solve(
-        milp_init_cpu::MILP,
+        milp_init_cpu::AbstractProgram,
         sol_init_cpu::PrimalDualSolution,
         algo::Algorithm
     )
     starting_time = time()
     milp, sol = preprocess(milp_init_cpu, sol_init_cpu, algo)
     state = initialize(milp, sol, algo; starting_time)
-    if nbcons(milp) == 0 && all(iszero, milp.c) # early exit for 0 obj/no cons
+    if milp isa LinearProgram && nbcons(milp) == 0 && all(iszero, milp.c) # early exit for 0 obj/no cons
         @. sol.x = proj_box(zero(eltype(milp.lv)), milp.lv, milp.uv)
         state.stats.termination_status = OPTIMAL
         return get_solution(state, milp), state.stats
@@ -201,7 +201,7 @@ function solve(
 end
 
 function solve(
-        milp_init_cpu::MILP,
+        milp_init_cpu::AbstractProgram,
         algo::Algorithm
     )
     sol_init_cpu = PrimalDualSolution(zero(milp_init_cpu.lv), zero(milp_init_cpu.lc))
@@ -217,7 +217,7 @@ function solve! end
 
 function termination_check!(
         state::AbstractState,
-        milp::MILP,
+        milp::AbstractProgram,
         algo::Algorithm
     )
     (; sol, scratch, stats) = state
@@ -230,6 +230,6 @@ function termination_check!(
     return stats.termination_status !== STILL_RUNNING
 end
 
-function get_solution(state::AbstractState, milp::MILP)
+function get_solution(state::AbstractState, milp::AbstractProgram)
     return unprecondition(state.sol, Preconditioner(milp))
 end
diff --git a/src/algorithms/pdhg.jl b/src/algorithms/pdhg.jl
index 009210a..f5db03f 100644
--- a/src/algorithms/pdhg.jl
+++ b/src/algorithms/pdhg.jl
@@ -31,6 +31,8 @@ $(TYPEDFIELDS)
     sol_last::PrimalDualSolution{T, V}
     "step sizes"
     step_sizes::StepSizes{T}
+    "step size parameters"
+    step_size_params::StepSizeParameters{T}
     "scratch space"
     scratch::Scratch{T, V}
     "convergence stats"
@@ -38,58 +40,58 @@ $(TYPEDFIELDS)
 end
 
 function initialize(
-        milp::MILP{T, V},
+        milp::AbstractProgram{T, V},
         sol::PrimalDualSolution{T, V},
         algo::Algorithm{:PDHG, T};
         starting_time::Float64
     ) where {T, V}
     sol_last = zero(sol)
-    η = fixed_stepsize(milp, algo.step_size)
-    ω = one(η)
-    step_sizes = StepSizes(; η, ω)
-    scratch = Scratch(; x = similar(sol.x), y = similar(sol.y), r = similar(sol.x))
+    η, ω, norm_A, norm_Q = init_stepsize(milp, algo.step_size)
+    step_sizes = StepSizes(; η, ω, norm_A, norm_Q)
+    step_size_params = algo.step_size
+    scratch = Scratch(sol, milp)
     stats = ConvergenceStats(T; starting_time)
-    state = PDHGState(; sol, sol_last, step_sizes, scratch, stats)
+    state = PDHGState(; sol, sol_last, step_sizes, step_size_params, scratch, stats)
     return state
 end
 
 function solve!(
         state::PDHGState,
-        milp::MILP,
+        milp::AbstractProgram,
         algo::Algorithm{:PDHG}
     )
-    prog = ProgressUnknown(desc = "PDHG iterations:", enabled = algo.generic.show_progress)
+    progress = ProgressUnknown(desc = "PDHG iterations:", enabled = algo.generic.show_progress)
     while true
         yield()
         for _ in 1:algo.generic.check_every
             step!(state, milp)
-            next!(prog; showvalues = prog_showvalues(state))
+            next!(progress; showvalues = prog_showvalues(state))
         end
         if termination_check!(state, milp, algo)
             break
         end
     end
-    finish!(prog)
+    finish!(progress)
     return state
 end
 
 function step!(
-        state::PDHGState{T, V},
-        milp::MILP{T, V},
-    ) where {T, V}
-    # switch pointers
+        state::PDHGState,
+        milp::AbstractProgram,
+    )
     state.sol, state.sol_last = state.sol_last, state.sol
-
-    (; sol, sol_last, step_sizes, scratch) = state
+    (; sol, sol_last, step_sizes, step_size_params, scratch) = state
     (; x, y) = sol_last
-    (; η, ω) = step_sizes
-    (; c, lv, uv, A, At, lc, uc) = milp
+    (; lv, uv, A, At, lc, uc) = milp
+
+    η, ω = update_step_size!(step_sizes, milp, step_size_params)
 
     τ, σ = η / ω, η * ω
 
-    # xp = proj_box.(x - τ * (c - At * y), lv, uv)
+    # xp = proj_box.(x - τ * (grad - At * y), lv, uv)
+    g = grad_x(scratch, x, milp)
     At_y = mul!(scratch.x, At, y)
-    @. sol.x = proj_box(x - τ * (c - At_y), lv, uv)
+    @. sol.x = proj_box(x - τ * (g - At_y), lv, uv)
     xdiff = @. scratch.x = 2sol.x - x
 
     # yp = y - σ * A * (2xp - x) - σ * proj_box.(inv(σ) * y - A * (2xp - x), -uc, -lc)
diff --git a/src/algorithms/pdlp.jl b/src/algorithms/pdlp.jl
index d2ba2c5..f1307cf 100644
--- a/src/algorithms/pdlp.jl
+++ b/src/algorithms/pdlp.jl
@@ -32,6 +32,8 @@ $(TYPEDFIELDS)
     sol_restart::PrimalDualSolution{T, V}
     "step sizes"
     step_sizes::StepSizes{T}
+    "step size parameters"
+    step_size_params::StepSizeParameters{T}
     "scratch space"
     scratch::Scratch{T, V}
     "iteration counter"
@@ -43,7 +45,7 @@ $(TYPEDFIELDS)
 end
 
 function initialize(
-        milp::MILP{T, V},
+        milp::AbstractProgram{T, V},
         sol::PrimalDualSolution{T, V},
         algo::Algorithm{:PDLP, T};
         starting_time::Float64
@@ -52,59 +54,59 @@ function initialize(
     sol_avg = copy(sol)
     sol_avg_last = zero(sol)
     sol_restart = copy(sol)
-    η = fixed_stepsize(milp, algo.step_size)
-    ω = primal_weight_init(milp, algo.step_size)
-    step_sizes = StepSizes(; η, ω)
-    scratch = Scratch(; x = similar(sol.x), y = similar(sol.y), r = similar(sol.x))
+    η, ω, norm_A, norm_Q = init_stepsize(milp, algo.step_size)
+    step_sizes = StepSizes(; η, ω, norm_A, norm_Q)
+    step_size_params = algo.step_size
+    scratch = Scratch(sol, milp)
     iteration = IterationCounter(0, 0, 0)
     restart_stats = RestartStats(T)
     stats = ConvergenceStats(T; starting_time)
     state = PDLPState(;
         sol, sol_last, sol_avg, sol_avg_last, sol_restart,
-        step_sizes, scratch, iteration, restart_stats, stats
+        step_sizes, step_size_params, scratch, iteration, restart_stats, stats
     )
     return state
 end
 
 function solve!(
         state::PDLPState,
-        milp::MILP,
+        milp::AbstractProgram,
         algo::Algorithm{:PDLP}
     )
-    prog = ProgressUnknown(desc = "PDLP iterations:", enabled = algo.generic.show_progress)
+    progress = ProgressUnknown(desc = "PDLP iterations:", enabled = algo.generic.show_progress)
     while true
         yield()
         for _ in 1:algo.generic.check_every
             step!(state, milp)
-            next!(prog; showvalues = prog_showvalues(state))
+            next!(progress; showvalues = prog_showvalues(state))
         end
         if termination_check!(state, milp, algo)
             break
         elseif restart_check!(state, milp, algo)
-            restart!(state, algo)
+            restart!(state, milp, algo)
         end
     end
-    finish!(prog)
+    finish!(progress)
     return state
 end
 
 function step!(
-        state::PDLPState{T, V},
-        milp::MILP{T, V},
-    ) where {T, V}
-    # switch pointers
+        state::PDLPState,
+        milp::AbstractProgram,
+    )
     state.sol, state.sol_last = state.sol_last, state.sol
-
-    (; sol, sol_last, step_sizes, scratch) = state
+    (; sol, sol_last, step_sizes, step_size_params, scratch) = state
     (; x, y) = sol_last
-    (; η, ω) = step_sizes
-    (; c, lv, uv, A, At, lc, uc) = milp
+    (; lv, uv, A, At, lc, uc) = milp
+
+    η, ω = update_step_size!(step_sizes, milp, step_size_params)
 
     τ, σ = η / ω, η * ω
 
-    # xp = proj_box.(x - τ * (c - At * y), lv, uv)
+    # xp = proj_box.(x - τ * (grad - At * y), lv, uv)
+    g = grad_x(scratch, x, milp)
     At_y = mul!(scratch.x, At, y)
-    @. sol.x = proj_box(x - τ * (c - At_y), lv, uv)
+    @. sol.x = proj_box(x - τ * (g - At_y), lv, uv)
     xdiff = @. scratch.x = 2sol.x - x
 
     # yp = y - σ * A * (2xp - x) - σ * proj_box.(inv(σ) * y - A * (2xp - x), -uc, -lc)
@@ -132,7 +134,7 @@ end
 
 function restart_check!(
         state::PDLPState,
-        milp::MILP,
+        milp::AbstractProgram,
         algo::Algorithm{:PDLP}
     )
     (;
@@ -164,10 +166,10 @@ function restart_check!(
     return should_restart(restart_stats, step_sizes, iteration, algo.restart)
 end
 
-function restart!(state::PDLPState{T}, algo::Algorithm{:PDLP}) where {T}
+function restart!(state::PDLPState{T}, milp::AbstractProgram, algo::Algorithm{:PDLP}) where {T}
     (;
         sol, sol_avg, sol_restart,
-        step_sizes, iteration, scratch, restart_stats,
+        step_sizes, step_size_params, iteration, scratch, restart_stats,
     ) = state
 
     # identify candidate
@@ -179,8 +181,9 @@ function restart!(state::PDLPState{T}, algo::Algorithm{:PDLP}) where {T}
     # update step sizes (must be done before losing previous restart)
     step_sizes.η_sum = zero(T)
     step_sizes.ω = primal_weight_update!(
-        scratch, step_sizes, sol_cand, sol_restart, algo.step_size
+        scratch, step_sizes, sol_cand, sol_restart, step_size_params
     )
+    update_step_size!(step_sizes, milp, step_size_params)
     # update solutions
     copy!(sol, sol_cand)
     zero!(sol_avg)
diff --git a/src/components/conversion.jl b/src/components/conversion.jl
index 962f666..488398e 100644
--- a/src/components/conversion.jl
+++ b/src/components/conversion.jl
@@ -35,7 +35,7 @@ function Base.show(io::IO, params::ConversionParameters{T, Ti, M}) where {T, Ti,
 end
 
 function perform_conversion(
-        milp::MILP,
+        milp::AbstractProgram,
         params::ConversionParameters{T, Ti, M},
     ) where {T, Ti, M}
     (; backend) = params
diff --git a/src/components/errors.jl b/src/components/errors.jl
index fc49be8..f33190c 100644
--- a/src/components/errors.jl
+++ b/src/components/errors.jl
@@ -62,28 +62,50 @@ function absolute(err::KKTErrors, ω::Number)
     return sqrt(ω^2 * primal^2 + inv(ω^2) * dual^2 + gap^2)
 end
 
+_half_xQx(cx, x, milp::LinearProgram{T}) where {T} = zero(T)
+_half_xQx(cx, x, milp::QuadraticProgram) = (cx - dot(milp.c, x)) / 2
+
+grad_x(scratch, x, milp::LinearProgram) = milp.c
+function grad_x(scratch, x, milp::QuadraticProgram)
+    mul!(scratch.r, milp.Q, x)
+    @. scratch.r += milp.c
+    return scratch.r
+end
+
 function kkt_errors!(
         scratch::Scratch,
         sol::PrimalDualSolution,
-        milp::MILP{T},
+        milp::AbstractProgram{T},
     ) where {T}
     (; x, y) = sol
-    (; c, lv, uv, A, At, lc, uc, D1, D2) = milp
+    (; lv, uv, A, At, lc, uc, D1, D2) = milp
+    g = grad_x(scratch, x, milp)
+
+    cx = dot(g, x)
+    rescaled_obj = @. scratch.x = inv(D2.diag) * g
+    dual_scale = one(T) + norm(rescaled_obj)
 
     A_x = mul!(scratch.y, A, x)
     At_y = mul!(scratch.x, At, y)
-    r = @. scratch.r = proj_multiplier(c - At_y, lv, uv)
+
+    h = @. scratch.r = g - At_y
+    dual_diff = @. scratch.x = inv(D2.diag) * (h - proj_multiplier(h, lv, uv))
+    dual = norm(dual_diff)
+    # h is no longer needed, reuse the scratch.r
+    r = @. scratch.r = proj_multiplier(scratch.r, lv, uv)
 
     primal_diff = @. scratch.y = inv(D1.diag) * (A_x - proj_box(A_x, lc, uc))
     primal = norm(primal_diff)
     rescaled_combined_bounds = @. scratch.y = inv(D1.diag) * combine(lc, uc)
     primal_scale = one(T) + norm(rescaled_combined_bounds)
 
-    dual_diff = @. scratch.x = inv(D2.diag) * (c - At_y - r)
-    dual = norm(dual_diff)
-    rescaled_obj = @. scratch.x = inv(D2.diag) * c
-    dual_scale = one(T) + norm(rescaled_obj)
-
+    # primal obj P = cᵀx + ½xᵀQx
+    # dual obj   D = lᵀy⁺ − uᵀy⁻ + lᵥᵀr⁺ − uᵥᵀr⁻ − ½xᵀQx
+    # gap = |P − D|
+    #     = |(cᵀx + ½xᵀQx) − (lᵀy⁺ − uᵀy⁻ + lᵥᵀr⁺ − uᵥᵀr⁻ − ½xᵀQx)|
+    #     = |cᵀx + xᵀQx + (uᵀy⁻ − lᵀy⁺) + (uᵥᵀr⁻ − lᵥᵀr⁺)|
+    #     = |gᵀx + pc_sum + pv_sum|;  g = c + Qx
+    # gap_scale = 1 + |P| + |D| = 1 + |gᵀx − ½xᵀQx| + |pc_sum + pv_sum + ½xᵀQx|
     pc = @. scratch.y = (
         safeprod_left(uc, positive_part(-y)) - safeprod_left(lc, negative_part(-y))
     )
@@ -92,10 +114,10 @@ function kkt_errors!(
     )
     pc_sum = sum(pc)
     pv_sum = sum(pv)
-    cx = dot(c, x)
 
     gap = abs(cx + pc_sum + pv_sum)
-    gap_scale = one(T) + abs(pc_sum + pv_sum) + abs(cx)
+    half_xQx = _half_xQx(cx, x, milp)
+    gap_scale = one(T) + abs(cx - half_xQx) + abs(pc_sum + pv_sum + half_xQx)
 
     err = KKTErrors(;
         primal,
diff --git a/src/components/permutation.jl b/src/components/permutation.jl
index f00a6a1..9844bd6 100644
--- a/src/components/permutation.jl
+++ b/src/components/permutation.jl
@@ -3,6 +3,7 @@ function increasing_column_order(A::SparseMatrixCSC)
     return sortperm(col_lengths)
 end
 
+permute_rows_columns(::Nothing; kwargs...) = nothing
 function permute_rows_columns(
         A::SparseMatrixCSC; perm_col::Vector{Int}, perm_row::Vector{Int}
     )
@@ -16,9 +17,9 @@ end
 """
     sort_rows_columns(milp)
 
-Return a new `MILP` where the constraint matrix has been permuted by order of increasing column and row density.
+Return a new program where the constraint matrix has been permuted by order of increasing column and row density.
 """
-function sort_rows_columns(milp::MILP)
+function sort_rows_columns(milp::AbstractProgram)
     (;
         c, lv, uv, A, At, lc, uc, D1, D2,
         int_var, var_names, dataset, name, path,
@@ -27,12 +28,15 @@ function sort_rows_columns(milp::MILP)
     perm_var = increasing_column_order(A)
     perm_cons = increasing_column_order(At)
 
-    return MILP(;
+    Q = get_Q(milp)
+    return rebuild(
+        milp;
         c = c[perm_var],
         lv = lv[perm_var],
         uv = uv[perm_var],
         A = permute_rows_columns(A; perm_col = perm_var, perm_row = perm_cons),
         At = permute_rows_columns(At; perm_col = perm_cons, perm_row = perm_var),
+        Q = permute_rows_columns(Q; perm_col = perm_var, perm_row = perm_var),
         lc = lc[perm_cons],
         uc = uc[perm_cons],
         D1 = Diagonal(diag(D1)[perm_cons]),
diff --git a/src/components/preconditioning.jl b/src/components/preconditioning.jl
index 5d81ff7..58ebe82 100644
--- a/src/components/preconditioning.jl
+++ b/src/components/preconditioning.jl
@@ -12,7 +12,7 @@ struct Preconditioner{T <: Number, V <: DenseVector{T}}
     D2::Diagonal{T, V}
 end
 
-Preconditioner(milp::MILP) = Preconditioner(milp.D1, milp.D2)
+Preconditioner(milp::AbstractProgram) = Preconditioner(milp.D1, milp.D2)
 
 function Base.:*(prec_out::Preconditioner, prec_in::Preconditioner)
     return Preconditioner(prec_out.D1 * prec_in.D1, prec_in.D2 * prec_out.D2)
@@ -58,7 +58,10 @@ function unprecondition(sol::PrimalDualSolution, prec::Preconditioner)
     return PrimalDualSolution(x, y)
 end
 
-function precondition(milp::MILP, prec::Preconditioner)
+precondition_Q(::Nothing, ::Diagonal) = nothing
+precondition_Q(Q::AbstractMatrix, D2::Diagonal) = D2 * Q * D2
+
+function precondition(milp::AbstractProgram, prec::Preconditioner)
     (;
         c, lv, uv, A, At, lc, uc,
         int_var, var_names, dataset, name, path,
@@ -71,12 +74,14 @@ function precondition(milp::MILP, prec::Preconditioner)
     A_p, At_p = cons_p.A, cons_p.At
     lc_p, uc_p = D1 * lc, D1 * uc
     new_prec = prec * Preconditioner(milp)
-    milp_p = MILP(;
+    return rebuild(
+        milp;
         c = c_p,
         lv = lv_p,
         uv = uv_p,
         A = A_p,
         At = At_p,
+        Q = precondition_Q(get_Q(milp), D2),
         lc = lc_p,
         uc = uc_p,
         D1 = new_prec.D1,
@@ -87,7 +92,6 @@ function precondition(milp::MILP, prec::Preconditioner)
         name,
         path
     )
-    return milp_p
 end
 
 # Preconditioner construction
@@ -100,10 +104,14 @@ function identity_preconditioner(cons::ConstraintMatrix{T}) where {T}
 end
 
 function diagonal_norm_preconditioner(
-        cons::ConstraintMatrix{T}; p_row::Number, p_col::Number
+        cons::ConstraintMatrix{T}, Q = nothing; p_row::Number, p_col::Number
     ) where {T}
     (; A, At) = cons
     col_norms = map(j -> column_norm(A, j, p_col), axes(A, 2))
+    if !isnothing(Q)
+        q_col_norms = map(j -> column_norm(Q, j, p_col), axes(Q, 2))
+        col_norms = combine_norms(col_norms, q_col_norms, p_col)
+    end
     row_norms = map(i -> column_norm(At, i, p_row), axes(A, 1))
     d1 = map(rn -> iszero(rn) ? one(T) : inv(sqrt(rn)), row_norms)
     d2 = map(cn -> iszero(cn) ? one(T) : inv(sqrt(cn)), col_norms)
@@ -114,14 +122,17 @@ function chambolle_pock_preconditioner(cons::ConstraintMatrix; alpha::Number)
     return diagonal_norm_preconditioner(cons; p_row = 2 - alpha, p_col = alpha)
 end
 
-function ruiz_preconditioner(cons::ConstraintMatrix; iterations::Integer)
+function ruiz_preconditioner(milp::AbstractProgram; iterations::Integer)
+    cons = ConstraintMatrix(milp.A, milp.At)
     prec = identity_preconditioner(cons)
+    Q_cur = get_Q(milp)
     for _ in 1:iterations
-        prec_next = diagonal_norm_preconditioner(cons; p_col = Inf, p_row = Inf)
+        prec_next = diagonal_norm_preconditioner(cons, Q_cur; p_col = Inf, p_row = Inf)
         cons = precondition(cons, prec_next)
+        Q_cur = precondition_Q(Q_cur, prec_next.D2)
         prec = prec_next * prec
     end
-    return prec
+    return prec, cons
 end
 
 """
@@ -143,12 +154,9 @@ function Base.show(io::IO, params::PreconditioningParameters)
     return print(io, "PreconditioningParameters: chambolle_pock_alpha=$chambolle_pock_alpha, ruiz_iter=$ruiz_iter")
 end
 
-function pdlp_preconditioner(milp::MILP, params::PreconditioningParameters)
-    (; A, At) = milp
+function pdlp_preconditioner(milp::AbstractProgram, params::PreconditioningParameters)
     (; chambolle_pock_alpha, ruiz_iter) = params
-    cons = ConstraintMatrix(A, At)
-    prec_r = ruiz_preconditioner(cons; iterations = ruiz_iter)
-    cons_r = precondition(cons, prec_r)
+    prec_r, cons_r = ruiz_preconditioner(milp; iterations = ruiz_iter)
     prec_cp = chambolle_pock_preconditioner(cons_r; alpha = chambolle_pock_alpha)
     prec = prec_r * prec_cp
     return prec
diff --git a/src/components/scratch.jl b/src/components/scratch.jl
index 834e0a3..3bfea42 100644
--- a/src/components/scratch.jl
+++ b/src/components/scratch.jl
@@ -7,4 +7,8 @@
     r::V
 end
 
-Scratch(sol::PrimalDualSolution) = Scratch(similar(sol.x), similar(sol.y), similar(sol.x))
+Scratch(sol::PrimalDualSolution) = Scratch(;
+    x = similar(sol.x), y = similar(sol.y), r = similar(sol.x),
+)
+
+Scratch(sol::PrimalDualSolution, ::AbstractProgram) = Scratch(sol)
diff --git a/src/components/step_size.jl b/src/components/step_size.jl
index 6cfa12f..ca8348d 100644
--- a/src/components/step_size.jl
+++ b/src/components/step_size.jl
@@ -19,14 +19,26 @@ function Base.show(io::IO, params::StepSizeParameters)
     return print(io, "StepSizeParameters: invnorm_scaling=$invnorm_scaling, primal_weight_damping=$primal_weight_damping, zero_tol=$zero_tol")
 end
 
-function fixed_stepsize(milp::MILP{T}, params::StepSizeParameters) where {T}
-    (; A, At) = milp
-    (; invnorm_scaling) = params
-    η = T(invnorm_scaling) * inv(spectral_norm(A, At))
-    return η
-end
+"""
+    compute_eta(norm_A, norm_Q, ω, invnorm_scaling)
 
-function primal_weight_init(milp::MILP{T}, params::StepSizeParameters) where {T}
+Compute step size `η` satisfying the Vu-Condat convergence condition:
+`η²‖A‖² + (η/(2ω))‖Q‖ ≤ 1`.
+"""
+function compute_eta(norm_A::T, norm_Q::T, ω::T, invnorm_scaling::T) where {T}
+    a = norm_A^2
+    b = norm_Q / (2ω)
+    if iszero(a) && iszero(b)
+        return one(T)
+    elseif iszero(a)
+        return invnorm_scaling * 2ω / norm_Q
+    else
+        η_max = 2 / (hypot(b, 2norm_A) + b)
+        return invnorm_scaling * η_max
+    end
+end
+primal_weight_init(::LinearProgram{T}, ::StepSizeParameters) where {T} = one(T)
+function primal_weight_init(milp::QuadraticProgram{T}, params::StepSizeParameters) where {T}
     (; c, lc, uc) = milp
     (; zero_tol) = params
     c_norm = norm(c)
@@ -39,6 +51,16 @@ function primal_weight_init(milp::MILP{T}, params::StepSizeParameters) where {T}
     end
 end
 
+function init_stepsize(milp::AbstractProgram{T}, params::StepSizeParameters) where {T}
+    (; A, At) = milp
+    Q = get_Q(milp)
+    norm_A = T(spectral_norm(A, At))
+    norm_Q = T(spectral_norm(Q, Q))
+    ω = primal_weight_init(milp, params)
+    η = compute_eta(norm_A, norm_Q, ω, T(params.invnorm_scaling))
+    return η, ω, norm_A, norm_Q
+end
+
 """
     StepSizes
 
@@ -53,6 +75,18 @@ $(TYPEDFIELDS)
     η_sum::T = zero(η)
     "primal weight"
     ω::T
+    "spectral norm of A"
+    norm_A::T = zero(η)
+    "spectral norm of Q"
+    norm_Q::T = zero(η)
+end
+
+update_step_size!(step_sizes::StepSizes, ::LinearProgram, ::StepSizeParameters) = (step_sizes.η, step_sizes.ω)
+function update_step_size!(step_sizes::StepSizes, ::QuadraticProgram, params::StepSizeParameters)
+    step_sizes.η = compute_eta(
+        step_sizes.norm_A, step_sizes.norm_Q, step_sizes.ω, params.invnorm_scaling,
+    )
+    return step_sizes.η, step_sizes.ω
 end
 
 function primal_weight_update!(
diff --git a/src/problems/milp.jl b/src/problems/milp.jl
deleted file mode 100644
index 657a226..0000000
--- a/src/problems/milp.jl
+++ /dev/null
@@ -1,218 +0,0 @@
-"""
-    MILP
-
-Represent a Mixed Integer Linear Program in "cuPDLPx form":
-
-    min cᵀx   s.t.   lv ≤ x ≤ uv
-                     lc ≤ A * x ≤ uc
-
-# Constructor
-
-    MILP(;
-        c, lv, uv, A, lc, uc,
-        [D1, D2, int_var, var_names, dataset, name, path]
-    )
-
-# Fields
-
-$(TYPEDFIELDS)
-"""
-struct MILP{
-        T <: Number,
-        V <: DenseVector{T},
-        M <: AbstractMatrix{T},
-        Vb <: DenseVector{Bool},
-    }
-    "objective vector"
-    c::V
-    "variable lower bound"
-    lv::V
-    "variable upper bound"
-    uv::V
-    "constraint matrix"
-    A::M
-    "transposed constraint matrix"
-    At::M
-    "constraint lower bound"
-    lc::V
-    "constraint upper bound"
-    uc::V
-    "left preconditioner"
-    D1::Diagonal{T, V}
-    "right preconditioner"
-    D2::Diagonal{T, V}
-    "which variables must be integers"
-    int_var::Vb
-    "variable names"
-    var_names::Vector{String}
-    "source dataset"
-    dataset::String
-    "instance name (last part of the path)"
-    name::String
-    "file path the MILP was read from"
-    path::String
-
-    function MILP(;
-            c,
-            lv,
-            uv,
-            A,
-            At = convert(typeof(A), transpose(A)),
-            lc,
-            uc,
-            D1 = Diagonal(one!(similar(lc))),
-            D2 = Diagonal(one!(similar(lv))),
-            int_var = zero!(similar(c, Bool)),
-            var_names = map(string, eachindex(c)),
-            dataset = "",
-            name = "",
-            path = ""
-        )
-        m, n = size(A)
-        if !(n == length(c) == length(lv) == length(uv) == size(D2, 1) == length(int_var) == length(var_names))
-            throw(DimensionMismatch("Variable size not consistent"))
-        elseif !(m == length(lc) == length(uc) == size(D1, 2))
-            throw(DimensionMismatch("Constraint size not consistent"))
-        end
-
-        T = Base.promote_eltype(c, lv, uv, A, At, lc, uc, D1, D2)
-        V = promote_type(typeof(c), typeof(lv), typeof(uv), typeof(lc), typeof(uc))
-        M = promote_type(typeof(A), typeof(At))
-        Vb = typeof(int_var)
-
-        if (
-                !isconcretetype(T) ||
-                    !isconcretetype(V) ||
-                    !isconcretetype(M) ||
-                    !isconcretetype(Vb)
-            )
-            throw(ArgumentError("Abstract type parameter"))
-        end
-
-        common_backend(c, lv, uv, A, At, lc, uc, D1, D2)
-
-        if isempty(name) && !isempty(path)
-            name = splitext(splitpath(path)[end])[1]
-        end
-
-        return new{T, V, M, Vb}(
-            c,
-            lv,
-            uv,
-            A,
-            At,
-            lc,
-            uc,
-            D1,
-            D2,
-            int_var,
-            var_names,
-            string(dataset),
-            string(name),
-            string(path)
-        )
-    end
-end
-
-"""
-    MILP(qps::QPSData; kwargs...)
-
-Construct a [`MILP`](@ref) from a `QPSData` object generated by [QPSReader.jl](https://github.com/JuliaSmoothOptimizers/QPSReader.jl).
-"""
-function MILP(qps::QPSData; kwargs...)
-    return MILP(;
-        c = qps.c,
-        lv = qps.lvar,
-        uv = qps.uvar,
-        A = sparse(qps.arows, qps.acols, qps.avals, length(qps.lcon), length(qps.lvar)),
-        At = sparse(qps.acols, qps.arows, qps.avals, length(qps.lvar), length(qps.lcon)),
-        lc = qps.lcon,
-        uc = qps.ucon,
-        D1 = Diagonal(ones(length(qps.lcon))),
-        D2 = Diagonal(ones(length(qps.lvar))),
-        int_var = convert(Vector{Bool}, (qps.vartypes .== VTYPE_Binary) .| (qps.vartypes .== VTYPE_Integer)),
-        var_names = qps.varnames,
-        kwargs...
-    )
-end
-
-function Base.show(io::IO, milp::MILP{T, V, M}) where {T, V, M}
-    return print(
-        io, """
-        MILP instance $(milp.name) from dataset $(milp.dataset):
-        - types:
-          - values $T
-          - vectors $V
-          - matrices $M
-        - variables: $(nbvar(milp))
-          - $(nbvar_cont(milp)) continuous
-          - $(nbvar_int(milp)) integer
-        - constraints: $(nbcons(milp))
-          - $(nbcons_ineq(milp)) inequalities
-          - $(nbcons_eq(milp)) equalities
-        - nonzeros: $(mynnz(milp.A))"""
-    )
-end
-
-KernelAbstractions.get_backend(milp::MILP) = get_backend(milp.c)
-
-"""
-    nbvar(milp)
-
-Return the number of variables in `milp`.
-"""
-nbvar(milp::MILP) = length(milp.c)
-
-"""
-    nbvar_int(milp)
-
-Return the number of integer variables in `milp`.
-"""
-nbvar_int(milp::MILP) = sum(milp.int_var)
-
-"""
-    nbvar_cont(milp)
-
-Return the number of integer variables in `milp`.
-"""
-nbvar_cont(milp::MILP) = nbvar(milp) - nbvar_int(milp)
-
-"""
-    nbcons(milp)
-
-Return the number of constraints in `milp`, not including variable bounds or integrality requirements.
-"""
-nbcons(milp::MILP) = size(milp.A, 1)
-
-"""
-    nbcons_eq(milp)
-
-Return the number of equality constraints in `milp`.
-"""
-nbcons_eq(milp::MILP) = mapreduce((l, u) -> (l == u), +, milp.lc, milp.uc)
-
-"""
-    nbcons_ineq(milp)
-
-Return the number of inequality constraints in `milp`, not including variable bounds.
-"""
-nbcons_ineq(milp::MILP) = nbcons(milp) - nbcons_eq(milp)
-
-function Base.isapprox(m1::MILP, m2::MILP; kwargs...)
-    return (
-        isapprox(m1.c, m2.c; kwargs...) &&
-            isapprox(m1.lv, m2.lv; kwargs...) &&
-            isapprox(m1.uv, m2.uv; kwargs...) &&
-            isapprox(m1.A, m2.A; kwargs...) &&
-            isapprox(m1.At, m2.At; kwargs...) &&
-            isapprox(m1.lc, m2.lc; kwargs...) &&
-            isapprox(m1.uc, m2.uc; kwargs...) &&
-            isapprox(m1.D1, m2.D1; kwargs...) &&
-            isapprox(m1.D2, m2.D2; kwargs...) &&
-            m1.int_var == m2.int_var &&
-            m1.var_names == m2.var_names &&
-            m1.dataset == m2.dataset &&
-            m1.name == m2.name &&
-            m1.path == m2.path
-    )
-end
diff --git a/src/problems/modify.jl b/src/problems/modify.jl
index 173f794..8ed93e2 100644
--- a/src/problems/modify.jl
+++ b/src/problems/modify.jl
@@ -4,22 +4,25 @@
 Change the element type of floating-point containers inside `milp` to `T`.
 """
 set_eltype(::Type{T}, A::AbstractArray{<:AbstractFloat}) where {T} = map(T, A)
+set_eltype(::Type, ::Nothing) = nothing
 
 function set_eltype(::Type{T}, sol::PrimalDualSolution) where {T}
     return PrimalDualSolution(set_eltype(T, sol.x), set_eltype(T, sol.y))
 end
 
-function set_eltype(::Type{T}, milp::MILP) where {T}
+function set_eltype(::Type{T}, milp::AbstractProgram) where {T}
     (;
         c, lv, uv, A, At, lc, uc, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
-    return MILP(;
+    return rebuild(
+        milp;
         c = set_eltype(T, c),
         lv = set_eltype(T, lv),
         uv = set_eltype(T, uv),
         A = set_eltype(T, A),
         At = set_eltype(T, At),
+        Q = set_eltype(T, get_Q(milp)),
         lc = set_eltype(T, lc),
         uc = set_eltype(T, uc),
         D1 = set_eltype(T, D1),
@@ -38,6 +41,7 @@ end
 Change the element type of integer containers inside `milp` to `Ti`.
 """
 set_indtype(::Type{Ti}, A::AbstractArray{<:Integer}) where {Ti} = map(Ti, A)
+set_indtype(::Type, ::Nothing) = nothing
 
 function set_indtype(::Type{Ti}, A::SparseMatrixCSC) where {Ti}
     return SparseMatrixCSC(
@@ -49,21 +53,20 @@ function set_indtype(::Type{Ti}, A::SparseMatrixCSC) where {Ti}
     )
 end
 
-function set_indtype(::Type{Ti}, milp::MILP) where {Ti}
+function set_indtype(::Type{Ti}, milp::AbstractProgram) where {Ti}
     (;
         c, lv, uv, A, At, lc, uc, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
-    return MILP(;
+    return rebuild(
+        milp;
         c,
         lv,
         uv,
         A = set_indtype(Ti, A),
         At = set_indtype(Ti, At),
-        lc,
-        uc,
-        D1,
-        D2,
+        Q = set_indtype(Ti, get_Q(milp)),
+        lc, uc, D1, D2,
         int_var,
         var_names,
         dataset,
@@ -77,48 +80,55 @@ end
 
 Convert the sparse matrices inside `milp` using constructor `M`.
 """
-function set_matrix_type(::Type{M}, milp::MILP) where {M}
+function set_matrix_type(::Type{M}, milp::AbstractProgram) where {M}
     (;
         c, lv, uv, A, At, lc, uc, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
-    A_M = M(A)
-    At_M = M(At)
-    backend = common_backend(A_M, At_M)
+    Q = get_Q(milp)
+    A_M = set_matrix_type(M, A)
+    At_M = set_matrix_type(M, At)
+    Q_M = set_matrix_type(M, Q)
+    backend = common_backend(A_M, At_M, Q_M)
 
-    return MILP(;
+    return rebuild(
+        milp;
         c = adapt(backend, c),
         lv = adapt(backend, lv),
         uv = adapt(backend, uv),
         A = A_M,
         At = At_M,
+        Q = Q_M,
         lc = adapt(backend, lc),
         uc = adapt(backend, uc),
         D1 = adapt(backend, D1),
         D2 = adapt(backend, D2),
         int_var = adapt(backend, int_var),
-        var_names,
-        dataset,
-        name,
-        path
+        var_names, dataset, name, path
     )
 end
 
+set_matrix_type(::Type{M}, mat) where {M} = M(mat)
+set_matrix_type(::Type, ::Nothing) = nothing
+
 function Adapt.adapt_structure(to, sol::PrimalDualSolution)
     return PrimalDualSolution(adapt(to, sol.x), adapt(to, sol.y))
 end
 
-function Adapt.adapt_structure(to, milp::MILP)
+function Adapt.adapt_structure(to, milp::AbstractProgram)
     (;
         c, lv, uv, A, At, lc, uc, D1, D2,
         int_var, var_names, dataset, name, path,
     ) = milp
-    return MILP(;
+    Q = get_Q(milp)
+    return rebuild(
+        milp;
         c = adapt(to, c),
         lv = adapt(to, lv),
         uv = adapt(to, uv),
         A = adapt(to, A),
         At = adapt(to, At),
+        Q = adapt(to, Q),
         lc = adapt(to, lc),
         uc = adapt(to, uc),
         D1 = adapt(to, D1),
diff --git a/src/problems/program.jl b/src/problems/program.jl
new file mode 100644
index 0000000..7c1d8ae
--- /dev/null
+++ b/src/problems/program.jl
@@ -0,0 +1,345 @@
+"""
+    AbstractProgram{T, V, M, Vb}
+
+Abstract supertype for optimization programs in "cuPDLPx form":
+
+    min (1/2)xᵀQx + cᵀx   s.t.   lv ≤ x ≤ uv
+                                 lc ≤ A * x ≤ uc
+
+Subtypes: [`LinearProgram`](@ref) (Q = 0) and [`QuadraticProgram`](@ref) (Q ≠ 0).
+"""
+abstract type AbstractProgram{
+    T <: Number,
+    V <: DenseVector{T},
+    M <: AbstractMatrix{T},
+    Vb <: DenseVector{Bool},
+}
+end
+
+function _validate_program_args(; c, lv, uv, A, At, lc, uc, D1, D2, int_var, Q = nothing)
+    m, n = size(A)
+    if !(n == length(c) == length(lv) == length(uv) == size(D2, 1) == length(int_var))
+        throw(DimensionMismatch("Variable size not consistent"))
+    elseif !(m == length(lc) == length(uc) == size(D1, 2))
+        throw(DimensionMismatch("Constraint size not consistent"))
+    end
+    if !isnothing(Q) && size(Q) != (n, n)
+        throw(DimensionMismatch("Q must be $n × $n, got $(size(Q))"))
+    end
+
+    eltype_args = isnothing(Q) ? (c, lv, uv, A, At, lc, uc, D1, D2) : (c, lv, uv, A, At, Q, lc, uc, D1, D2)
+    T = Base.promote_eltype(eltype_args...)
+    V = promote_type(typeof(c), typeof(lv), typeof(uv), typeof(lc), typeof(uc))
+    M = promote_type(typeof(A), typeof(At))
+    Vb = typeof(int_var)
+
+    if !isconcretetype(T) || !isconcretetype(V) || !isconcretetype(M) || !isconcretetype(Vb)
+        throw(ArgumentError("Abstract type parameter"))
+    end
+
+    return T, V, M, Vb
+end
+
+"""
+    LinearProgram
+
+Represent a Mixed Integer Linear Program in "cuPDLPx form":
+
+    min cᵀx   s.t.   lv ≤ x ≤ uv
+                     lc ≤ A * x ≤ uc
+
+# Constructor
+
+    LinearProgram(;
+        c, lv, uv, A, lc, uc,
+        [D1, D2, int_var, var_names, dataset, name, path]
+    )
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct LinearProgram{
+        T <: Number,
+        V <: DenseVector{T},
+        M <: AbstractMatrix{T},
+        Vb <: DenseVector{Bool},
+    } <: AbstractProgram{T, V, M, Vb}
+    "objective vector"
+    c::V
+    "variable lower bound"
+    lv::V
+    "variable upper bound"
+    uv::V
+    "constraint matrix"
+    A::M
+    "transposed constraint matrix"
+    At::M
+    "constraint lower bound"
+    lc::V
+    "constraint upper bound"
+    uc::V
+    "left preconditioner"
+    D1::Diagonal{T, V}
+    "right preconditioner"
+    D2::Diagonal{T, V}
+    "which variables must be integers"
+    int_var::Vb
+    "variable names"
+    var_names::Vector{String}
+    "source dataset"
+    dataset::String
+    "instance name (last part of the path)"
+    name::String
+    "file path the program was read from"
+    path::String
+
+    function LinearProgram(;
+            c,
+            lv,
+            uv,
+            A,
+            At = convert(typeof(A), transpose(A)),
+            lc,
+            uc,
+            D1 = Diagonal(one!(similar(lc))),
+            D2 = Diagonal(one!(similar(lv))),
+            int_var = zero!(similar(c, Bool)),
+            var_names = map(string, eachindex(c)),
+            dataset = "",
+            name = "",
+            path = ""
+        )
+        T, V, M, Vb = _validate_program_args(; c, lv, uv, A, At, lc, uc, D1, D2, int_var)
+
+        common_backend(c, lv, uv, A, At, lc, uc, D1, D2)
+
+        if isempty(name) && !isempty(path)
+            name = splitext(splitpath(path)[end])[1]
+        end
+
+        return new{T, V, M, Vb}(
+            c, lv, uv, A, At, lc, uc, D1, D2,
+            int_var, var_names,
+            string(dataset), string(name), string(path)
+        )
+    end
+end
+
+"""
+    QuadraticProgram
+
+Represent a Mixed Integer Quadratic Program in "cuPDLPx form":
+
+    min (1/2)xᵀQx + cᵀx   s.t.   lv ≤ x ≤ uv
+                                 lc ≤ A * x ≤ uc
+
+# Constructor
+
+    QuadraticProgram(;
+        c, lv, uv, A, Q, lc, uc,
+        [D1, D2, int_var, var_names, dataset, name, path]
+    )
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct QuadraticProgram{
+        T <: Number,
+        V <: DenseVector{T},
+        M <: AbstractMatrix{T},
+        Vb <: DenseVector{Bool},
+    } <: AbstractProgram{T, V, M, Vb}
+    "objective vector"
+    c::V
+    "variable lower bound"
+    lv::V
+    "variable upper bound"
+    uv::V
+    "constraint matrix"
+    A::M
+    "transposed constraint matrix"
+    At::M
+    "objective Hessian (symmetric positive semidefinite)"
+    Q::M
+    "constraint lower bound"
+    lc::V
+    "constraint upper bound"
+    uc::V
+    "left preconditioner"
+    D1::Diagonal{T, V}
+    "right preconditioner"
+    D2::Diagonal{T, V}
+    "which variables must be integers"
+    int_var::Vb
+    "variable names"
+    var_names::Vector{String}
+    "source dataset"
+    dataset::String
+    "instance name (last part of the path)"
+    name::String
+    "file path the program was read from"
+    path::String
+
+    function QuadraticProgram(;
+            c,
+            lv,
+            uv,
+            A,
+            At = convert(typeof(A), transpose(A)),
+            Q,
+            lc,
+            uc,
+            D1 = Diagonal(one!(similar(lc))),
+            D2 = Diagonal(one!(similar(lv))),
+            int_var = zero!(similar(c, Bool)),
+            var_names = map(string, eachindex(c)),
+            dataset = "",
+            name = "",
+            path = ""
+        )
+        T, V, M, Vb = _validate_program_args(; c, lv, uv, A, At, lc, uc, D1, D2, int_var, Q)
+
+        Q_M = convert(M, Q)
+
+        common_backend(c, lv, uv, A, At, Q_M, lc, uc, D1, D2)
+
+        if isempty(name) && !isempty(path)
+            name = splitext(splitpath(path)[end])[1]
+        end
+
+        return new{T, V, M, Vb}(
+            c, lv, uv, A, At, Q_M, lc, uc, D1, D2,
+            int_var, var_names,
+            string(dataset), string(name), string(path)
+        )
+    end
+end
+
+"""
+    LinearProgram(qps::QPSData; kwargs...)
+
+Construct a [`LinearProgram`](@ref) from a `QPSData` object
+generated by [QPSReader.jl](https://github.com/JuliaSmoothOptimizers/QPSReader.jl).
+
+See also [`QuadraticProgram(::QPSData)`](@ref).
+"""
+function LinearProgram(qps::QPSData; kwargs...)
+    return LinearProgram(;
+        c = qps.c,
+        lv = qps.lvar,
+        uv = qps.uvar,
+        A = sparse(qps.arows, qps.acols, qps.avals, length(qps.lcon), length(qps.lvar)),
+        At = sparse(qps.acols, qps.arows, qps.avals, length(qps.lvar), length(qps.lcon)),
+        lc = qps.lcon,
+        uc = qps.ucon,
+        D1 = Diagonal(ones(length(qps.lcon))),
+        D2 = Diagonal(ones(length(qps.lvar))),
+        int_var = convert(Vector{Bool}, (qps.vartypes .== VTYPE_Binary) .| (qps.vartypes .== VTYPE_Integer)),
+        var_names = qps.varnames,
+        kwargs...
+    )
+end
+
+"""
+    QuadraticProgram(qps::QPSData; kwargs...)
+
+Construct a [`QuadraticProgram`](@ref) from a `QPSData` object
+generated by [QPSReader.jl](https://github.com/JuliaSmoothOptimizers/QPSReader.jl).
+
+See also [`LinearProgram(::QPSData)`](@ref).
+"""
+function QuadraticProgram(qps::QPSData; kwargs...)
+    n = length(qps.lvar)
+    Q_upper = sparse(qps.qrows, qps.qcols, qps.qvals, n, n)
+    Q = Q_upper + Q_upper' - Diagonal(diag(Q_upper))
+    return QuadraticProgram(;
+        c = qps.c,
+        lv = qps.lvar,
+        uv = qps.uvar,
+        A = sparse(qps.arows, qps.acols, qps.avals, length(qps.lcon), length(qps.lvar)),
+        At = sparse(qps.acols, qps.arows, qps.avals, length(qps.lvar), length(qps.lcon)),
+        Q,
+        lc = qps.lcon,
+        uc = qps.ucon,
+        D1 = Diagonal(ones(length(qps.lcon))),
+        D2 = Diagonal(ones(length(qps.lvar))),
+        int_var = convert(Vector{Bool}, (qps.vartypes .== VTYPE_Binary) .| (qps.vartypes .== VTYPE_Integer)),
+        var_names = qps.varnames,
+        kwargs...
+    )
+end
+
+# Show methods
+
+function Base.show(io::IO, milp::AbstractProgram{T, V, M}) where {T, V, M}
+    type_name = nameof(typeof(milp))
+    quadratic_line = milp isa QuadraticProgram ? "\n        - quadratic: true" : ""
+    return print(
+        io, """
+        $type_name instance $(milp.name) from dataset $(milp.dataset):
+        - types:
+          - values $T
+          - vectors $V
+          - matrices $M
+        - variables: $(nbvar(milp))
+          - $(nbvar_cont(milp)) continuous
+          - $(nbvar_int(milp)) integer
+        - constraints: $(nbcons(milp))
+          - $(nbcons_ineq(milp)) inequalities
+          - $(nbcons_eq(milp)) equalities
+        - nonzeros: $(mynnz(milp.A))$quadratic_line"""
+    )
+end
+
+KernelAbstractions.get_backend(milp::AbstractProgram) = get_backend(milp.c)
+
+
+get_Q(::LinearProgram) = nothing
+get_Q(milp::QuadraticProgram) = milp.Q
+
+rebuild(::LinearProgram; Q::Nothing = nothing, kwargs...) = LinearProgram(; kwargs...)
+rebuild(::QuadraticProgram; kwargs...) = QuadraticProgram(; kwargs...)
+
+"""
+    nbvar(milp)
+
+Return the number of variables in `milp`.
+"""
+nbvar(milp::AbstractProgram) = length(milp.c)
+
+"""
+    nbvar_int(milp)
+
+Return the number of integer variables in `milp`.
+"""
+nbvar_int(milp::AbstractProgram) = sum(milp.int_var)
+
+"""
+    nbvar_cont(milp)
+
+Return the number of continuous variables in `milp`.
+"""
+nbvar_cont(milp::AbstractProgram) = nbvar(milp) - nbvar_int(milp)
+
+"""
+    nbcons(milp)
+
+Return the number of constraints in `milp`, not including variable bounds or integrality requirements.
+"""
+nbcons(milp::AbstractProgram) = size(milp.A, 1)
+
+"""
+    nbcons_eq(milp)
+
+Return the number of equality constraints in `milp`.
+"""
+nbcons_eq(milp::AbstractProgram) = mapreduce((l, u) -> (l == u), +, milp.lc, milp.uc)
+
+"""
+    nbcons_ineq(milp)
+
+Return the number of inequality constraints in `milp`, not including variable bounds.
+"""
+nbcons_ineq(milp::AbstractProgram) = nbcons(milp) - nbcons_eq(milp)
diff --git a/src/problems/solution.jl b/src/problems/solution.jl
index 587ec3d..8b07650 100644
--- a/src/problems/solution.jl
+++ b/src/problems/solution.jl
@@ -10,7 +10,7 @@ Check whether solution vector `x` is feasible for `milp`.
 - `verbose`: whether to display warnings
 """
 function is_feasible(
-        x, milp::MILP;
+        x, milp::AbstractProgram;
         cons_tol = 1.0e-6, int_tol = 1.0e-5, verbose::Bool = true
     )
     (; lv, uv, A, lc, uc, int_var) = milp
@@ -35,9 +35,13 @@ end
 """
     objective_value(x, milp)
 
-Compute the value of the linear objective of `milp` at solution vector `x`.
+Compute the objective value of `milp` at solution vector `x`.
 """
-objective_value(x, milp::MILP) = dot(x, milp.c)
+objective_value(x, milp::LinearProgram) = dot(x, milp.c)
+
+function objective_value(x, milp::QuadraticProgram)
+    return dot(x, milp.c) + dot(x, milp.Q * x) / 2
+end
 
 """
     PrimalDualSolution
diff --git a/src/utils/device.jl b/src/utils/device.jl
index 3c02471..a0797b9 100644
--- a/src/utils/device.jl
+++ b/src/utils/device.jl
@@ -4,9 +4,12 @@
 Return the common GPU backend of several arguments, if it exists, and throw an error otherwise.
 """
 function common_backend(args::Vararg{Any, N}) where {N}
-    backends = map(get_backend, args)
-    if !all(==(backends[1]), backends)
+    backends = map(_get_backend, args)
+    if !all(x -> isnothing(x) || (x == backends[1]), backends)
         throw(ArgumentError("There are several different backends among the arguments: $(unique(backends))"))
     end
     return backends[1]
 end
+
+_get_backend(::Nothing) = nothing
+_get_backend(x) = get_backend(x)
diff --git a/src/utils/linalg.jl b/src/utils/linalg.jl
index 9207d00..95ac38c 100644
--- a/src/utils/linalg.jl
+++ b/src/utils/linalg.jl
@@ -94,9 +94,11 @@ function spectral_norm(
     λ, _ = powm!(KᵀK, x0; kwargs...)
     return sqrt(λ)
 end
+spectral_norm(::Nothing, ::Nothing; kwargs...) = 0
 
 column_norm(A::AbstractMatrix, j::Integer, p) = norm(view(A, :, j), p)
 column_norm(A::SparseMatrixCSC, j::Integer, p) = norm(view(nonzeros(A), nzrange(A, j)), p)
+combine_norms(a, b, p) = isinf(p) ? max.(a, b) : (a .^ p .+ b .^ p) .^ (1 / p)
 
 mynnz(A::AbstractSparseMatrix) = nnz(A)
 mynnz(A::AbstractMatrix) = prod(size(A))
diff --git a/src/utils/test.jl b/src/utils/test.jl
index f5cdb8b..09c5291 100644
--- a/src/utils/test.jl
+++ b/src/utils/test.jl
@@ -46,5 +46,5 @@ function random_milp_and_sol(m::Int, n::Int, p::Float64)
     int_var = rand(Bool, length(c))
     x = proj_box.(randn(n), lv, uv)
     y = proj_multiplier.(randn(m), lc, uc)
-    return MILP(; c, lv, uv, A, lc, uc, int_var), PrimalDualSolution(x, y)
+    return LinearProgram(; c, lv, uv, A, lc, uc, int_var), PrimalDualSolution(x, y)
 end
diff --git a/test/algorithms/all.jl b/test/algorithms/all.jl
index 6ddd7bb..008d882 100644
--- a/test/algorithms/all.jl
+++ b/test/algorithms/all.jl
@@ -10,7 +10,7 @@ using SparseArrays
 using Test
 
 netlib_milps = map(list_instances(Netlib)) do name
-    MILP(read_instance(dataset, name)[1]; dataset, name)
+    LinearProgram(read_instance(dataset, name)[1]; dataset, name)
 end
 sort!(netlib_milps, by = milp -> nbvar(milp))
 small_names = filter(map(milp -> milp.name, netlib_milps[1:3])) do name
@@ -22,7 +22,7 @@ function test_optimizer(
         obj_rtol::Float64 = 1.0e-2, cons_tol::Float64 = 1.0e-2, int_tol::Float64 = Inf,
     )
     qps, path = read_instance(dataset, name)
-    milp = MILP(qps; dataset, path)
+    milp = LinearProgram(qps; dataset, path)
 
     jump_model = JuMP.read_from_file(path; format = MOI.FileFormats.FORMAT_MPS)
     JuMP.set_optimizer(jump_model, HiGHS.Optimizer)
diff --git a/test/algorithms/qp.jl b/test/algorithms/qp.jl
new file mode 100644
index 0000000..1b11fb3
--- /dev/null
+++ b/test/algorithms/qp.jl
@@ -0,0 +1,175 @@
+using CoolPDLP
+using Random
+using JLArrays
+using KernelAbstractions
+using LinearAlgebra
+using SparseArrays
+using Test
+
+# min (1/2)(x₁² + x₂²) s.t. x₁ + x₂ = 1, x ≥ 0
+function simple_equality_qp()
+    n = 2
+    c = zeros(n)
+    Q = sparse(1.0I, n, n)  # Identity Hessian
+    A = sparse([1.0 1.0])
+    lv = zeros(n)
+    uv = fill(Inf, n)
+    lc = [1.0]
+    uc = [1.0]
+    return QuadraticProgram(; c, Q, lv, uv, A, lc, uc)
+end
+
+# min (1/2)(2x₁² + x₂²) + x₁ + 2x₂ s.t. x₁ + x₂ ≥ 1, x ≥ 0
+function linear_quadratic_qp()
+    n = 2
+    c = [1.0, 2.0]
+    Q = sparse([2.0 0.0; 0.0 1.0])
+    A = sparse([1.0 1.0])
+    lv = zeros(n)
+    uv = fill(Inf, n)
+    lc = [1.0]
+    uc = [Inf]
+    return QuadraticProgram(; c, Q, lv, uv, A, lc, uc)
+end
+
+function random_qp(n, m; seed = 42)
+    rng = Random.MersenneTwister(seed)
+    B = randn(rng, n, n)
+    Q_dense = B' * B + 0.1I
+    Q = sparse(Q_dense)
+    c = randn(rng, n)
+    A = sprandn(rng, m, n, 0.3)
+    x_feas = abs.(randn(rng, n))
+    b = A * x_feas
+    lv = zeros(n)
+    uv = fill(Inf, n)
+    lc = b
+    uc = b
+    return QuadraticProgram(; c, Q, lv, uv, A, lc, uc)
+end
+
+function high_linear_weight_qp()
+    n = 2
+    c = [100.0, 100.0]
+    Q = sparse([10.0 0.0; 0.0 10.0])
+    A = sparse([1.0 1.0])
+    lv = zeros(n)
+    uv = fill(Inf, n)
+    lc = [1.0]
+    uc = [1.0]
+    return QuadraticProgram(; c, Q, lv, uv, A, lc, uc)
+end
+
+@testset "QP - PDLP" begin
+    @testset "Simple equality QP" begin
+        milp = simple_equality_qp()
+        @test milp isa QuadraticProgram
+        algo = PDLP(; termination_reltol = 1.0e-6, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test sol.x[1] ≈ 1 / 2 atol = 1.0e-4
+        @test sol.x[2] ≈ 1 / 2 atol = 1.0e-4
+        @test objective_value(sol.x, milp) ≈ 0.25 atol = 1.0e-4
+    end
+
+    @testset "Linear + quadratic QP" begin
+        milp = linear_quadratic_qp()
+        algo = PDLP(; termination_reltol = 1.0e-6, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test sol.x[1] ≈ 2 / 3 atol = 1.0e-3
+        @test sol.x[2] ≈ 1 / 3 atol = 1.0e-3
+    end
+
+    @testset "Random QP" begin
+        milp = random_qp(20, 5)
+        algo = PDLP(; termination_reltol = 1.0e-4, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test is_feasible(sol.x, milp; cons_tol = 1.0e-2)
+    end
+end
+
+@testset "QP - PDHG" begin
+    @testset "Simple equality QP" begin
+        milp = simple_equality_qp()
+        algo = PDHG(; termination_reltol = 1.0e-6, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test sol.x[1] ≈ 1 / 2 atol = 1.0e-4
+        @test sol.x[2] ≈ 1 / 2 atol = 1.0e-4
+    end
+
+    @testset "Linear + quadratic QP" begin
+        milp = linear_quadratic_qp()
+        algo = PDHG(; termination_reltol = 1.0e-6, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test sol.x[1] ≈ 2 / 3 atol = 1.0e-3
+        @test sol.x[2] ≈ 1 / 3 atol = 1.0e-3
+    end
+
+    @testset "High linear weight QP" begin
+        milp = high_linear_weight_qp()
+        algo = PDHG(; termination_reltol = 1.0e-6, max_kkt_passes = 10^6, show_progress = false)
+        sol, stats = solve(milp, algo)
+        @test stats.termination_status == CoolPDLP.OPTIMAL
+        @test sol.x[1] ≈ 1 / 2 atol = 1.0e-4
+        @test sol.x[2] ≈ 1 / 2 atol = 1.0e-4
+    end
+end
+
+@testset "QP step-size numerics" begin
+    @testset "compute_eta" begin
+        T = Float64
+        # norm_A and norm_Q nonzero: η = invnorm_scaling * 2 / (hypot(b, 2*norm_A) + b), b = norm_Q/(2ω)
+        let norm_A = 2.0, norm_Q = 4.0, ω = 1.5, s = 0.9
+            b = norm_Q / (2ω)
+            @test CoolPDLP.compute_eta(norm_A, norm_Q, ω, s) ≈ s * 2 / (hypot(b, 2norm_A) + b)
+        end
+
+        # norm_A = 0 (Q-only): η = invnorm_scaling * 2ω / norm_Q
+        let norm_Q = 4.0, ω = 2.0, s = 0.5
+            @test CoolPDLP.compute_eta(zero(T), norm_Q, ω, s) ≈ s * 2ω / norm_Q
+        end
+
+        # norm_Q = 0 (LP): η = invnorm_scaling / norm_A
+        let norm_A = 2.0, s = 0.9
+            @test CoolPDLP.compute_eta(norm_A, zero(T), 1.0, s) ≈ s / norm_A
+        end
+
+        # both zero: fallback to 1.0 regardless of ω and invnorm_scaling
+        let
+            @test CoolPDLP.compute_eta(zero(T), zero(T), 2.5, 0.7) == 1.0
+        end
+    end
+
+    @testset "primal_weight_init" begin
+        lp, _ = CoolPDLP.random_milp_and_sol(5, 10, 0.4)
+        params = CoolPDLP.StepSizeParameters(; invnorm_scaling = 0.9, primal_weight_damping = 0.5, zero_tol = 1.0e-10)
+        @test CoolPDLP.primal_weight_init(lp, params) == 1.0
+
+        qp = linear_quadratic_qp()
+        ω_qp = CoolPDLP.primal_weight_init(qp, params)
+        expected_ω = norm(qp.c) / norm(CoolPDLP.combine.(qp.lc, qp.uc))
+        @test ω_qp ≈ expected_ω
+    end
+
+    @testset "update_step_size! LP is no-op" begin
+        lp, _ = CoolPDLP.random_milp_and_sol(5, 10, 0.4)
+        params = CoolPDLP.StepSizeParameters(; invnorm_scaling = 0.9, primal_weight_damping = 0.5, zero_tol = 1.0e-10)
+        step_sizes = CoolPDLP.StepSizes(; η = 0.1, ω = 2.0, norm_A = 1.0, norm_Q = 0.0)
+        η0, ω0 = step_sizes.η, step_sizes.ω
+        CoolPDLP.update_step_size!(step_sizes, lp, params)
+        @test step_sizes.η == η0
+        @test step_sizes.ω == ω0
+    end
+
+    @testset "update_step_size! QP recomputes η" begin
+        qp = simple_equality_qp()
+        params = CoolPDLP.StepSizeParameters(; invnorm_scaling = 0.9, primal_weight_damping = 0.5, zero_tol = 1.0e-10)
+        step_sizes = CoolPDLP.StepSizes(; η = 999.0, ω = 1.0, norm_A = 1.0, norm_Q = 1.0)
+        CoolPDLP.update_step_size!(step_sizes, qp, params)
+        @test step_sizes.η ≈ CoolPDLP.compute_eta(1.0, 1.0, 1.0, 0.9)
+    end
+end
diff --git a/test/components/conversion.jl b/test/components/conversion.jl
index 333fd5e..1038d46 100644
--- a/test/components/conversion.jl
+++ b/test/components/conversion.jl
@@ -6,7 +6,7 @@ milp, sol = CoolPDLP.random_milp_and_sol(10, 20, 0.4)
 params = CoolPDLP.ConversionParameters(Float32, Int32, GPUSparseMatrixCSR; backend = JLBackend())
 
 milp_gpu = CoolPDLP.perform_conversion(milp, params)
-@test milp_gpu isa MILP{Float32, JLVector{Float32}, GPUSparseMatrixCSR{Float32, Int32, JLVector{Float32}, JLVector{Int32}}}
+@test milp_gpu isa LinearProgram{Float32, JLVector{Float32}, GPUSparseMatrixCSR{Float32, Int32, JLVector{Float32}, JLVector{Int32}}}
 
 sol_gpu = CoolPDLP.perform_conversion(sol, params)
 @test sol_gpu isa PrimalDualSolution{Float32, JLVector{Float32}}
diff --git a/test/components/errors.jl b/test/components/errors.jl
index 5d416a4..a8e2282 100644
--- a/test/components/errors.jl
+++ b/test/components/errors.jl
@@ -1,5 +1,6 @@
 using CoolPDLP
 using LinearAlgebra
+using SparseArrays
 using Test
 
 function p(y, l, u)
@@ -36,3 +37,47 @@ end
 @testset "Invariance by preconditioning" begin
     @test err_p ≈ err
 end
+
+n, m = 20, 10
+c = randn(n)
+lv = zeros(n)
+uv = fill(Inf, n)
+A = sprand(m, n, 0.4)
+At = sparse(A')
+lc = randn(m)
+uc = lc + rand(m)
+H = sprand(n, n, 0.3)
+Q = Matrix(H' * H)
+
+qp = QuadraticProgram(; c, lv, uv, A, Q = sparse(Q), lc, uc)
+x = abs.(randn(n))
+y = randn(m)
+sol_qp = PrimalDualSolution(x, y)
+scratch_qp = CoolPDLP.Scratch(sol_qp)
+
+err_qp = CoolPDLP.kkt_errors!(scratch_qp, sol_qp, qp)
+
+g = c + Q * x
+r_qp = CoolPDLP.proj_multiplier.(g - At * y, lv, uv)
+half_xQx = dot(x, Q * x) / 2
+cx = dot(g, x)  # = cᵀx + xᵀQx
+pc_sum = p(-y, lc, uc)
+pv_sum = p(-r_qp, lv, uv)
+
+@testset "Correct QP KKT errors" begin
+    @test err_qp.primal ≈ norm(A * x - CoolPDLP.proj_box.(A * x, lc, uc))
+    @test err_qp.dual ≈ norm(g - At * y - r_qp)
+    @test err_qp.gap ≈ abs(cx + pc_sum + pv_sum)
+    @test err_qp.dual_scale ≈ 1 + norm(g)
+    @test err_qp.gap_scale ≈ 1 + abs(cx - half_xQx) + abs(pc_sum + pv_sum + half_xQx)
+end
+
+@testset "QP KKT invariance by preconditioning" begin
+    params_qp = CoolPDLP.PreconditioningParameters(; chambolle_pock_alpha = 1.0, ruiz_iter = 5)
+    prec_qp = CoolPDLP.pdlp_preconditioner(qp, params_qp)
+    qp_p = CoolPDLP.precondition(qp, prec_qp)
+    sol_qp_p = CoolPDLP.precondition(sol_qp, prec_qp)
+    scratch_qp_p = CoolPDLP.Scratch(sol_qp_p)
+    err_qp_p = CoolPDLP.kkt_errors!(scratch_qp_p, sol_qp_p, qp_p)
+    @test err_qp_p ≈ err_qp
+end
diff --git a/test/components/permutation.jl b/test/components/permutation.jl
index 6d993da..8bb2ee0 100644
--- a/test/components/permutation.jl
+++ b/test/components/permutation.jl
@@ -1,7 +1,12 @@
 using CoolPDLP
+using LinearAlgebra
 using SparseArrays
 using Test
 
+@testset "permute_rows_columns(nothing)" begin
+    @test isnothing(CoolPDLP.permute_rows_columns(nothing; perm_col = [1, 2], perm_row = [1, 2]))
+end
+
 @testset "Sort columns" begin
     for m in (10, 20, 30), n in (10, 20, 30), p in (0.01, 0.1, 0.3)
         A = sprand(n, n, p)
@@ -13,3 +18,34 @@ using Test
         @test issorted(map(row -> count(!iszero, row), eachrow(A_sorted)))
     end
 end
+
+@testset "sort_rows_columns(LinearProgram)" begin
+    n, m = 15, 8
+    A = sprand(m, n, 0.3)
+    lp = LinearProgram(;
+        c = rand(n), lv = zeros(n), uv = ones(n),
+        A, lc = zeros(m), uc = ones(m),
+    )
+    lp_sorted = CoolPDLP.sort_rows_columns(lp)
+    @test lp_sorted isa LinearProgram
+    @test nbvar(lp_sorted) == nbvar(lp)
+    @test nbcons(lp_sorted) == nbcons(lp)
+    @test lp_sorted.A != lp.A
+end
+
+@testset "sort_rows_columns(QuadraticProgram)" begin
+    n, m = 15, 8
+    A = sprand(m, n, 0.3)
+    H = sprand(n, n, 0.3)
+    Q = H' * H
+    qp = QuadraticProgram(;
+        c = rand(n), lv = zeros(n), uv = ones(n),
+        A, Q, lc = zeros(m), uc = ones(m),
+    )
+    qp_sorted = CoolPDLP.sort_rows_columns(qp)
+    @test qp_sorted isa QuadraticProgram
+    @test nbvar(qp_sorted) == nbvar(qp)
+    @test nbcons(qp_sorted) == nbcons(qp)
+    @test qp_sorted.A != qp.A
+    @test qp_sorted.Q != qp.Q
+end
diff --git a/test/components/preconditioning.jl b/test/components/preconditioning.jl
index 8df8374..12f25cd 100644
--- a/test/components/preconditioning.jl
+++ b/test/components/preconditioning.jl
@@ -4,6 +4,35 @@ using SparseArrays
 using Random: Xoshiro
 using Test
 
+function _isapprox(m1::LinearProgram, m2::LinearProgram; kwargs...)
+    return (
+        isapprox(m1.c, m2.c; kwargs...) &&
+            isapprox(m1.lv, m2.lv; kwargs...) &&
+            isapprox(m1.uv, m2.uv; kwargs...) &&
+            isapprox(m1.A, m2.A; kwargs...) &&
+            isapprox(m1.At, m2.At; kwargs...) &&
+            isapprox(m1.lc, m2.lc; kwargs...) &&
+            isapprox(m1.uc, m2.uc; kwargs...) &&
+            isapprox(m1.D1, m2.D1; kwargs...) &&
+            isapprox(m1.D2, m2.D2; kwargs...)
+    )
+end
+
+function _isapprox(m1::QuadraticProgram, m2::QuadraticProgram; kwargs...)
+    return (
+        isapprox(m1.c, m2.c; kwargs...) &&
+            isapprox(m1.lv, m2.lv; kwargs...) &&
+            isapprox(m1.uv, m2.uv; kwargs...) &&
+            isapprox(m1.A, m2.A; kwargs...) &&
+            isapprox(m1.At, m2.At; kwargs...) &&
+            isapprox(m1.Q, m2.Q; kwargs...) &&
+            isapprox(m1.lc, m2.lc; kwargs...) &&
+            isapprox(m1.uc, m2.uc; kwargs...) &&
+            isapprox(m1.D1, m2.D1; kwargs...) &&
+            isapprox(m1.D2, m2.D2; kwargs...)
+    )
+end
+
 @testset "Composition" begin
     A = sprand(10, 20, 0.4)
     cons = CoolPDLP.ConstraintMatrix(A, sparse(transpose(A)))
@@ -34,14 +63,14 @@ end
         @test id_prec.D2 == I
     end
     @testset "Ruiz" begin
-        ruiz_prec = CoolPDLP.ruiz_preconditioner(cons; iterations = 10000)
-        cons_p = CoolPDLP.precondition(cons, ruiz_prec)
+        milp = LinearProgram(; c = zeros(20), lv = zeros(20), uv = ones(20), A, lc = zeros(10), uc = ones(10))
+        ruiz_prec, cons_p = CoolPDLP.ruiz_preconditioner(milp; iterations = 10000)
         @test all(≈(1; rtol = 1.0e-2), map(col -> norm(col, Inf), eachcol(cons_p.A)))
         @test all(≈(1; rtol = 1.0e-2), map(col -> norm(col, Inf), eachcol(cons_p.At)))
     end
 end
 
-@testset "Effect on MILP" begin
+@testset "Effect on LinearProgram" begin
     m, n, p = 10, 20, 0.4
 
     c = rand(n)
@@ -50,7 +79,7 @@ end
     A = sprand(m, n, p)
     lc = rand(m)
     uc = lc + rand(m)
-    milp = MILP(; c, lv, uv, A, lc, uc)
+    milp = LinearProgram(; c, lv, uv, A, lc, uc)
     x = randn(n)
     y = randn(m)
     sol = PrimalDualSolution(x, y)
@@ -60,8 +89,8 @@ end
 
     milp_p = CoolPDLP.precondition(milp, prec)
     milp_unp = CoolPDLP.precondition(milp_p, inv(prec))
-    @test isapprox(milp, milp_unp)
-    @test !isapprox(milp, milp_p)
+    @test _isapprox(milp, milp_unp)
+    @test !_isapprox(milp, milp_p)
 
     sol_p = CoolPDLP.precondition(sol, prec)
     sol_unp = CoolPDLP.unprecondition(sol_p, prec)
@@ -73,3 +102,40 @@ end
     @test CoolPDLP.proj_box.(sol.x, milp.lv, milp.uv) ≈ prec.D2 * CoolPDLP.proj_box.(sol_p.x, milp_p.lv, milp_p.uv)
     @test CoolPDLP.proj_box.(milp.A * sol.x, milp.lc, milp.uc) ≈ prec.D1 \ CoolPDLP.proj_box.(milp_p.A * sol_p.x, milp_p.lc, milp_p.uc)
 end
+
+@testset "Effect on QuadraticProgram" begin
+    m, n, p = 10, 20, 0.4
+
+    c = rand(n)
+    lv = rand(n)
+    uv = lv + rand(n)
+    A = sprand(m, n, p)
+    Q = let H = sprand(n, n, p)
+        H' * H
+    end
+    lc = rand(m)
+    uc = lc + rand(m)
+    qp = QuadraticProgram(; c, lv, uv, A, Q, lc, uc)
+    x = randn(n)
+    y = randn(m)
+    sol = PrimalDualSolution(x, y)
+
+    params = CoolPDLP.PreconditioningParameters(; chambolle_pock_alpha = 1, ruiz_iter = 10)
+    prec = CoolPDLP.pdlp_preconditioner(qp, params)
+
+    qp_p = CoolPDLP.precondition(qp, prec)
+    qp_unp = CoolPDLP.precondition(qp_p, inv(prec))
+    @test _isapprox(qp, qp_unp)
+    @test !_isapprox(qp, qp_p)
+    @test !(qp.Q ≈ qp_p.Q)
+
+    sol_p = CoolPDLP.precondition(sol, prec)
+    sol_unp = CoolPDLP.unprecondition(sol_p, prec)
+    @test isapprox(sol, sol_unp)
+    @test !isapprox(sol, sol_p)
+
+    @test objective_value(sol.x, qp) ≈ objective_value(sol_p.x, qp_p)
+    @test dot(sol.y, qp.A, sol.x) ≈ dot(sol_p.y, qp_p.A, sol_p.x)
+    @test CoolPDLP.proj_box.(sol.x, qp.lv, qp.uv) ≈ prec.D2 * CoolPDLP.proj_box.(sol_p.x, qp_p.lv, qp_p.uv)
+    @test CoolPDLP.proj_box.(qp.A * sol.x, qp.lc, qp.uc) ≈ prec.D1 \ CoolPDLP.proj_box.(qp_p.A * sol_p.x, qp_p.lc, qp_p.uc)
+end
diff --git a/test/components/termination.jl b/test/components/termination.jl
index dc49dbf..7bcc8f8 100644
--- a/test/components/termination.jl
+++ b/test/components/termination.jl
@@ -3,7 +3,7 @@
     lc, A, uc = [1.0], sparse([1.0 1.0]), [Inf]
     lv, uv = [0.0, 0.0], [Inf, Inf]
 
-    milp = CoolPDLP.MILP(; c, lv, uv, A, lc, uc)
+    milp = CoolPDLP.LinearProgram(; c, lv, uv, A, lc, uc)
     algo = CoolPDLP.PDLP()
     sol, stats = CoolPDLP.solve(milp, algo)
     @test stats.termination_status == CoolPDLP.OPTIMAL
diff --git a/test/moi.jl b/test/moi.jl
index f58a17d..03d2e9e 100644
--- a/test/moi.jl
+++ b/test/moi.jl
@@ -83,6 +83,19 @@ end
     @test JuMP.value(x) isa Float32
 end
 
+@testset "QuadraticProgram via JuMP" begin
+    # min (1/2)(x₁ + x₂)² s.t. x₁ + x₂ = 1, x ≥ 0
+    model = JuMP.Model(CoolPDLP.Optimizer)
+    JuMP.set_silent(model)
+    JuMP.@variable(model, x[1:2] >= 0)
+    JuMP.@constraint(model, x[1] + x[2] == 1)
+    JuMP.@objective(model, Min, 0.5 * (x[1] + x[2])^2)
+    JuMP.optimize!(model)
+    @test JuMP.termination_status(model) == MOI.OPTIMAL
+    @test JuMP.objective_value(model) ≈ 0.5 atol = 1.0e-3
+    @test JuMP.value(x[1]) + JuMP.value(x[2]) ≈ 1.0 atol = 1.0e-3
+end
+
 if CUDA.functional()
     @info "Running CUDA tests"
     CUDA.versioninfo()
diff --git a/test/problems/modify.jl b/test/problems/modify.jl
index 8cdfcea..a5b2691 100644
--- a/test/problems/modify.jl
+++ b/test/problems/modify.jl
@@ -7,11 +7,11 @@ milp, sol = CoolPDLP.random_milp_and_sol(10, 20, 0.4)
 
 @testset "Set types" begin
     milp_f32 = CoolPDLP.set_eltype(Float32, milp)
-    @test milp_f32 isa MILP{Float32, Vector{Float32}, SparseMatrixCSC{Float32, Int}}
+    @test milp_f32 isa LinearProgram{Float32, Vector{Float32}, SparseMatrixCSC{Float32, Int}}
     milp_i32 = CoolPDLP.set_indtype(Int32, milp)
-    @test milp_i32 isa MILP{Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int32}}
+    @test milp_i32 isa LinearProgram{Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int32}}
     milp_dense = CoolPDLP.set_matrix_type(Matrix, milp)
-    @test milp_dense isa MILP{Float64, Vector{Float64}, Matrix{Float64}}
+    @test milp_dense isa LinearProgram{Float64, Vector{Float64}, Matrix{Float64}}
 
     sol_f32 = CoolPDLP.set_eltype(Float32, sol)
     @test sol_f32 isa PrimalDualSolution{Float32, Vector{Float32}}
@@ -19,9 +19,9 @@ end
 
 @testset "Change backend" begin
     milp_flexible = CoolPDLP.set_matrix_type(GPUSparseMatrixCSR, milp)
-    @test milp_flexible isa MILP{Float64, Vector{Float64}, GPUSparseMatrixCSR{Float64, Int, Vector{Float64}, Vector{Int}}, Vector{Bool}}
+    @test milp_flexible isa LinearProgram{Float64, Vector{Float64}, GPUSparseMatrixCSR{Float64, Int, Vector{Float64}, Vector{Int}}, Vector{Bool}}
     milp_gpu = adapt(JLBackend(), milp_flexible)
-    @test milp_gpu isa MILP{Float64, JLVector{Float64}, GPUSparseMatrixCSR{Float64, Int, JLVector{Float64}, JLVector{Int}}, JLVector{Bool}}
+    @test milp_gpu isa LinearProgram{Float64, JLVector{Float64}, GPUSparseMatrixCSR{Float64, Int, JLVector{Float64}, JLVector{Int}}, JLVector{Bool}}
 
     sol_gpu = adapt(JLBackend(), sol)
     @test sol_gpu isa PrimalDualSolution{Float64, JLVector{Float64}}
diff --git a/test/problems/milp.jl b/test/problems/program.jl
similarity index 51%
rename from test/problems/milp.jl
rename to test/problems/program.jl
index 2066528..c230223 100644
--- a/test/problems/milp.jl
+++ b/test/problems/program.jl
@@ -10,33 +10,33 @@ using Test
     milp, _ = CoolPDLP.random_milp_and_sol(100, 200, 0.4)
     (; c, lv, uv, A, At, lc, uc, D1, D2, int_var) = milp
 
-    @test_nowarn MILP(;
+    @test_nowarn LinearProgram(;
         c, lv, uv, A, At, lc, uc, D1, D2,
     )
     # Type issues
-    @test_throws ArgumentError MILP(;
+    @test_throws ArgumentError LinearProgram(;
         c = Vector{Any}(c), lv, uv, A, At, lc, uc,
     )
-    @test_throws ArgumentError MILP(;
+    @test_throws ArgumentError LinearProgram(;
         c = jl(c), lv, uv, A, At, lc, uc,
     )
     # Dimension issues
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c = lc, lv, uv, A, At, lc, uc,
     )
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c, lv = lc, uv, A, At, lc, uc,
     )
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c, lv, uv, A = At, At, lc, uc,
     )
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c, lv, uv, A, At, lc = lv, uc,
     )
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c, lv, uv, A, At, lc, uc, D1 = D2, D2,
     )
-    @test_throws DimensionMismatch MILP(;
+    @test_throws DimensionMismatch LinearProgram(;
         c, lv, uv, A, At, lc, uc, int_var = vcat(int_var, false)
     )
 end
@@ -60,7 +60,7 @@ end
     netlib = list_instances(Netlib)
     @testset for name in netlib[randperm(length(netlib))[1:20]]
         qps, path = read_instance(Netlib, name)
-        milp = MILP(qps; path, name, dataset = "Netlib")
+        milp = LinearProgram(qps; path, name, dataset = "Netlib")
         if name in ["agg", "blend", "dfl001", "forplan", "gfrd-pnc", "sierra"]
             @test_skip JuMP.read_from_file(path; format = MOI.FileFormats.FORMAT_MPS)
         else
@@ -74,13 +74,49 @@ end;
 
 @testset "Show" begin
     qps, path = read_instance(Netlib, "seba")
-    milp = MILP(qps; path, name = "seba")
-    @test startswith(string(milp), "MILP instance seba")
+    milp = LinearProgram(qps; path, name = "seba")
+    @test startswith(string(milp), "LinearProgram instance seba")
 end
 
-@testset "Approx" begin
-    netlib = list_instances(Netlib)
-    qps, path = read_instance(Netlib, netlib[1])
-    milp = MILP(qps; path, dataset = "Netlib")
-    @test milp ≈ milp
+@testset "QuadraticProgram Checks" begin
+    n, m = 20, 10
+    c = rand(n)
+    lv = rand(n)
+    uv = lv + rand(n)
+    A = sprand(m, n, 0.4)
+    H = sprand(n, n, 0.3)
+    Q = H' * H
+    lc = rand(m)
+    uc = lc + rand(m)
+
+    @test_nowarn QuadraticProgram(; c, lv, uv, A, Q, lc, uc)
+    @test_throws DimensionMismatch QuadraticProgram(; c, lv, uv, A, Q = Q[1:(end - 1), :], lc, uc)
+    @test_throws DimensionMismatch QuadraticProgram(; c, lv, uv, A, Q = Q[:, 1:(end - 1)], lc, uc)
+    @test_throws ArgumentError QuadraticProgram(; c = Vector{Any}(c), lv, uv, A, Q, lc, uc)
+end
+
+@testset "QuadraticProgram Show" begin
+    n, m = 5, 3
+    c = zeros(n)
+    Q = sparse(1.0I, n, n)
+    A = spzeros(m, n)
+    qp = QuadraticProgram(; c, Q, lv = zeros(n), uv = ones(n), A, lc = zeros(m), uc = ones(m))
+    s = string(qp)
+    @test startswith(s, "QuadraticProgram")
+    @test contains(s, "quadratic: true")  # FIXME: redundant with type name
+end
+
+@testset "get_Q" begin
+    n, m = 10, 5
+    c = rand(n)
+    lv = zeros(n)
+    uv = ones(n)
+    A = sprand(m, n, 0.4)
+    lc = zeros(m)
+    uc = ones(m)
+    lp = LinearProgram(; c, lv, uv, A, lc, uc)
+    @test isnothing(CoolPDLP.get_Q(lp))
+    Q = sparse(1.0I, n, n)
+    qp = QuadraticProgram(; c, lv, uv, A, Q, lc, uc)
+    @test CoolPDLP.get_Q(qp) === qp.Q
 end
diff --git a/test/problems/solution.jl b/test/problems/solution.jl
index f3b509a..ab96788 100644
--- a/test/problems/solution.jl
+++ b/test/problems/solution.jl
@@ -5,7 +5,7 @@ using LinearAlgebra
 using JuMP: JuMP, MOI
 using Test
 
-@testset "Cube MILP" begin
+@testset "Cube LinearProgram" begin
     c = [1.0, 2.0]
     lv = zeros(2)
     uv = 2 .* ones(2)
@@ -14,7 +14,7 @@ using Test
     uc = [1.0]
     int_var = [true, false]
 
-    milp = MILP(; c, lv, uv, A, lc, uc, int_var)
+    milp = LinearProgram(; c, lv, uv, A, lc, uc, int_var)
     @test is_feasible([1.0, 0.0], milp)
     @test @test_warn "Integrality not satisfied" !is_feasible([0.5, 0.5], milp)
     @test @test_warn "Constraints not satisfied" !is_feasible([0.0, 0.0], milp)
@@ -25,7 +25,7 @@ end
 @testset "Comparison with JuMP" begin
     name = "afiro"
     qps, path = read_instance(Netlib, name)
-    milp = MILP(qps; path, name, dataset = "Netlib")
+    milp = LinearProgram(qps; path, name, dataset = "Netlib")
 
     jump_model = JuMP.read_from_file(milp.path; format = MOI.FileFormats.FORMAT_MPS)
     JuMP.set_optimizer(jump_model, HiGHS.Optimizer)
diff --git a/test/tutorial.jl b/test/tutorial.jl
index 664eea3..5e5f898 100644
--- a/test/tutorial.jl
+++ b/test/tutorial.jl
@@ -7,18 +7,18 @@ using JuMP: JuMP, MOI
 using MathOptBenchmarkInstances: Netlib, list_instances, read_instance
 using Test  #src
 
-# ## Creating a MILP
+# ## Creating a LinearProgram
 
-# You can use [QPSReader.jl](https://github.com/JuliaSmoothOptimizers/QPSReader.jl) to read a MILP from a local MPS file, or [MathOptBenchmarkInstances.jl](https://github.com/JuliaDecisionFocusedLearning/MathOptBenchmarkInstances.jl) to automatically download standard benchmark sets (which we do here).
+# You can use [QPSReader.jl](https://github.com/JuliaSmoothOptimizers/QPSReader.jl) to read a LinearProgram from a local MPS file, or [MathOptBenchmarkInstances.jl](https://github.com/JuliaDecisionFocusedLearning/MathOptBenchmarkInstances.jl) to automatically download standard benchmark sets (which we do here).
 
 dataset = Netlib
 list = list_instances(dataset)
 name = list[4]
 qps, path = read_instance(dataset, name);
 
-# A [`MILP`](@ref) object can be constructed from there:
+# A [`LinearProgram`](@ref) object can be constructed from there:
 
-milp = MILP(qps; dataset, name, path)
+milp = LinearProgram(qps; dataset, name, path)
 
 # Its attributes can be queried:
 
@@ -30,9 +30,9 @@ nbcons(milp)
 
 # Note that manual construction is also an option if you provide the constraints, variable bounds and objectives as arrays.
 
-# ## Solving a MILP
+# ## Solving a LinearProgram
 
-# You can use the PDLP algortithm to solve the continuous relaxation of a MILP.
+# You can use the PDLP algortithm to solve the continuous relaxation of a LinearProgram.
 # The first thing to do is define parameters inside a [`PDLP`](@ref) struct.
 
 algo = PDLP(;
diff --git a/test/utils/device.jl b/test/utils/device.jl
index 4496969..771bcae 100644
--- a/test/utils/device.jl
+++ b/test/utils/device.jl
@@ -7,4 +7,6 @@ using Test
     @test CoolPDLP.common_backend(rand(2), rand(4)) == CPU()
     @test CoolPDLP.common_backend(jl(rand(2)), jl(rand(4)), jl(rand(6))) == JLBackend()
     @test_throws ArgumentError CoolPDLP.common_backend(rand(2), jl(rand(4)), jl(rand(6)))
+    @test CoolPDLP.common_backend(rand(2), nothing, rand(4)) == CPU()
+    @test CoolPDLP.common_backend(nothing, nothing) == nothing
 end
diff --git a/test/utils/linalg.jl b/test/utils/linalg.jl
index 9c497fb..4d8bd0e 100644
--- a/test/utils/linalg.jl
+++ b/test/utils/linalg.jl
@@ -1,5 +1,6 @@
 using CoolPDLP
 using LinearAlgebra
+using Random: Xoshiro
 using SparseArrays
 using Test
 
@@ -67,3 +68,15 @@ end
         end
     end
 end
+
+@testset "spectral_norm(nothing, nothing)" begin
+    @test CoolPDLP.spectral_norm(nothing, nothing) == 0
+end
+
+@testset "combine_norms" begin
+    a = [1.0, 2.0, 3.0]
+    b = [4.0, 5.0, 6.0]
+    @test CoolPDLP.combine_norms(a, b, 1) ≈ a + b
+    @test CoolPDLP.combine_norms(a, b, 2) ≈ sqrt.(a .^ 2 + b .^ 2)
+    @test CoolPDLP.combine_norms(a, b, Inf) == max.(a, b)
+end