From 3d93ebc38f33692e49bec9fbe1202d84f1dc040f Mon Sep 17 00:00:00 2001 From: mtanneau <9593025+mtanneau@users.noreply.github.com> Date: Sun, 27 Sep 2020 23:08:42 -0400 Subject: [PATCH] Standardize naming conventions (#63) --- docs/src/manual/linear_systems.md | 20 +- docs/src/reference/kkt_solvers.md | 14 +- docs/src/reference/parameters.md | 2 +- src/Interfaces/MOI/MOI_wrapper.jl | 36 +-- src/Interfaces/MOI/attributes.jl | 12 +- src/Interfaces/MOI/constraints.jl | 320 +++++++++++++------------- src/Interfaces/MOI/objective.jl | 36 +-- src/Interfaces/MOI/variables.jl | 8 +- src/Interfaces/tulip_julia_api.jl | 38 +-- src/KKT/KKT.jl | 22 +- src/KKT/cholmod.jl | 40 ++-- src/KKT/krylov.jl | 144 ++++++------ src/KKT/lapack.jl | 76 +++--- src/KKT/ldlfact.jl | 52 ++--- src/KKT/test.jl | 24 +- src/LinearAlgebra/LinearAlgebra.jl | 10 +- src/Presolve/Presolve.jl | 148 ++++++------ src/Presolve/dominated_column.jl | 52 ++--- src/Presolve/empty_column.jl | 48 ++-- src/Presolve/empty_row.jl | 42 ++-- src/Presolve/fixed_variable.jl | 12 +- src/Presolve/forcing_row.jl | 44 ++-- src/Presolve/free_column_singleton.jl | 30 +-- src/Presolve/row_singleton.jl | 20 +- src/model.jl | 20 +- src/parameters.jl | 24 +- src/problemData.jl | 134 +++++------ src/solution.jl | 28 +-- test/KKT/krylov.jl | 6 +- test/KKT/lapack.jl | 8 +- test/KKT/ldlfact.jl | 8 +- test/Presolve/empty_column.jl | 18 +- test/Presolve/empty_row.jl | 60 ++--- test/examples.jl | 12 +- 34 files changed, 784 insertions(+), 784 deletions(-) diff --git a/docs/src/manual/linear_systems.md b/docs/src/manual/linear_systems.md index ad766cf6..8b2b3905 100644 --- a/docs/src/manual/linear_systems.md +++ b/docs/src/manual/linear_systems.md @@ -49,14 +49,14 @@ Here is a list of currently supported linear solvers: | Linear solver type | System | Backend | Method | |:-------------------|:-------|:--------|:-------| -| [`Dense_SymPosDef`](@ref) | Normal equations | Dense / LAPACK | Cholesky -| [`Cholmod_SymQuasDef`](@ref) | Augmented system | CHOLMOD | LDLᵀ -| [`Cholmod_SymPosDef`](@ref) | Normal equations | CHOLMOD | Cholesky -| [`LDLFact_SymQuasDef`](@ref) | Augmented system | [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl) | LDLᵀ -| [`KrylovSPDSolver`](@ref) | Normal equations | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov -| [`KrylovSIDSolver`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov -| [`KrylovSQDSolver`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov - -[^1]: [`KrylovSIDSolver`](@ref)s view the augmented system as a symmetric indefinite system, - while [`KrylovSQDSolver`](@ref)s exploit its 2x2 structure and quasi-definite property. +| [`DenseSPD`](@ref) | Normal equations | Dense / LAPACK | Cholesky +| [`CholmodSQD`](@ref) | Augmented system | CHOLMOD | LDLᵀ +| [`CholmodSPD`](@ref) | Normal equations | CHOLMOD | Cholesky +| [`LDLFactSQD`](@ref) | Augmented system | [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl) | LDLᵀ +| [`KrylovSPD`](@ref) | Normal equations | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov +| [`KrylovSID`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov +| [`KrylovSQD`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov + +[^1]: [`KrylovSID`](@ref)s view the augmented system as a symmetric indefinite system, + while [`KrylovSQD`](@ref)s exploit its 2x2 structure and quasi-definite property. See the reference documentation for more details. \ No newline at end of file diff --git a/docs/src/reference/kkt_solvers.md b/docs/src/reference/kkt_solvers.md index d8377449..e1c0d3da 100644 --- a/docs/src/reference/kkt_solvers.md +++ b/docs/src/reference/kkt_solvers.md @@ -44,7 +44,7 @@ CurrentModule = Tulip.KKT ## Dense/LAPACK ```@docs -Dense_SymPosDef +DenseSPD ``` ## CHOLMOD @@ -54,17 +54,17 @@ CholmodSolver ``` ```@docs -Cholmod_SymQuasDef +CholmodSQD ``` ```@docs -Cholmod_SymPosDef +CholmodSPD ``` ## LDLFactorizations ```@docs -LDLFact_SymQuasDef +LDLFactSQD ``` ## Krylov @@ -75,13 +75,13 @@ LDLFact_SymQuasDef ```@docs -KrylovSPDSolver +KrylovSPD ``` ```@docs -KrylovSIDSolver +KrylovSID ``` ```@docs -KrylovSQDSolver +KrylovSQD ``` \ No newline at end of file diff --git a/docs/src/reference/parameters.md b/docs/src/reference/parameters.md index ba11daff..fd3ebdd0 100644 --- a/docs/src/reference/parameters.md +++ b/docs/src/reference/parameters.md @@ -42,7 +42,7 @@ Numerical tolerances for the interior-point algorithm. | Parameter | Description | Default | |:----------|:------------|:--------| | `MatrixOptions` | See [`TLA.MatrixOptions`](@ref) | `SparseMatrixCSC` -| `KKTOptions` | See [`KKT.SolverOptions`](@ref) | [`KKT.Cholmod_SymQuasDef`](@ref) for `Float64`, [`KKT.LDLFact_SymQuasDef`](@ref) otherwise | +| `KKTOptions` | See [`KKT.SolverOptions`](@ref) | [`KKT.CholmodSQD`](@ref) for `Float64`, [`KKT.LDLFactSQD`](@ref) otherwise | | `BarrierPRegMin` | Minimum value of primal regularization | ``\sqrt{\epsilon}`` | | `BarrierDregMin` | Minimum value of dual regularization | ``\sqrt{\epsilon}`` diff --git a/src/Interfaces/MOI/MOI_wrapper.jl b/src/Interfaces/MOI/MOI_wrapper.jl index 705e26a4..aee1c0d8 100644 --- a/src/Interfaces/MOI/MOI_wrapper.jl +++ b/src/Interfaces/MOI/MOI_wrapper.jl @@ -53,17 +53,17 @@ end _bounds(s) """ -_bounds(s::MOI.EqualTo{Tv}) where{Tv} = s.value, s.value -_bounds(s::MOI.LessThan{Tv}) where{Tv} = Tv(-Inf), s.upper -_bounds(s::MOI.GreaterThan{Tv}) where{Tv} = s.lower, Tv(Inf) -_bounds(s::MOI.Interval{Tv}) where{Tv} = s.lower, s.upper +_bounds(s::MOI.EqualTo{T}) where{T} = s.value, s.value +_bounds(s::MOI.LessThan{T}) where{T} = T(-Inf), s.upper +_bounds(s::MOI.GreaterThan{T}) where{T} = s.lower, T(Inf) +_bounds(s::MOI.Interval{T}) where{T} = s.lower, s.upper -const SCALAR_SETS{Tv} = Union{ - MOI.LessThan{Tv}, - MOI.GreaterThan{Tv}, - MOI.EqualTo{Tv}, - MOI.Interval{Tv} -} where{Tv} +const SCALAR_SETS{T} = Union{ + MOI.LessThan{T}, + MOI.GreaterThan{T}, + MOI.EqualTo{T}, + MOI.Interval{T} +} where{T} @@ -76,12 +76,12 @@ const SCALAR_SETS{Tv} = Union{ # ============================================================================== """ - Optimizer{Tv} + Optimizer{T} Wrapper for MOI. """ -mutable struct Optimizer{Tv} <: MOI.AbstractOptimizer - inner::Model{Tv} +mutable struct Optimizer{T} <: MOI.AbstractOptimizer + inner::Model{T} is_feas::Bool # Model is feasibility problem if true @@ -91,7 +91,7 @@ mutable struct Optimizer{Tv} <: MOI.AbstractOptimizer var_indices_moi::Vector{MOI.VariableIndex} var_indices::Dict{MOI.VariableIndex, Int} con_indices_moi::Vector{MOI.ConstraintIndex} - con_indices::Dict{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, <:SCALAR_SETS{Tv}}, Int} + con_indices::Dict{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, <:SCALAR_SETS{T}}, Int} # Variable and constraint names name2var::Dict{String, MOI.VariableIndex} @@ -104,14 +104,14 @@ mutable struct Optimizer{Tv} <: MOI.AbstractOptimizer # Keep track of bound constraints var2bndtype::Dict{MOI.VariableIndex, Set{Type{<:MOI.AbstractScalarSet}}} - function Optimizer{Tv}(;kwargs...) where{Tv} - m = new{Tv}( - Model{Tv}(), false, + function Optimizer{T}(;kwargs...) where{T} + m = new{T}( + Model{T}(), false, # Variable and constraint counters 0, 0, # Index mapping MOI.VariableIndex[], Dict{MOI.VariableIndex, Int}(), - MOI.ConstraintIndex[], Dict{MOI.ConstraintIndex{MOI.ScalarAffineFunction, <:SCALAR_SETS{Tv}}, Int}(), + MOI.ConstraintIndex[], Dict{MOI.ConstraintIndex{MOI.ScalarAffineFunction, <:SCALAR_SETS{T}}, Int}(), # Name -> index mapping Dict{String, MOI.VariableIndex}(), Dict{String, MOI.ConstraintIndex}(), Dict{MOI.ConstraintIndex, String}(), # Variable bounds tracking diff --git a/src/Interfaces/MOI/attributes.jl b/src/Interfaces/MOI/attributes.jl index c2811da3..89febcb7 100644 --- a/src/Interfaces/MOI/attributes.jl +++ b/src/Interfaces/MOI/attributes.jl @@ -111,8 +111,8 @@ MOI.get(m::Optimizer, ::MOI.NumberOfVariables) = m.inner.pbdata.nvar # ObjectiveFunctionType # MOI.get( - m::Optimizer{Tv}, ::MOI.ObjectiveFunctionType -) where{Tv} = MOI.ScalarAffineFunction{Tv} + m::Optimizer{T}, ::MOI.ObjectiveFunctionType +) where{T} = MOI.ScalarAffineFunction{T} # # ObjectiveSense @@ -140,7 +140,7 @@ end # # ObjectiveValue # -function MOI.get(m::Optimizer{Tv}, attr::MOI.ObjectiveValue) where{Tv} +function MOI.get(m::Optimizer{T}, attr::MOI.ObjectiveValue) where{T} MOI.check_result_index_bounds(m, attr) return get_attribute(m.inner, ObjectiveValue()) end @@ -148,7 +148,7 @@ end # # DualObjectiveValue # -function MOI.get(m::Optimizer{Tv}, attr::MOI.DualObjectiveValue) where{Tv} +function MOI.get(m::Optimizer{T}, attr::MOI.DualObjectiveValue) where{T} MOI.check_result_index_bounds(m, attr) return get_attribute(m.inner, DualObjectiveValue()) end @@ -161,11 +161,11 @@ MOI.get(m::Optimizer, ::MOI.RawSolver) = m.inner # # RelativeGap # -function MOI.get(m::Optimizer{Tv}, ::MOI.RelativeGap) where{Tv} +function MOI.get(m::Optimizer{T}, ::MOI.RelativeGap) where{T} # TODO: dispatch a function call on m.inner zp = m.inner.solver.primal_bound_scaled zd = m.inner.solver.dual_bound_scaled - return (abs(zp - zd) / (Tv(1 // 10^6)) + abs(zd)) + return (abs(zp - zd) / (T(1 // 10^6)) + abs(zd)) end # diff --git a/src/Interfaces/MOI/constraints.jl b/src/Interfaces/MOI/constraints.jl index 51176122..06c94717 100644 --- a/src/Interfaces/MOI/constraints.jl +++ b/src/Interfaces/MOI/constraints.jl @@ -22,23 +22,23 @@ MOI.supports(::Optimizer, ::A, ::Type{<:MOI.ConstraintIndex}) where{A<:SUPPORTED # Variable bounds function MOI.supports_constraint( - ::Optimizer{Tv}, ::Type{MOI.SingleVariable}, ::Type{S} -) where {Tv, S<:SCALAR_SETS{Tv}} + ::Optimizer{T}, ::Type{MOI.SingleVariable}, ::Type{S} +) where {T, S<:SCALAR_SETS{T}} return true end # Linear constraints function MOI.supports_constraint( - ::Optimizer{Tv}, ::Type{MOI.ScalarAffineFunction{Tv}}, ::Type{S} -) where {Tv, S<:SCALAR_SETS{Tv}} + ::Optimizer{T}, ::Type{MOI.ScalarAffineFunction{T}}, ::Type{S} +) where {T, S<:SCALAR_SETS{T}} return true end function MOI.is_valid( - m::Optimizer{Tv}, + m::Optimizer{T}, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S <:SCALAR_SETS{Tv}} +) where{T, S <:SCALAR_SETS{T}} v = MOI.VariableIndex(c.value) MOI.is_valid(m, v) || return false res = S ∈ m.var2bndtype[v] @@ -46,9 +46,9 @@ function MOI.is_valid( end function MOI.is_valid( - m::Optimizer{Tv}, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} return haskey(m.con_indices, c) end @@ -59,21 +59,21 @@ end # TODO: make it clear that only finite bounds can be given in input. # To relax variable bounds, one should delete the associated bound constraint. function MOI.add_constraint( - m::Optimizer{Tv}, + m::Optimizer{T}, f::MOI.SingleVariable, - s::MOI.LessThan{Tv} -) where{Tv} + s::MOI.LessThan{T} +) where{T} # Check that variable exists v = f.variable MOI.throw_if_not_valid(m, v) # Check if upper bound already exists - if MOI.LessThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Tv}, MOI.LessThan{Tv}}(v.value)) - elseif MOI.EqualTo{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Tv}, MOI.LessThan{Tv}}(v.value)) - elseif MOI.Interval{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.Interval{Tv}, MOI.LessThan{Tv}}(v.value)) + if MOI.LessThan{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.LessThan{T}, MOI.LessThan{T}}(v.value)) + elseif MOI.EqualTo{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{T}, MOI.LessThan{T}}(v.value)) + elseif MOI.Interval{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.Interval{T}, MOI.LessThan{T}}(v.value)) end # Update inner model @@ -81,27 +81,27 @@ function MOI.add_constraint( set_attribute(m.inner, VariableUpperBound(), j, s.upper) # Update bound tracking - push!(m.var2bndtype[v], MOI.LessThan{Tv}) + push!(m.var2bndtype[v], MOI.LessThan{T}) - return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Tv}}(v.value) + return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{T}}(v.value) end function MOI.add_constraint( - m::Optimizer{Tv}, + m::Optimizer{T}, f::MOI.SingleVariable, - s::MOI.GreaterThan{Tv} -) where{Tv} + s::MOI.GreaterThan{T} +) where{T} # Check that variable exists v = f.variable MOI.throw_if_not_valid(m, v) # Check if lower bound already exists - if MOI.GreaterThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Tv}, MOI.GreaterThan{Tv}}(v.value)) - elseif MOI.EqualTo{Tv} ∈ m.var2bndtype[v] - throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Tv}, MOI.GreaterThan{Tv}}(v.value)) - elseif MOI.Interval{Tv} ∈ m.var2bndtype[v] - throw(MOI.LowerBoundAlreadySet{MOI.Interval{Tv}, MOI.GreaterThan{Tv}}(v.value)) + if MOI.GreaterThan{T} ∈ m.var2bndtype[v] + throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{T}, MOI.GreaterThan{T}}(v.value)) + elseif MOI.EqualTo{T} ∈ m.var2bndtype[v] + throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{T}, MOI.GreaterThan{T}}(v.value)) + elseif MOI.Interval{T} ∈ m.var2bndtype[v] + throw(MOI.LowerBoundAlreadySet{MOI.Interval{T}, MOI.GreaterThan{T}}(v.value)) end # Update inner model @@ -109,29 +109,29 @@ function MOI.add_constraint( set_attribute(m.inner, VariableLowerBound(), j, s.lower) # Update upper-bound - push!(m.var2bndtype[v], MOI.GreaterThan{Tv}) + push!(m.var2bndtype[v], MOI.GreaterThan{T}) - return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Tv}}(v.value) + return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{T}}(v.value) end function MOI.add_constraint( - m::Optimizer{Tv}, + m::Optimizer{T}, f::MOI.SingleVariable, - s::MOI.EqualTo{Tv} -) where{Tv} + s::MOI.EqualTo{T} +) where{T} # Check that variable exists v = f.variable MOI.throw_if_not_valid(m, v) # Check if a bound already exists - if MOI.LessThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Tv}, MOI.EqualTo{Tv}}(v.value)) - elseif MOI.GreaterThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Tv}, MOI.EqualTo{Tv}}(v.value)) - elseif MOI.EqualTo{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Tv}, MOI.EqualTo{Tv}}(v.value)) - elseif MOI.Interval{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.Interval{Tv}, MOI.EqualTo{Tv}}(v.value)) + if MOI.LessThan{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.LessThan{T}, MOI.EqualTo{T}}(v.value)) + elseif MOI.GreaterThan{T} ∈ m.var2bndtype[v] + throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{T}, MOI.EqualTo{T}}(v.value)) + elseif MOI.EqualTo{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{T}, MOI.EqualTo{T}}(v.value)) + elseif MOI.Interval{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.Interval{T}, MOI.EqualTo{T}}(v.value)) end # Update inner model @@ -140,29 +140,29 @@ function MOI.add_constraint( set_attribute(m.inner, VariableUpperBound(), j, s.value) # Update bound tracking - push!(m.var2bndtype[v], MOI.EqualTo{Tv}) + push!(m.var2bndtype[v], MOI.EqualTo{T}) - return MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Tv}}(v.value) + return MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{T}}(v.value) end function MOI.add_constraint( - m::Optimizer{Tv}, + m::Optimizer{T}, f::MOI.SingleVariable, - s::MOI.Interval{Tv} -) where{Tv} + s::MOI.Interval{T} +) where{T} # Check that variable exists v = f.variable MOI.throw_if_not_valid(m, v) # Check if a bound already exists - if MOI.LessThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Tv}, MOI.Interval{Tv}}(v.value)) - elseif MOI.GreaterThan{Tv} ∈ m.var2bndtype[v] - throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Tv}, MOI.Interval{Tv}}(v.value)) - elseif MOI.EqualTo{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Tv}, MOI.Interval{Tv}}(v.value)) - elseif MOI.Interval{Tv} ∈ m.var2bndtype[v] - throw(MOI.UpperBoundAlreadySet{MOI.Interval{Tv}, MOI.Interval{Tv}}(v.value)) + if MOI.LessThan{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.LessThan{T}, MOI.Interval{T}}(v.value)) + elseif MOI.GreaterThan{T} ∈ m.var2bndtype[v] + throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{T}, MOI.Interval{T}}(v.value)) + elseif MOI.EqualTo{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{T}, MOI.Interval{T}}(v.value)) + elseif MOI.Interval{T} ∈ m.var2bndtype[v] + throw(MOI.UpperBoundAlreadySet{MOI.Interval{T}, MOI.Interval{T}}(v.value)) end # Update variable bounds @@ -171,20 +171,20 @@ function MOI.add_constraint( set_attribute(m.inner, VariableUpperBound(), j, s.upper) # Update bound tracking - push!(m.var2bndtype[v], MOI.Interval{Tv}) + push!(m.var2bndtype[v], MOI.Interval{T}) - return MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Tv}}(v.value) + return MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{T}}(v.value) end # General linear constraints function MOI.add_constraint( - m::Optimizer{Tv}, - f::MOI.ScalarAffineFunction{Tv}, - s::SCALAR_SETS{Tv} -) where{Tv} + m::Optimizer{T}, + f::MOI.ScalarAffineFunction{T}, + s::SCALAR_SETS{T} +) where{T} # Check that constant term is zero if !iszero(f.constant) - throw(MOI.ScalarFunctionConstantNotZero{Tv, typeof(f), typeof(s)}(f.constant)) + throw(MOI.ScalarFunctionConstantNotZero{T, typeof(f), typeof(s)}(f.constant)) end # Convert to canonical form @@ -193,7 +193,7 @@ function MOI.add_constraint( # Extract row nz = length(fc.terms) rind = Vector{Int}(undef, nz) - rval = Vector{Tv}(undef, nz) + rval = Vector{T}(undef, nz) lb, ub = _bounds(s) for (k, t) in enumerate(fc.terms) rind[k] = m.var_indices[t.variable_index] @@ -218,9 +218,9 @@ end # 3. Delete constraints # ============================================= function MOI.delete( - m::Optimizer{Tv}, + m::Optimizer{T}, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} # Sanity check MOI.throw_if_not_valid(m, c) @@ -228,16 +228,16 @@ function MOI.delete( # Update inner model j = m.var_indices[v] - if S == MOI.LessThan{Tv} + if S == MOI.LessThan{T} # Remove upper-bound - set_attribute(m.inner, VariableUpperBound(), j, Tv(Inf)) - elseif S == MOI.GreaterThan{Tv} + set_attribute(m.inner, VariableUpperBound(), j, T(Inf)) + elseif S == MOI.GreaterThan{T} # Remove lower bound - set_attribute(m.inner, VariableLowerBound(), j, Tv(-Inf)) + set_attribute(m.inner, VariableLowerBound(), j, T(-Inf)) else # Set variable to free - set_attribute(m.inner, VariableLowerBound(), j, Tv(-Inf)) - set_attribute(m.inner, VariableUpperBound(), j, Tv(Inf)) + set_attribute(m.inner, VariableLowerBound(), j, T(-Inf)) + set_attribute(m.inner, VariableUpperBound(), j, T(Inf)) end # Update name tracking @@ -252,9 +252,9 @@ function MOI.delete( end function MOI.delete( - m::Optimizer{Tv}, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Update inner model @@ -279,10 +279,10 @@ end # 4. Modify constraints # ============================================= function MOI.modify( - m::Optimizer{Tv}, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S}, - chg::MOI.ScalarCoefficientChange{Tv} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}, + chg::MOI.ScalarCoefficientChange{T} +) where{T, S<:SCALAR_SETS{T}} MOI.is_valid(m, c) || throw(MOI.InvalidIndex(c)) MOI.is_valid(m, chg.variable) || throw(MOI.InvalidIndex(chg.variable)) @@ -303,9 +303,9 @@ end # ListOfConstraintIndices # function MOI.get( - m::Optimizer{Tv}, + m::Optimizer{T}, ::MOI.ListOfConstraintIndices{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} indices = MOI.ConstraintIndex{MOI.SingleVariable, S}[] for (var, bounds_set) in m.var2bndtype @@ -315,13 +315,13 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, - ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, + ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} return [ cidx for cidx in keys(m.con_indices) if isa(cidx, - MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} + MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} ) ] end @@ -330,9 +330,9 @@ end # NumberOfConstraints # function MOI.get( - m::Optimizer{Tv}, + m::Optimizer{T}, ::MOI.NumberOfConstraints{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} ncon = 0 for (v, bound_sets) in m.var2bndtype ncon += S ∈ bound_sets @@ -341,13 +341,13 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, - ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, + ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} ncon = 0 for cidx in keys(m.con_indices) - ncon += isa(cidx, MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S}) + ncon += isa(cidx, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}) end return ncon @@ -358,18 +358,18 @@ end # function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintName, + m::Optimizer{T}, ::MOI.ConstraintName, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) return get(m.bnd2name, c, "") end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintName, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, ::MOI.ConstraintName, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Get name from inner model @@ -378,10 +378,10 @@ function MOI.get( end function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintName, + m::Optimizer{T}, ::MOI.ConstraintName, c::MOI.ConstraintIndex{MOI.SingleVariable, S}, name::String -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} # Sanity checks MOI.throw_if_not_valid(m, c) c_ = get(m.name2con, name, nothing) @@ -402,10 +402,10 @@ function MOI.set( end function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintName, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S}, + m::Optimizer{T}, ::MOI.ConstraintName, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}, name::String -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Check for dupplicate name @@ -434,26 +434,26 @@ end # ConstraintFunction # function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintFunction, + m::Optimizer{T}, ::MOI.ConstraintFunction, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Sanity check return MOI.SingleVariable(MOI.VariableIndex(c.value)) end function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintFunction, + m::Optimizer{T}, ::MOI.ConstraintFunction, c::MOI.ConstraintIndex{MOI.SingleVariable, S}, ::MOI.SingleVariable -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} return throw(MOI.SettingSingleVariableFunctionNotAllowed()) end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintFunction, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, ::MOI.ConstraintFunction, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Sanity check # Get row from inner model @@ -462,22 +462,22 @@ function MOI.get( nz = length(row.nzind) # Map inner indices to MOI indices - terms = Vector{MOI.ScalarAffineTerm{Tv}}(undef, nz) + terms = Vector{MOI.ScalarAffineTerm{T}}(undef, nz) for (k, (j, v)) in enumerate(zip(row.nzind, row.nzval)) - terms[k] = MOI.ScalarAffineTerm{Tv}(v, m.var_indices_moi[j]) + terms[k] = MOI.ScalarAffineTerm{T}(v, m.var_indices_moi[j]) end - return MOI.ScalarAffineFunction(terms, zero(Tv)) + return MOI.ScalarAffineFunction(terms, zero(T)) end # TODO function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintFunction, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S}, - f::MOI.ScalarAffineFunction{Tv} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, ::MOI.ConstraintFunction, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}, + f::MOI.ScalarAffineFunction{T} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) - iszero(f.constant) || throw(MOI.ScalarFunctionConstantNotZero{Tv, typeof(f), S}(f.constant)) + iszero(f.constant) || throw(MOI.ScalarFunctionConstantNotZero{T, typeof(f), S}(f.constant)) fc = MOI.Utilities.canonical(f) @@ -488,7 +488,7 @@ function MOI.set( f_old = MOI.get(m, MOI.ConstraintFunction(), c) for term in f_old.terms j = m.var_indices[term.variable_index] - set_coefficient!(m.inner.pbdata, i, j, zero(Tv)) + set_coefficient!(m.inner.pbdata, i, j, zero(T)) end # Set new row coefficients for term in fc.terms @@ -505,9 +505,9 @@ end # ConstraintSet # function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Tv}} -) where{Tv} + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{T}} +) where{T} # Sanity check MOI.throw_if_not_valid(m, c) v = MOI.VariableIndex(c.value) @@ -520,9 +520,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Tv}} -) where{Tv} + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{T}} +) where{T} # Sanity check MOI.throw_if_not_valid(m, c) v = MOI.VariableIndex(c.value) @@ -535,9 +535,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Tv}} -) where{Tv} + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{T}} +) where{T} # Sanity check MOI.throw_if_not_valid(m, c) v = MOI.VariableIndex(c.value) @@ -550,9 +550,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Tv}} -) where{Tv} + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{T}} +) where{T} # Sanity check MOI.throw_if_not_valid(m, c) v = MOI.VariableIndex(c.value) @@ -566,9 +566,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Sanity check # Get inner bounds @@ -576,22 +576,22 @@ function MOI.get( lb = m.inner.pbdata.lcon[i] ub = m.inner.pbdata.ucon[i] - if S == MOI.LessThan{Tv} + if S == MOI.LessThan{T} return MOI.LessThan(ub) - elseif S == MOI.GreaterThan{Tv} + elseif S == MOI.GreaterThan{T} return MOI.GreaterThan(lb) - elseif S == MOI.EqualTo{Tv} + elseif S == MOI.EqualTo{T} return MOI.EqualTo(lb) - elseif S == MOI.Interval{Tv} + elseif S == MOI.Interval{T} return MOI.Interval(lb, ub) end end function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintSet, + m::Optimizer{T}, ::MOI.ConstraintSet, c::MOI.ConstraintIndex{MOI.SingleVariable, S}, s::S -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} # Sanity check MOI.throw_if_not_valid(m, c) v = MOI.VariableIndex(c.value) @@ -599,14 +599,14 @@ function MOI.set( # Update inner bounds # Bound key does not need to be updated j = m.var_indices[v] - if S == MOI.LessThan{Tv} + if S == MOI.LessThan{T} set_attribute(m.inner, VariableUpperBound(), j, s.upper) - elseif S == MOI.GreaterThan{Tv} + elseif S == MOI.GreaterThan{T} set_attribute(m.inner, VariableLowerBound(), j, s.lower) - elseif S == MOI.EqualTo{Tv} + elseif S == MOI.EqualTo{T} set_attribute(m.inner, VariableLowerBound(), j, s.value) set_attribute(m.inner, VariableUpperBound(), j, s.value) - elseif S == MOI.Interval{Tv} + elseif S == MOI.Interval{T} set_attribute(m.inner, VariableLowerBound(), j, s.lower) set_attribute(m.inner, VariableUpperBound(), j, s.upper) else @@ -617,22 +617,22 @@ function MOI.set( end function MOI.set( - m::Optimizer{Tv}, ::MOI.ConstraintSet, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S}, + m::Optimizer{T}, ::MOI.ConstraintSet, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S}, s::S -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) # Update inner bounds i = m.con_indices[c] - if S == MOI.LessThan{Tv} + if S == MOI.LessThan{T} set_attribute(m.inner, ConstraintUpperBound(), i, s.upper) - elseif S == MOI.GreaterThan{Tv} + elseif S == MOI.GreaterThan{T} set_attribute(m.inner, ConstraintLowerBound(), i, s.lower) - elseif S == MOI.EqualTo{Tv} + elseif S == MOI.EqualTo{T} set_attribute(m.inner, ConstraintLowerBound(), i, s.value) set_attribute(m.inner, ConstraintUpperBound(), i, s.value) - elseif S == MOI.Interval{Tv} + elseif S == MOI.Interval{T} set_attribute(m.inner, ConstraintLowerBound(), i, s.lower) set_attribute(m.inner, ConstraintUpperBound(), i, s.upper) else @@ -646,9 +646,9 @@ end # ConstraintPrimal # function MOI.get( - m::Optimizer{Tv}, attr::MOI.ConstraintPrimal, + m::Optimizer{T}, attr::MOI.ConstraintPrimal, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) MOI.check_result_index_bounds(m, attr) @@ -658,9 +658,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, attr::MOI.ConstraintPrimal, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, attr::MOI.ConstraintPrimal, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) MOI.check_result_index_bounds(m, attr) @@ -673,9 +673,9 @@ end # ConstraintDual # function MOI.get( - m::Optimizer{Tv}, attr::MOI.ConstraintDual, + m::Optimizer{T}, attr::MOI.ConstraintDual, c::MOI.ConstraintIndex{MOI.SingleVariable, S} -) where{Tv, S<:SCALAR_SETS{Tv}} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) MOI.check_result_index_bounds(m, attr) @@ -683,9 +683,9 @@ function MOI.get( j = m.var_indices[MOI.VariableIndex(c.value)] # Extract reduced cost - if S == MOI.LessThan{Tv} + if S == MOI.LessThan{T} return -m.inner.solution.s_upper[j] - elseif S == MOI.GreaterThan{Tv} + elseif S == MOI.GreaterThan{T} return m.inner.solution.s_lower[j] else return m.inner.solution.s_lower[j] - m.inner.solution.s_upper[j] @@ -693,9 +693,9 @@ function MOI.get( end function MOI.get( - m::Optimizer{Tv}, attr::MOI.ConstraintDual, - c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Tv}, S} -) where{Tv, S<:SCALAR_SETS{Tv}} + m::Optimizer{T}, attr::MOI.ConstraintDual, + c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, S} +) where{T, S<:SCALAR_SETS{T}} MOI.throw_if_not_valid(m, c) MOI.check_result_index_bounds(m, attr) diff --git a/src/Interfaces/MOI/objective.jl b/src/Interfaces/MOI/objective.jl index 2ea47b2b..9e835908 100644 --- a/src/Interfaces/MOI/objective.jl +++ b/src/Interfaces/MOI/objective.jl @@ -1,17 +1,17 @@ # ============================================= # 1. Supported objectives # ============================================= -MOI.supports(::Optimizer{Tv}, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Tv}}) where{Tv} = true +MOI.supports(::Optimizer{T}, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}) where{T} = true # ============================================= # 2. Get/set objective function # ============================================= function MOI.get( - m::Optimizer{Tv}, - ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Tv}} -) where{Tv} + m::Optimizer{T}, + ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}} +) where{T} # Objective coeffs - terms = MOI.ScalarAffineTerm{Tv}[] + terms = MOI.ScalarAffineTerm{T}[] for (j, cj) in enumerate(m.inner.pbdata.obj) !iszero(cj) && push!(terms, MOI.ScalarAffineTerm(cj, m.var_indices_moi[j])) end @@ -24,10 +24,10 @@ end # TODO: use inner API function MOI.set( - m::Optimizer{Tv}, - ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Tv}}, - f::MOI.ScalarAffineFunction{Tv} -) where{Tv} + m::Optimizer{T}, + ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}, + f::MOI.ScalarAffineFunction{T} +) where{T} # Sanity checks isfinite(f.constant) || error("Objective constant term must be finite") @@ -36,7 +36,7 @@ function MOI.set( end # Update inner model - m.inner.pbdata.obj .= zero(Tv) # Reset inner objective to zero + m.inner.pbdata.obj .= zero(T) # Reset inner objective to zero for t in f.terms j = m.var_indices[t.variable_index] m.inner.pbdata.obj[j] += t.coefficient # there may be dupplicates @@ -50,10 +50,10 @@ end # 3. Modify objective # ============================================= function MOI.modify( - m::Optimizer{Tv}, - c::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Tv}}, - chg::MOI.ScalarCoefficientChange{Tv} -) where{Tv} + m::Optimizer{T}, + c::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}, + chg::MOI.ScalarCoefficientChange{T} +) where{T} # Sanity checks v = chg.variable MOI.throw_if_not_valid(m, v) @@ -65,10 +65,10 @@ function MOI.modify( end function MOI.modify( - m::Optimizer{Tv}, - ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Tv}}, - chg::MOI.ScalarConstantChange{Tv} -) where{Tv} + m::Optimizer{T}, + ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}, + chg::MOI.ScalarConstantChange{T} +) where{T} isfinite(chg.new_constant) || error("Objective constant term must be finite") m.inner.pbdata.obj0 = chg.new_constant return nothing diff --git a/src/Interfaces/MOI/variables.jl b/src/Interfaces/MOI/variables.jl index 9f55cce5..8876bc62 100644 --- a/src/Interfaces/MOI/variables.jl +++ b/src/Interfaces/MOI/variables.jl @@ -24,11 +24,11 @@ function MOI.is_valid(m::Optimizer, x::MOI.VariableIndex) return haskey(m.var_indices, x) end -function MOI.add_variable(m::Optimizer{Tv}) where{Tv} +function MOI.add_variable(m::Optimizer{T}) where{T} # TODO: dispatch a function call to m.inner instead of m.inner.pbdata m.var_counter += 1 x = MOI.VariableIndex(m.var_counter) - j = Tulip.add_variable!(m.inner.pbdata, Int[], Tv[], zero(Tv), Tv(-Inf), Tv(Inf)) + j = Tulip.add_variable!(m.inner.pbdata, Int[], T[], zero(T), T(-Inf), T(Inf)) # Update tracking of variables m.var_indices[x] = j @@ -110,10 +110,10 @@ function MOI.set(m::Optimizer, ::MOI.VariableName, v::MOI.VariableIndex, name::S return nothing end -function MOI.get(m::Optimizer{Tv}, +function MOI.get(m::Optimizer{T}, attr::MOI.VariablePrimal, x::MOI.VariableIndex -) where{Tv} +) where{T} MOI.throw_if_not_valid(m, x) MOI.check_result_index_bounds(m, attr) diff --git a/src/Interfaces/tulip_julia_api.jl b/src/Interfaces/tulip_julia_api.jl index 1fa65d40..9956636c 100644 --- a/src/Interfaces/tulip_julia_api.jl +++ b/src/Interfaces/tulip_julia_api.jl @@ -8,27 +8,27 @@ using QPSReader """ - load_problem!(m::Model{Tv}, fname::String) + load_problem!(m::Model{T}, fname::String) Read a model from file `fname` and load it into model `m`. Only free MPS files are currently supported. """ -function load_problem!(m::Model{Tv}, fname::String) where{Tv} +function load_problem!(m::Model{T}, fname::String) where{T} Base.empty!(m) dat = with_logger(Logging.NullLogger()) do readqps(fname, mpsformat=:free) end - # TODO: avoid allocations when Tv is Float64 + # TODO: avoid allocations when T is Float64 objsense = !(dat.objsense == :max) load_problem!(m.pbdata, dat.name, - objsense, Tv.(dat.c), Tv(dat.c0), - sparse(dat.arows, dat.acols, Tv.(dat.avals), dat.ncon, dat.nvar), - Tv.(dat.lcon), Tv.(dat.ucon), - Tv.(dat.lvar), Tv.(dat.uvar), + objsense, T.(dat.c), T(dat.c0), + sparse(dat.arows, dat.acols, T.(dat.avals), dat.ncon, dat.nvar), + T.(dat.lcon), T.(dat.ucon), + T.(dat.lvar), T.(dat.uvar), dat.connames, dat.varnames ) @@ -72,11 +72,11 @@ function get_attribute(m::Model, ::BarrierIterations) end """ - set_attribute(m::Model{Tv}, ::VariableLowerBound, j::Int, lb::Tv) + set_attribute(m::Model{T}, ::VariableLowerBound, j::Int, lb::T) Set the lower bound of variable `j` in model `m` to `lb`. """ -function set_attribute(m::Model{Tv}, ::VariableLowerBound, j::Int, lb::Tv) where{Tv} +function set_attribute(m::Model{T}, ::VariableLowerBound, j::Int, lb::T) where{T} # sanity checks 1 <= j <= m.pbdata.nvar || error("Invalid variable index $j") @@ -86,7 +86,7 @@ function set_attribute(m::Model{Tv}, ::VariableLowerBound, j::Int, lb::Tv) where end """ - get_attribute(m::Model{Tv}, ::VariableLowerBound, j::Int) + get_attribute(m::Model{T}, ::VariableLowerBound, j::Int) Query the lower bound of variable `j` in model `m`. """ @@ -99,11 +99,11 @@ function get_attribute(m::Model, ::VariableLowerBound, j::Int) end """ - set_attribute(m::Model{Tv}, ::VariableUpperBound, j::Int, ub::Tv) + set_attribute(m::Model{T}, ::VariableUpperBound, j::Int, ub::T) Set the upper bound of variable `j` in model `m` to `ub`. """ -function set_attribute(m::Model{Tv}, ::VariableUpperBound, j::Int, ub::Tv) where{Tv} +function set_attribute(m::Model{T}, ::VariableUpperBound, j::Int, ub::T) where{T} # sanity checks 1 <= j <= m.pbdata.nvar || error("Invalid variable index $j") @@ -113,11 +113,11 @@ function set_attribute(m::Model{Tv}, ::VariableUpperBound, j::Int, ub::Tv) where end """ - set_attribute(m::Model{Tv}, ::ConstraintLowerBound, i::Int, lb::Tv) + set_attribute(m::Model{T}, ::ConstraintLowerBound, i::Int, lb::T) Set the lower bound of constraint `i` in model `m` to `lb`. """ -function set_attribute(m::Model{Tv}, ::ConstraintLowerBound, i::Int, lb::Tv) where{Tv} +function set_attribute(m::Model{T}, ::ConstraintLowerBound, i::Int, lb::T) where{T} # sanity checks 1 <= i <= m.pbdata.nvar || error("Invalid constraint index $i") @@ -127,11 +127,11 @@ function set_attribute(m::Model{Tv}, ::ConstraintLowerBound, i::Int, lb::Tv) whe end """ - set_attribute(m::Model{Tv}, ::ConstraintUpperBound, i::Int, ub::Tv) + set_attribute(m::Model{T}, ::ConstraintUpperBound, i::Int, ub::T) Set the upper bound of constraint `i` in model `m` to `ub`. """ -function set_attribute(m::Model{Tv}, ::ConstraintUpperBound, i::Int, ub::Tv) where{Tv} +function set_attribute(m::Model{T}, ::ConstraintUpperBound, i::Int, ub::T) where{T} # sanity checks 1 <= i <= m.pbdata.nvar || error("Invalid constraint index $i") @@ -206,14 +206,14 @@ end # TODO: Query solution value get_attribute(m::Model, ::ObjectiveConstant) = m.pbdata.obj0 -set_attribute(m::Model{Tv}, ::ObjectiveConstant, obj0::Tv) where{Tv} = (m.pbdata.obj0 = obj0; return nothing) +set_attribute(m::Model{T}, ::ObjectiveConstant, obj0::T) where{T} = (m.pbdata.obj0 = obj0; return nothing) """ get_attribute(model::Model, ::ObjectiveValue) Query the `ObjectiveValue` attribute from `model` """ -function get_attribute(m::Model{Tv}, ::ObjectiveValue) where{Tv} +function get_attribute(m::Model{T}, ::ObjectiveValue) where{T} if isnothing(m.solution) error("Model has no solution") else @@ -227,7 +227,7 @@ end Query the `DualObjectiveValue` attribute from `model` """ -function get_attribute(m::Model{Tv}, ::DualObjectiveValue) where{Tv} +function get_attribute(m::Model{T}, ::DualObjectiveValue) where{T} if isnothing(m.solution) error("Model has no solution") else diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index c5bdef54..f193cc2a 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -5,7 +5,7 @@ using LinearAlgebra export AbstractKKTSolver """ - AbstractKKTSolver{Tv} + AbstractKKTSolver{T} Abstract container for solving an augmented system ``` @@ -14,7 +14,7 @@ Abstract container for solving an augmented system ``` where `ξd` and `ξp` are given right-hand side. """ -abstract type AbstractKKTSolver{Tv} end +abstract type AbstractKKTSolver{T} end """ SolverOptions @@ -54,10 +54,10 @@ After this call, `kkt` can be used to solve the augmented system for given right-hand sides `ξd` and `ξp`. # Arguments -* `kkt::AbstractKKTSolver{Tv}`: the KKT solver object -* `θinv::AbstractVector{Tv}`: ``θ⁻¹`` -* `regP::AbstractVector{Tv}`: primal regularizations -* `regD::AbstractVector{Tv}`: dual regularizations +* `kkt::AbstractKKTSolver{T}`: the KKT solver object +* `θinv::AbstractVector{T}`: ``θ⁻¹`` +* `regP::AbstractVector{T}`: primal regularizations +* `regD::AbstractVector{T}`: dual regularizations """ function update! end @@ -84,7 +84,7 @@ function solve! end Return the arithmetic used by the solver. """ -arithmetic(kkt::AbstractKKTSolver{Tv}) where Tv = Tv +arithmetic(kkt::AbstractKKTSolver{T}) where T = T """ backend(kkt) @@ -110,15 +110,15 @@ include("ldlfact.jl") include("krylov.jl") """ - default_options(Tv) + default_options(T) Use CHOLMOD for `Float64` and LDLFactorizations otherwise. """ -function default_options(::Type{Tv}) where{Tv} - if Tv == Float64 +function default_options(::Type{T}) where{T} + if T == Float64 return SolverOptions(CholmodSolver; normal_equations=false) else - return SolverOptions(LDLFact_SymQuasDef) + return SolverOptions(LDLFactSQD) end end diff --git a/src/KKT/cholmod.jl b/src/KKT/cholmod.jl index 5a0412ce..12b10492 100644 --- a/src/KKT/cholmod.jl +++ b/src/KKT/cholmod.jl @@ -7,13 +7,13 @@ using SuiteSparse.CHOLMOD Uses CHOLMOD's factorization routines to solve the augmented system. To use an LDLᵀ factorization of the augmented system -(see [`Cholmod_SymQuasDef`](@ref)) +(see [`CholmodSQD`](@ref)) ```julia model.params.KKTOptions = Tulip.KKT.SolverOptions(CholmodSolver, normal_equations=false) ``` To use a Cholesky factorization of the normal equations system -(see [`Cholmod_SymPosDef`](@ref)) +(see [`CholmodSPD`](@ref)) ```julia model.params.KKTOptions = Tulip.KKT.SolverOptions(CholmodSolver, normal_equations=true) ``` @@ -29,25 +29,25 @@ function CholmodSolver( kwargs... ) if normal_equations - return Cholmod_SymPosDef(A; kwargs...) + return CholmodSPD(A; kwargs...) else - return Cholmod_SymQuasDef(A; kwargs...) + return CholmodSQD(A; kwargs...) end end # ============================================================================== -# Cholmod_SymQuasDef +# CholmodSQD # ============================================================================== """ - Cholmod_SymQuasDef + CholmodSQD Linear solver for the 2x2 augmented system with ``A`` sparse. Uses an LDLᵀ factorization of the quasi-definite augmented system. """ -mutable struct Cholmod_SymQuasDef <: CholmodSolver +mutable struct CholmodSQD <: CholmodSolver m::Int # Number of rows n::Int # Number of columns @@ -67,7 +67,7 @@ mutable struct Cholmod_SymQuasDef <: CholmodSolver F::CHOLMOD.Factor{Float64} # TODO: constructor with initial memory allocation - function Cholmod_SymQuasDef(A::AbstractMatrix{Float64}) + function CholmodSQD(A::AbstractMatrix{Float64}) m, n = size(A) θ = ones(Float64, n) @@ -84,8 +84,8 @@ mutable struct Cholmod_SymQuasDef <: CholmodSolver end -backend(::Cholmod_SymQuasDef) = "CHOLMOD" -linear_system(::Cholmod_SymQuasDef) = "Augmented system" +backend(::CholmodSQD) = "CHOLMOD" +linear_system(::CholmodSQD) = "Augmented system" """ update!(kkt, θ, regP, regD) @@ -97,7 +97,7 @@ Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute Throws a `PosDefException` if matrix is not quasi-definite. """ function update!( - kkt::Cholmod_SymQuasDef, + kkt::CholmodSQD, θ::AbstractVector{Float64}, regP::AbstractVector{Float64}, regD::AbstractVector{Float64} @@ -143,7 +143,7 @@ Solve the augmented system, overwriting `dx, dy` with the result. """ function solve!( dx::Vector{Float64}, dy::Vector{Float64}, - kkt::Cholmod_SymQuasDef, + kkt::CholmodSQD, ξp::Vector{Float64}, ξd::Vector{Float64} ) m, n = kkt.m, kkt.n @@ -178,11 +178,11 @@ function solve!( end # ============================================================================== -# Cholmod_SymPosDef +# CholmodSPD # ============================================================================== """ - Cholmod_SymPosDef + CholmodSPD Linear solver for the 2x2 augmented system ``` @@ -197,7 +197,7 @@ Uses a Cholesky factorization of the positive definite normal equations system dx = (Θ⁻¹ + Rp)⁻¹ * (Aᵀ * dy - xi_d) ``` """ -mutable struct Cholmod_SymPosDef <: CholmodSolver +mutable struct CholmodSPD <: CholmodSolver m::Int # Number of rows n::Int # Number of columns @@ -216,7 +216,7 @@ mutable struct Cholmod_SymPosDef <: CholmodSolver # Constructor and initial memory allocation # TODO: symbolic only + allocation - function Cholmod_SymPosDef(A::AbstractMatrix{Float64}) + function CholmodSPD(A::AbstractMatrix{Float64}) m, n = size(A) θ = ones(Float64, n) @@ -230,8 +230,8 @@ mutable struct Cholmod_SymPosDef <: CholmodSolver end -backend(::Cholmod_SymPosDef) = "CHOLMOD - Cholesky" -linear_system(::Cholmod_SymPosDef) = "Normal equations" +backend(::CholmodSPD) = "CHOLMOD - Cholesky" +linear_system(::CholmodSPD) = "Normal equations" """ update!(kkt, θ, regP, regD) @@ -239,7 +239,7 @@ linear_system(::Cholmod_SymPosDef) = "Normal equations" Compute normal equation system matrix, and update the factorization. """ function update!( - kkt::Cholmod_SymPosDef, + kkt::CholmodSPD, θ::AbstractVector{Float64}, regP::AbstractVector{Float64}, regD::AbstractVector{Float64} @@ -280,7 +280,7 @@ Solve the augmented system, overwriting `dx, dy` with the result. """ function solve!( dx::Vector{Float64}, dy::Vector{Float64}, - kkt::Cholmod_SymPosDef, + kkt::CholmodSPD, ξp::Vector{Float64}, ξd::Vector{Float64} ) m, n = kkt.m, kkt.n diff --git a/src/KKT/krylov.jl b/src/KKT/krylov.jl index 952476d2..c6821839 100644 --- a/src/KKT/krylov.jl +++ b/src/KKT/krylov.jl @@ -9,11 +9,11 @@ Abstract type for Kyrlov-based linear solvers. abstract type AbstractKrylovSolver{T} <: AbstractKKTSolver{T} end # ============================================================================== -# KrylovSPDSolver: +# KrylovSPD: # ============================================================================== """ - KrylovSPDSolver{T, F, Tv, Ta} + KrylovSPD{T, F, Tv, Ta} Krylov-based **S**ymmetric **P**ositive-**D**efinite (SPD) linear solver. @@ -22,10 +22,10 @@ to the augmented system. The selected Krylov method must therefore handle symmetric positive-definite systems. Suitable methods are CG or CR. -A `KrylovSPDSolver` is selected as follows +A `KrylovSPD` is selected as follows ```julia model.params.KKTOptions = Tulip.KKT.SolverOptions( - KrylovSPDSolver, method=Krylov.cg, + KrylovSPD, method=Krylov.cg, atol=1e-12, rtol=1e-12 ) ``` @@ -38,7 +38,7 @@ where `S`, `ξ` and `dy` are the normal equations system's left- and right-hand side, and solution vector, respectively. `S` may take the form of a matrix, or of a suitable linear operator. """ -mutable struct KrylovSPDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} +mutable struct KrylovSPD{T, F, Tv, Ta} <: AbstractKrylovSolver{T} m::Int n::Int @@ -52,7 +52,7 @@ mutable struct KrylovSPDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} regP::Tv # primal regularization regD::Tv # dual regularization - function KrylovSPDSolver(f::Function, A::Ta; + function KrylovSPD(f::Function, A::Ta; atol::T=eps(T)^(3 // 4), rtol::T=eps(T)^(3 // 4) ) where{T, Ta <: AbstractMatrix{T}} @@ -67,23 +67,23 @@ mutable struct KrylovSPDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} end end -setup(::Type{KrylovSPDSolver}, A; method, kwargs...) = KrylovSPDSolver(method, A; kwargs...) +setup(::Type{KrylovSPD}, A; method, kwargs...) = KrylovSPD(method, A; kwargs...) -backend(kkt::KrylovSPDSolver) = "Krylov ($(kkt.f))" -linear_system(::KrylovSPDSolver) = "Normal equations" +backend(kkt::KrylovSPD) = "Krylov ($(kkt.f))" +linear_system(::KrylovSPD) = "Normal equations" """ - update!(kkt::KrylovSPDSolver, θ, regP, regD) + update!(kkt::KrylovSPD, θ, regP, regD) Update diagonal scaling ``θ`` and primal-dual regularizations. """ function update!( - kkt::KrylovSPDSolver{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:Real} + kkt::KrylovSPD{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -105,19 +105,19 @@ function update!( end """ - solve!(dx, dy, kkt::KrylovSPDSolver, ξp, ξd) + solve!(dx, dy, kkt::KrylovSPD, ξp, ξd) Solve the normal equations system using the selected Krylov method. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::KrylovSPDSolver{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:Real} + dx::Vector{T}, dy::Vector{T}, + kkt::KrylovSPD{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T} m, n = kkt.m, kkt.n # Setup - d = one(Tv) ./ (kkt.θ .+ kkt.regP) + d = one(T) ./ (kkt.θ .+ kkt.regP) D = Diagonal(d) # Set-up right-hand side @@ -125,9 +125,9 @@ function solve!( # Form linear operator S = A * D * A' + Rd # Here we pre-allocate the intermediary and final vectors - v1 = zeros(Tv, n) - v2 = zeros(Tv, m) - opS = LO.LinearOperator(Tv, m, m, true, true, + v1 = zeros(T, n) + v2 = zeros(T, m) + opS = LO.LinearOperator(T, m, m, true, true, w -> begin mul!(v1, kkt.A', w) v1 .*= d @@ -149,11 +149,11 @@ end # ============================================================================== -# KrylovSIDSolver: +# KrylovSID: # ============================================================================== """ - KrylovSIDSolver{T, F, Tv, Ta} + KrylovSID{T, F, Tv, Ta} Krylov-based **S**ymmetric **I**n**D**efinite (SID) linear solver. @@ -162,10 +162,10 @@ without exploiting its quasi-definiteness properties. The selected Krylov method must therefore handle symmetric indefinite systems. Suitable methods are MINRES or MINRES-QLP. -A `KrylovSIDSolver` is selected as follows +A `KrylovSID` is selected as follows ```julia model.params.KKTOptions = Tulip.KKT.SolverOptions( - KrylovSIDSolver, method=Krylov.minres, + KrylovSID, method=Krylov.minres, atol=1e-12, rtol=1e-12 ) ``` @@ -178,7 +178,7 @@ where `K`, `ξ` and `Δ` are the augmented system's left- and right-hand side, and solution vector, respectively. `K` may take the form of a matrix, or of a suitable linear operator. """ -mutable struct KrylovSIDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} +mutable struct KrylovSID{T, F, Tv, Ta} <: AbstractKrylovSolver{T} m::Int n::Int @@ -192,7 +192,7 @@ mutable struct KrylovSIDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} regP::Tv # primal regularization regD::Tv # dual regularization - function KrylovSIDSolver(f::Function, A::Ta; + function KrylovSID(f::Function, A::Ta; atol::T=eps(T)^(3 // 4), rtol::T=eps(T)^(3 // 4) ) where{T, Ta <: AbstractMatrix{T}} @@ -207,23 +207,23 @@ mutable struct KrylovSIDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} end end -setup(::Type{KrylovSIDSolver}, A; method, kwargs...) = KrylovSIDSolver(method, A; kwargs...) +setup(::Type{KrylovSID}, A; method, kwargs...) = KrylovSID(method, A; kwargs...) -backend(kkt::KrylovSIDSolver) = "Krylov ($(kkt.f))" -linear_system(::KrylovSIDSolver) = "Augmented system" +backend(kkt::KrylovSID) = "Krylov ($(kkt.f))" +linear_system(::KrylovSID) = "Augmented system" """ - update!(kkt::KrylovSIDSolver, θ, regP, regD) + update!(kkt::KrylovSID, θ, regP, regD) Update diagonal scaling ``θ`` and primal-dual regularizations. """ function update!( - kkt::KrylovSIDSolver{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:Real} + kkt::KrylovSID{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -245,15 +245,15 @@ function update!( end """ - solve!(dx, dy, kkt::KrylovSIDSolver, ξp, ξd) + solve!(dx, dy, kkt::KrylovSID, ξp, ξd) Solve the augmented system using the selected Krylov method. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::KrylovSIDSolver{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:Real} + dx::Vector{T}, dy::Vector{T}, + kkt::KrylovSID{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T} m, n = kkt.m, kkt.n # Setup @@ -265,8 +265,8 @@ function solve!( # Form linear operator K = [-(Θ⁻¹ + Rp) Aᵀ; A Rd] # Here we pre-allocate the final vector - z = zeros(Tv, m+n) - opK = LO.LinearOperator(Tv, n+m, n+m, true, false, + z = zeros(T, m+n) + opK = LO.LinearOperator(T, n+m, n+m, true, false, w -> begin @views z1 = z[1:n] @views z2 = z[n+1:n+m] @@ -275,12 +275,12 @@ function solve!( @views w2 = w[n+1:n+m] # z1 = -(Θ⁻¹ + Rp) w1 + A'w2 - mul!(z1, D, w1, one(Tv), zero(Tv)) # z1 = -(Θ⁻¹ + Rp) w1 - mul!(z1, A', w2, one(Tv), one(Tv)) # z1 += A'w2 + mul!(z1, D, w1, one(T), zero(T)) # z1 = -(Θ⁻¹ + Rp) w1 + mul!(z1, A', w2, one(T), one(T)) # z1 += A'w2 # z2 = A w1 + Rd w2 - mul!(z2, A, w1, one(Tv), zero(Tv)) # z2 = A w1 - mul!(z2, Diagonal(kkt.regD), w2, one(Tv), one(Tv)) # z2 += Rd * w2 + mul!(z2, A, w1, one(T), zero(T)) # z2 = A w1 + mul!(z2, Diagonal(kkt.regD), w2, one(T), one(T)) # z2 += Rd * w2 z end @@ -298,11 +298,11 @@ end # ============================================================================== -# KrylovSQDSolver: +# KrylovSQD: # ============================================================================== """ - KrylovSQDSolver{T, F, Tv, Ta} + KrylovSQD{T, F, Tv, Ta} Krylov-based **S**ymmetric **Q**uasi-**D**efinite (SQD) linear solver. @@ -311,10 +311,10 @@ its 2x2 block structure and quasi-definiteness. The selected Krylov method must therefore handle 2x2 symmetric quasi-definite systems. Suitable methods are TriCG and TriMR. -A `KrylovSIDSolver` is selected as follows +A `KrylovSID` is selected as follows ```julia model.params.KKTOptions = Tulip.KKT.SolverOptions( - KrylovSQDSolver, method=Krylov.trimr, + KrylovSQD, method=Krylov.trimr, atol=1e-12, rtol=1e-12 ) ``` @@ -330,7 +330,7 @@ where the augmented system is of the form ``` i.e., ``N = (Θ^{-1} + R_{p})^{-1}`` and ``M = R_{d}^{-1}``. """ -mutable struct KrylovSQDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} +mutable struct KrylovSQD{T, F, Tv, Ta} <: AbstractKrylovSolver{T} m::Int n::Int @@ -344,7 +344,7 @@ mutable struct KrylovSQDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} regP::Tv # primal regularization regD::Tv # dual regularization - function KrylovSQDSolver(f::Function, A::Ta; + function KrylovSQD(f::Function, A::Ta; atol::T=eps(T)^(3 // 4), rtol::T=eps(T)^(3 // 4) ) where{T, Ta <: AbstractMatrix{T}} @@ -359,23 +359,23 @@ mutable struct KrylovSQDSolver{T, F, Tv, Ta} <: AbstractKrylovSolver{T} end end -setup(::Type{KrylovSQDSolver}, A; method, kwargs...) = KrylovSQDSolver(method, A; kwargs...) +setup(::Type{KrylovSQD}, A; method, kwargs...) = KrylovSQD(method, A; kwargs...) -backend(kkt::KrylovSQDSolver) = "Krylov ($(kkt.f))" -linear_system(::KrylovSQDSolver) = "Augmented system" +backend(kkt::KrylovSQD) = "Krylov ($(kkt.f))" +linear_system(::KrylovSQD) = "Augmented system" """ - update!(kkt::KrylovSQDSolver, θ, regP, regD) + update!(kkt::KrylovSQD, θ, regP, regD) Update diagonal scaling ``θ`` and primal-dual regularizations. """ function update!( - kkt::KrylovSQDSolver{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:Real} + kkt::KrylovSQD{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -397,22 +397,22 @@ function update!( end """ - solve!(dx, dy, kkt::KrylovSQDSolver, ξp, ξd) + solve!(dx, dy, kkt::KrylovSQD, ξp, ξd) Solve the augmented system using the selected Kyrlov method. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::KrylovSQDSolver{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:Real} + dx::Vector{T}, dy::Vector{T}, + kkt::KrylovSQD{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T} m, n = kkt.m, kkt.n # Setup - d = one(Tv) ./ (kkt.θ .+ kkt.regP) + d = one(T) ./ (kkt.θ .+ kkt.regP) N = Diagonal(d) - M = Diagonal(one(Tv) ./ kkt.regD) + M = Diagonal(one(T) ./ kkt.regD) # Solve augmented system _dy, _dx, stats = kkt.f(kkt.A, ξp, ξd, diff --git a/src/KKT/lapack.jl b/src/KKT/lapack.jl index 3b9e8639..06140dff 100644 --- a/src/KKT/lapack.jl +++ b/src/KKT/lapack.jl @@ -3,7 +3,7 @@ using LinearAlgebra.LAPACK const BlasReal = LinearAlgebra.BlasReal """ - Dense_SymPosDef{Tv} + DenseSPD{T} KKT solver for dense linear systems. @@ -12,58 +12,58 @@ Uses a Cholesky factorization of the normal equations system. Not recommended ```julia model = Tulip.Model{Float64}() -model.params.KKTOptions = Tulip.KKT.SolverOptions(Tulip.KKT.Dense_SymPosDef) +model.params.KKTOptions = Tulip.KKT.SolverOptions(Tulip.KKT.DenseSPD) ``` !!! info Dispatches to BLAS/LAPACK in `Float32` and `Float64` arithmetic. """ -mutable struct Dense_SymPosDef{Tv<:Real} <: AbstractKKTSolver{Tv} +mutable struct DenseSPD{T} <: AbstractKKTSolver{T} m::Int # Number of rows in A n::Int # Number of columns in A # Problem data - A::Matrix{Tv} - θ::Vector{Tv} + A::Matrix{T} + θ::Vector{T} # Regularizations - regP::Vector{Tv} # primal - regD::Vector{Tv} # dual + regP::Vector{T} # primal + regD::Vector{T} # dual - _A::Matrix{Tv} # place-holder to avoid allocating memory afterwards + _A::Matrix{T} # place-holder to avoid allocating memory afterwards # Factorization - S::Matrix{Tv} # Normal equations matrix, that also holds the factorization + S::Matrix{T} # Normal equations matrix, that also holds the factorization # Constructor - function Dense_SymPosDef(A::AbstractMatrix{Tv}) where{Tv<:Real} + function DenseSPD(A::AbstractMatrix{T}) where{T} m, n = size(A) - θ = ones(Tv, n) + θ = ones(T, n) # We just need to allocate memory here, # so no factorization yet - S = zeros(Tv, m, m) + S = zeros(T, m, m) - return new{Tv}(m, n, A, θ, zeros(Tv, n), zeros(Tv, m), Matrix(A), S) + return new{T}(m, n, A, θ, zeros(T, n), zeros(T, m), Matrix(A), S) end end -backend(::Dense_SymPosDef) = "LAPACK" -linear_system(::Dense_SymPosDef) = "Normal equations" +backend(::DenseSPD) = "LAPACK" +linear_system(::DenseSPD) = "Normal equations" """ - update!(kkt::Dense_SymPosDef{<:Real}, θ, regP, regD) + update!(kkt::DenseSPD{T}, θ, regP, regD) Compute normal equations system matrix and update Cholesky factorization. Uses Julia's generic linear algebra. """ function update!( - kkt::Dense_SymPosDef{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:Real} + kkt::DenseSPD{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -82,7 +82,7 @@ function update!( # Re-compute normal equations matrix # There's no function that does S = A*D*A', so we cache a copy of A copyto!(kkt._A, kkt.A) - D = Diagonal(one(Tv) ./ (kkt.θ .+ kkt.regP)) + D = Diagonal(one(T) ./ (kkt.θ .+ kkt.regP)) rmul!(kkt._A, D) # B = A * D mul!(kkt.S, kkt._A, transpose(kkt.A)) # Now S = A*D*A' # TODO: do not re-compute S if only dual regularization changes @@ -97,18 +97,18 @@ function update!( end """ - update!(kkt::Dense_SymPosDef{<:BlasReal}, θ, regP, regD) + update!(kkt::DenseSPD{<:BlasReal}, θ, regP, regD) Compute normal equations system matrix and update Cholesky factorization. Uses BLAS and LAPACK routines. """ function update!( - kkt::Dense_SymPosDef{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:BlasReal} + kkt::DenseSPD{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T<:BlasReal} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -127,9 +127,9 @@ function update!( # Re-compute normal equations matrix # There's no function that does S = A*D*A', so we cache a copy of A copyto!(kkt._A, kkt.A) - D = Diagonal(sqrt.(one(Tv) ./ (kkt.θ .+ kkt.regP))) + D = Diagonal(sqrt.(one(T) ./ (kkt.θ .+ kkt.regP))) rmul!(kkt._A, D) # B = A * √D - BLAS.syrk!('U', 'N', one(Tv), kkt._A, zero(Tv), kkt.S) + BLAS.syrk!('U', 'N', one(T), kkt._A, zero(T), kkt.S) # Regularization # TODO: only update diagonal of S if only dual regularization changes @@ -152,10 +152,10 @@ Solve the augmented system, overwriting `dx, dy` with the result. Uses two generic triangular solves for solving the normal equations system. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::Dense_SymPosDef{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:Real} + dx::Vector{T}, dy::Vector{T}, + kkt::DenseSPD{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T} m, n = kkt.m, kkt.n # Set-up right-hand side @@ -181,10 +181,10 @@ Solve the augmented system, overwriting `dx, dy` with the result. Uses one LAPACK call for solving the normal equations system. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::Dense_SymPosDef{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:BlasReal} + dx::Vector{T}, dy::Vector{T}, + kkt::DenseSPD{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T<:BlasReal} m, n = kkt.m, kkt.n # Set-up right-hand side diff --git a/src/KKT/ldlfact.jl b/src/KKT/ldlfact.jl index 353eee2c..b0ca60aa 100644 --- a/src/KKT/ldlfact.jl +++ b/src/KKT/ldlfact.jl @@ -2,21 +2,21 @@ import LDLFactorizations const LDLF = LDLFactorizations # ============================================================================== -# LDLFact_SymQuasDef +# LDLFactSQD # ============================================================================== """ - LDLFact_SymQuasDef{Tv} + LDLFactSQD{T} Linear solver for the 2x2 augmented system with ``A`` sparse. Uses LDLFactorizations.jl to compute an LDLᵀ factorization of the quasi-definite augmented system. ```julia -model.params.KKTOptions = Tulip.KKT.SolverOptions(LDLFact_SymQuasDef) +model.params.KKTOptions = Tulip.KKT.SolverOptions(LDLFactSQD) ``` """ -mutable struct LDLFact_SymQuasDef{Tv<:Real} <: AbstractKKTSolver{Tv} +mutable struct LDLFactSQD{T<:Real} <: AbstractKKTSolver{T} m::Int # Number of rows n::Int # Number of columns @@ -24,40 +24,40 @@ mutable struct LDLFact_SymQuasDef{Tv<:Real} <: AbstractKKTSolver{Tv} # and add flag (default to false) to know whether user ordering should be used # Problem data - A::SparseMatrixCSC{Tv, Int} - θ::Vector{Tv} - regP::Vector{Tv} # primal regularization - regD::Vector{Tv} # dual regularization + A::SparseMatrixCSC{T, Int} + θ::Vector{T} + regP::Vector{T} # primal regularization + regD::Vector{T} # dual regularization # Left-hand side matrix - S::SparseMatrixCSC{Tv, Int} + S::SparseMatrixCSC{T, Int} # Factorization - F::LDLF.LDLFactorization{Tv} + F::LDLF.LDLFactorization{T} # TODO: constructor with initial memory allocation - function LDLFact_SymQuasDef(A::AbstractMatrix{Tv}) where{Tv<:Real} + function LDLFactSQD(A::AbstractMatrix{T}) where{T<:Real} m, n = size(A) - θ = ones(Tv, n) + θ = ones(T, n) S = [ spdiagm(0 => -θ) A'; - spzeros(Tv, m, n) spdiagm(0 => ones(m)) + spzeros(T, m, n) spdiagm(0 => ones(m)) ] # TODO: PSD-ness checks # TODO: symbolic factorization only F = LDLF.ldl_analyze(Symmetric(S)) - return new{Tv}(m, n, A, θ, ones(Tv, n), ones(Tv, m), S, F) + return new{T}(m, n, A, θ, ones(T, n), ones(T, m), S, F) end end -setup(::Type{LDLFact_SymQuasDef}, A) = LDLFact_SymQuasDef(A) +setup(::Type{LDLFactSQD}, A) = LDLFactSQD(A) -backend(::LDLFact_SymQuasDef) = "LDLFactorizations.jl" -linear_system(::LDLFact_SymQuasDef) = "Augmented system" +backend(::LDLFactSQD) = "LDLFactorizations.jl" +linear_system(::LDLFactSQD) = "Augmented system" """ update!(kkt, θ, regP, regD) @@ -69,11 +69,11 @@ Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute Throws a `PosDefException` if matrix is not quasi-definite. """ function update!( - kkt::LDLFact_SymQuasDef{Tv}, - θ::AbstractVector{Tv}, - regP::AbstractVector{Tv}, - regD::AbstractVector{Tv} -) where{Tv<:Real} + kkt::LDLFactSQD{T}, + θ::AbstractVector{T}, + regP::AbstractVector{T}, + regD::AbstractVector{T} +) where{T<:Real} # Sanity checks length(θ) == kkt.n || throw(DimensionMismatch( "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." @@ -119,10 +119,10 @@ end Solve the augmented system, overwriting `dx, dy` with the result. """ function solve!( - dx::Vector{Tv}, dy::Vector{Tv}, - kkt::LDLFact_SymQuasDef{Tv}, - ξp::Vector{Tv}, ξd::Vector{Tv} -) where{Tv<:Real} + dx::Vector{T}, dy::Vector{T}, + kkt::LDLFactSQD{T}, + ξp::Vector{T}, ξd::Vector{T} +) where{T<:Real} m, n = kkt.m, kkt.n # Set-up right-hand side diff --git a/src/KKT/test.jl b/src/KKT/test.jl index 2737056d..56079ba3 100644 --- a/src/KKT/test.jl +++ b/src/KKT/test.jl @@ -7,15 +7,15 @@ using LinearAlgebra """ function run_ls_tests( - A::AbstractMatrix{Tv}, - kkt::AbstractKKTSolver{Tv}; - atol::Tv=sqrt(eps(Tv)) -) where{Tv<:Real} + A::AbstractMatrix{T}, + kkt::AbstractKKTSolver{T}; + atol::T=sqrt(eps(T)) +) where{T} # Check that required methods are implemented @testset "Required methods" begin Tls = typeof(kkt) - V = Vector{Tv} + V = Vector{T} @test hasmethod(update!, Tuple{Tls, V, V, V}) @test hasmethod(solve!, Tuple{V, V, Tls, V, V}) end @@ -23,16 +23,16 @@ function run_ls_tests( m, n = size(A) # Factorization/pre-conditionner update - θ = ones(Tv, n) - regP = ones(Tv, n) - regD = ones(Tv, m) + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) update!(kkt, θ, regP, regD) # Solve linear system - ξp = ones(Tv, m) - ξd = ones(Tv, n) - dx = zeros(Tv, n) - dy = zeros(Tv, m) + ξp = ones(T, m) + ξd = ones(T, n) + dx = zeros(T, n) + dy = zeros(T, m) solve!(dx, dy, kkt, ξp, ξd) # Check residuals diff --git a/src/LinearAlgebra/LinearAlgebra.jl b/src/LinearAlgebra/LinearAlgebra.jl index 11e10945..12b3f54a 100644 --- a/src/LinearAlgebra/LinearAlgebra.jl +++ b/src/LinearAlgebra/LinearAlgebra.jl @@ -33,9 +33,9 @@ function construct_matrix end function construct_matrix( ::Type{Matrix}, m::Int, n::Int, - aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{Tv} -) where{Tv<:Real} - A = zeros(Tv, m, n) + aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{T} +) where{T} + A = zeros(T, m, n) # TODO: may be more efficient to first sort indices so that # A is accessed in column-major order. for(i, j, v) in zip(aI, aJ, aV) @@ -46,7 +46,7 @@ end construct_matrix( ::Type{SparseMatrixCSC}, m::Int, n::Int, - aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{Tv} -) where{Tv<:Real} = sparse(aI, aJ, aV, m, n) + aI::Vector{Int}, aJ::Vector{Int}, aV::Vector{T} +) where{T} = sparse(aI, aJ, aV, m, n) end # module \ No newline at end of file diff --git a/src/Presolve/Presolve.jl b/src/Presolve/Presolve.jl index ae7417db..c4fadca4 100644 --- a/src/Presolve/Presolve.jl +++ b/src/Presolve/Presolve.jl @@ -1,12 +1,12 @@ """ - PresolveTransformation{Tv} + PresolveTransformation{T} Abstract type for pre-solve transformations. """ -abstract type PresolveTransformation{Tv} end +abstract type PresolveTransformation{T} end """ - PresolveData{Tv} + PresolveData{T} Stores information about an LP in the form ``` @@ -21,16 +21,16 @@ whose dual writes y⁺, y⁻, s⁺, s⁻ ⩾ 0 ``` """ -mutable struct PresolveData{Tv} +mutable struct PresolveData{T} updated::Bool status::TerminationStatus # Original problem - pb0::ProblemData{Tv} + pb0::ProblemData{T} # Reduced problem # Nothing until the reduced problem is extracted - pb_red::Union{Nothing, ProblemData{Tv}} - solution::Solution{Tv} # only used if presolve solves the problem + pb_red::Union{Nothing, ProblemData{T}} + solution::Solution{T} # only used if presolve solves the problem # Presolved data @@ -44,28 +44,28 @@ mutable struct PresolveData{Tv} # Objective objsense::Bool - obj::Vector{Tv} - obj0::Tv + obj::Vector{T} + obj0::T # Current number of constraints/variables in presolved problem nrow::Int ncol::Int # Primal bounds - lrow::Vector{Tv} - urow::Vector{Tv} - lcol::Vector{Tv} - ucol::Vector{Tv} + lrow::Vector{T} + urow::Vector{T} + lcol::Vector{T} + ucol::Vector{T} # Dual bounds - ly::Vector{Tv} - uy::Vector{Tv} - ls::Vector{Tv} - us::Vector{Tv} + ly::Vector{T} + uy::Vector{T} + ls::Vector{T} + us::Vector{T} # Scaling - row_scaling::Vector{Tv} - col_scaling::Vector{Tv} + row_scaling::Vector{T} + col_scaling::Vector{T} # TODO: objective and RHS scaling # Old <-> new index mapping @@ -80,17 +80,17 @@ mutable struct PresolveData{Tv} free_col_singletons::Vector{Int} # (implied) free column singletons # TODO: set of transformations for pre-post crush - ops::Vector{PresolveTransformation{Tv}} + ops::Vector{PresolveTransformation{T}} - function PresolveData(pb::ProblemData{Tv}) where{Tv} - ps = new{Tv}() + function PresolveData(pb::ProblemData{T}) where{T} + ps = new{T}() ps.updated = false ps.status = Trm_Unknown ps.pb0 = pb ps.pb_red = nothing - ps.solution = Solution{Tv}(pb.ncon, pb.nvar) + ps.solution = Solution{T}(pb.ncon, pb.nvar) ps.nrow = pb.ncon ps.ncol = pb.nvar @@ -128,22 +128,22 @@ mutable struct PresolveData{Tv} ps.ucol = copy(pb.uvar) # Set dual bounds - ps.ly = Vector{Tv}(undef, ps.nrow) - ps.uy = Vector{Tv}(undef, ps.nrow) - ps.ls = Vector{Tv}(undef, ps.ncol) - ps.us = Vector{Tv}(undef, ps.ncol) + ps.ly = Vector{T}(undef, ps.nrow) + ps.uy = Vector{T}(undef, ps.nrow) + ps.ls = Vector{T}(undef, ps.ncol) + ps.us = Vector{T}(undef, ps.ncol) for (i, (lc, uc)) in enumerate(zip(ps.lrow, ps.urow)) - ps.ly[i] = (uc == Tv( Inf)) ? zero(Tv) : Tv(-Inf) - ps.uy[i] = (lc == Tv(-Inf)) ? zero(Tv) : Tv( Inf) + ps.ly[i] = (uc == T( Inf)) ? zero(T) : T(-Inf) + ps.uy[i] = (lc == T(-Inf)) ? zero(T) : T( Inf) end for (j, (lv, uv)) in enumerate(zip(ps.lcol, ps.ucol)) - ps.ls[j] = (uv == Tv( Inf)) ? zero(Tv) : Tv(-Inf) - ps.us[j] = (lv == Tv(-Inf)) ? zero(Tv) : Tv( Inf) + ps.ls[j] = (uv == T( Inf)) ? zero(T) : T(-Inf) + ps.us[j] = (lv == T(-Inf)) ? zero(T) : T( Inf) end # Scalings - ps.row_scaling = ones(Tv, ps.nrow) - ps.col_scaling = ones(Tv, ps.ncol) + ps.row_scaling = ones(T, ps.nrow) + ps.col_scaling = ones(T, ps.ncol) # Index mappings ps.new_con_idx = Int[] @@ -155,16 +155,16 @@ mutable struct PresolveData{Tv} ps.row_singletons = Int[] ps.free_col_singletons = Int[] - ps.ops = PresolveTransformation{Tv}[] + ps.ops = PresolveTransformation{T}[] return ps end end # Extract pre-solved problem data, to be passed to the IPM solver -function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} +function extract_reduced_problem!(ps::PresolveData{T}) where{T} - pb = ProblemData{Tv}() + pb = ProblemData{T}() pb.ncon = sum(ps.rowflag) pb.nvar = sum(ps.colflag) @@ -185,7 +185,7 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} pb.ucon = ps.urow[ps.rowflag] # Extract new rows - pb.arows = Vector{Row{Tv}}(undef, pb.ncon) + pb.arows = Vector{Row{T}}(undef, pb.ncon) inew = 0 for (iold, row) in enumerate(ps.pb0.arows) ps.rowflag[iold] || continue @@ -193,7 +193,7 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} inew += 1 # Compute new row rind = Vector{Int}(undef, ps.nzrow[iold]) - rval = Vector{Tv}(undef, ps.nzrow[iold]) + rval = Vector{T}(undef, ps.nzrow[iold]) k = 0 for (jold, aij) in zip(row.nzind, row.nzval) @@ -207,11 +207,11 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} end # Set new row - pb.arows[inew] = Row{Tv}(rind, rval) + pb.arows[inew] = Row{T}(rind, rval) end # Extract new columns - pb.acols = Vector{Col{Tv}}(undef, pb.nvar) + pb.acols = Vector{Col{T}}(undef, pb.nvar) jnew = 0 for (jold, col) in enumerate(ps.pb0.acols) ps.colflag[jold] || continue @@ -219,7 +219,7 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} jnew += 1 # Compute new row cind = Vector{Int}(undef, ps.nzcol[jold]) - cval = Vector{Tv}(undef, ps.nzcol[jold]) + cval = Vector{T}(undef, ps.nzcol[jold]) k = 0 for (iold, aij) in zip(col.nzind, col.nzval) @@ -233,7 +233,7 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} end # Set new column - pb.acols[jnew] = Col{Tv}(cind, cval) + pb.acols[jnew] = Col{T}(cind, cval) end # Variable and constraint names @@ -242,19 +242,19 @@ function extract_reduced_problem!(ps::PresolveData{Tv}) where{Tv<:Real} pb.con_names = ps.pb0.con_names[ps.rowflag] # Scaling - rscale = zeros(Tv, ps.nrow) - cscale = zeros(Tv, ps.ncol) + rscale = zeros(T, ps.nrow) + cscale = zeros(T, ps.ncol) # Compute norm of each row and column # TODO: use a parameter p and do norm(.., p) p = 2 for (i, row) in enumerate(pb.arows) r = norm(row.nzval, p) - rscale[i] = r > zero(Tv) ? r : one(Tv) + rscale[i] = r > zero(T) ? r : one(T) end for (j, col) in enumerate(pb.acols) r = norm(col.nzval, p) - cscale[j] = r > zero(Tv) ? r : one(Tv) + cscale[j] = r > zero(T) ? r : one(T) end map!(sqrt, cscale, cscale) @@ -306,7 +306,7 @@ include("dominated_column.jl") Perform post-solve. """ -function postsolve!(sol::Solution{Tv}, sol_::Solution{Tv}, ps::PresolveData{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, sol_::Solution{T}, ps::PresolveData{T}) where{T} # Check dimensions (sol_.m, sol_.n) == (ps.nrow, ps.ncol) || error( @@ -343,7 +343,7 @@ function postsolve!(sol::Solution{Tv}, sol_::Solution{Tv}, ps::PresolveData{Tv}) # Compute row primals for (i, row) in enumerate(ps.pb0.arows) - sol.Ax[i] = zero(Tv) + sol.Ax[i] = zero(T) for (j, aij) in zip(row.nzind, row.nzval) sol.Ax[i] += aij * sol.x[j] end @@ -359,7 +359,7 @@ end Perform pre-solve. """ -function presolve!(ps::PresolveData{Tv}) where{Tv<:Real} +function presolve!(ps::PresolveData{T}) where{T} # Check bound consistency on all rows/columns st = bounds_consistency_checks!(ps) @@ -487,7 +487,7 @@ Check that all primal & dual bounds are consistent. TODO: If not, declare primal/dual infeasibility and extract ray. """ -function bounds_consistency_checks!(ps::PresolveData{Tv}) where{Tv} +function bounds_consistency_checks!(ps::PresolveData{T}) where{T} # Check primal bounds for (i, (l, u)) in enumerate(zip(ps.lrow, ps.urow)) if ps.rowflag[i] && l > u @@ -499,21 +499,21 @@ function bounds_consistency_checks!(ps::PresolveData{Tv}) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Farkas ray: y⁺_i = y⁻_i = 1 (any > 0 value works) ps.solution.primal_status = Sln_Unknown ps.solution.dual_status = Sln_InfeasibilityCertificate ps.solution.is_primal_ray = false ps.solution.is_dual_ray = true - ps.solution.z_primal = ps.solution.z_dual = Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = T(Inf) i_ = ps.new_con_idx[i] - ps.solution.y_lower[i_] = one(Tv) - ps.solution.y_upper[i_] = one(Tv) + ps.solution.y_lower[i_] = one(T) + ps.solution.y_upper[i_] = one(T) return end @@ -528,21 +528,21 @@ function bounds_consistency_checks!(ps::PresolveData{Tv}) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Farkas ray: y⁺_i = y⁻_i = 1 (any > 0 value works) ps.solution.primal_status = Sln_Unknown ps.solution.dual_status = Sln_InfeasibilityCertificate ps.solution.is_primal_ray = false ps.solution.is_dual_ray = true - ps.solution.z_primal = ps.solution.z_dual = Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = T(Inf) j_ = ps.new_var_idx[j] - ps.solution.s_lower[j_] = one(Tv) - ps.solution.s_upper[j_] = one(Tv) + ps.solution.s_lower[j_] = one(T) + ps.solution.s_upper[j_] = one(T) return end @@ -561,7 +561,7 @@ Remove all empty rows. Called once at the beginning of the presolve procedure. If an empty row is created later, it is removed on the spot. """ -function remove_empty_rows!(ps::PresolveData{Tv}) where{Tv} +function remove_empty_rows!(ps::PresolveData{T}) where{T} nempty = 0 for i in 1:ps.pb0.ncon (ps.rowflag[i] && (ps.nzrow[i] == 0)) || continue @@ -580,7 +580,7 @@ Remove all empty columns. Called once at the beginning of the presolve procedure. If an empty column is created later, it is removed on the spot. """ -function remove_empty_columns!(ps::PresolveData{Tv}) where{Tv} +function remove_empty_columns!(ps::PresolveData{T}) where{T} for j in 1:ps.pb0.nvar remove_empty_column!(ps, j) ps.status == Trm_Unknown || break @@ -593,7 +593,7 @@ end Remove all fixed variables. """ -function remove_fixed_variables!(ps::PresolveData{Tv}) where{Tv} +function remove_fixed_variables!(ps::PresolveData{T}) where{T} for (j, flag) in enumerate(ps.colflag) flag || continue remove_fixed_variable!(ps, j) @@ -601,7 +601,7 @@ function remove_fixed_variables!(ps::PresolveData{Tv}) where{Tv} return nothing end -function remove_row_singletons!(ps::PresolveData{Tv}) where{Tv} +function remove_row_singletons!(ps::PresolveData{T}) where{T} nsingletons = 0 for i in ps.row_singletons remove_row_singleton!(ps, i) @@ -633,7 +633,7 @@ function remove_free_column_singletons!(ps::PresolveData) return nothing end -function remove_dominated_columns!(ps::PresolveData{Tv}) where{Tv} +function remove_dominated_columns!(ps::PresolveData{T}) where{T} # Strengthen dual bounds with column singletons for (j, (l, u)) in enumerate(zip(ps.lcol, ps.ucol)) (ps.colflag[j] && ps.nzcol[j] == 1) || continue @@ -641,7 +641,7 @@ function remove_dominated_columns!(ps::PresolveData{Tv}) where{Tv} col = ps.pb0.acols[j] # Find non-zero index nz = 0 - i, aij = 0, zero(Tv) + i, aij = 0, zero(T) for (i_, a_) in zip(col.nzind, col.nzval) if ps.rowflag[i_] && !iszero(a_) nz += 1; nz <= 1 || break @@ -664,7 +664,7 @@ function remove_dominated_columns!(ps::PresolveData{Tv}) where{Tv} # TODO elseif isfinite(l) && !isfinite(u) # Lower-bounded variable: `aij * yi ≤ cj` - if aij > zero(Tv) + if aij > zero(T) # yi ≤ cj / aij @debug "Col $j forces y$i <= $y_" ps.uy[i] = min(ps.uy[i], y_) @@ -676,7 +676,7 @@ function remove_dominated_columns!(ps::PresolveData{Tv}) where{Tv} elseif !isfinite(l) && isfinite(u) # Upper-bounded variable: `aij * yi ≥ cj` - if aij > zero(Tv) + if aij > zero(T) # yi ≥ cj / aij @debug "Col $j forces y$i >= $y_" ps.ly[i] = max(ps.ly[i], y_) diff --git a/src/Presolve/dominated_column.jl b/src/Presolve/dominated_column.jl index 0968c86c..d6d81a74 100644 --- a/src/Presolve/dominated_column.jl +++ b/src/Presolve/dominated_column.jl @@ -1,21 +1,21 @@ -struct DominatedColumn{Tv} <: PresolveTransformation{Tv} +struct DominatedColumn{T} <: PresolveTransformation{T} j::Int - x::Tv # Primal value - cj::Tv # Objective - col::Col{Tv} # Column + x::T # Primal value + cj::T # Objective + col::Col{T} # Column end -function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt(eps(Tv))) where{Tv} +function remove_dominated_column!(ps::PresolveData{T}, j::Int; tol::T=100*sqrt(eps(T))) where{T} ps.colflag[j] || return nothing # Compute implied bounds on reduced cost: `ls ≤ s ≤ us` - ls = us = zero(Tv) + ls = us = zero(T) col = ps.pb0.acols[j] for (i, aij) in zip(col.nzind, col.nzval) (ps.rowflag[i] && !iszero(aij)) || continue - ls += aij * ( (aij >= zero(Tv)) ? ps.ly[i] : ps.uy[i] ) - us += aij * ( (aij >= zero(Tv)) ? ps.uy[i] : ps.ly[i] ) + ls += aij * ( (aij >= zero(T)) ? ps.ly[i] : ps.uy[i] ) + us += aij * ( (aij >= zero(T)) ? ps.uy[i] : ps.ly[i] ) end # Check if column is dominated @@ -34,20 +34,20 @@ function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Unbounded ray: xj = -1 ps.solution.primal_status = Sln_InfeasibilityCertificate ps.solution.dual_status = Sln_Unknown ps.solution.is_primal_ray = true ps.solution.is_dual_ray = false - ps.solution.z_primal = ps.solution.z_dual = -Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = -T(Inf) j_ = ps.new_var_idx[j] - ps.solution.x[j_] = -one(Tv) + ps.solution.x[j_] = -one(T) return nothing end @@ -56,7 +56,7 @@ function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt ps.obj0 += cj * lb # Extract column and update rows - col_ = Col{Tv}(Int[], Tv[]) + col_ = Col{T}(Int[], T[]) for (i, aij) in zip(col.nzind, col.nzval) ps.rowflag[i] || continue @@ -91,20 +91,20 @@ function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt # Resize solution compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Unbounded ray: xj = -1 ps.solution.primal_status = Sln_InfeasibilityCertificate ps.solution.dual_status = Sln_Unknown ps.solution.is_primal_ray = true ps.solution.is_dual_ray = false - ps.solution.z_primal = ps.solution.z_dual = -Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = -T(Inf) j_ = ps.new_var_idx[j] - ps.solution.x[j_] = one(Tv) + ps.solution.x[j_] = one(T) return nothing end @@ -115,7 +115,7 @@ function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt ps.obj0 += cj * ub # Extract column and update rows - col_ = Col{Tv}(Int[], Tv[]) + col_ = Col{T}(Int[], T[]) for (i, aij) in zip(col.nzind, col.nzval) ps.rowflag[i] || continue @@ -140,12 +140,12 @@ function remove_dominated_column!(ps::PresolveData{Tv}, j::Int; tol::Tv=100*sqrt return nothing end -function postsolve!(sol::Solution{Tv}, op::DominatedColumn{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::DominatedColumn{T}) where{T} # Primal value sol.x[op.j] = op.x # Reduced cost - s = sol.is_dual_ray ? zero(Tv) : op.cj + s = sol.is_dual_ray ? zero(T) : op.cj for (i, aij) in zip(op.col.nzind, op.col.nzval) s -= aij * (sol.y_lower[i] - sol.y_upper[i]) end diff --git a/src/Presolve/empty_column.jl b/src/Presolve/empty_column.jl index f84cb03f..ec89edf4 100644 --- a/src/Presolve/empty_column.jl +++ b/src/Presolve/empty_column.jl @@ -1,10 +1,10 @@ -struct EmptyColumn{Tv} <: PresolveTransformation{Tv} +struct EmptyColumn{T} <: PresolveTransformation{T} j::Int # variable index - x::Tv # Variable value - s::Tv # Reduced cost + x::T # Variable value + s::T # Reduced cost end -function remove_empty_column!(ps::PresolveData{Tv}, j::Int) where{Tv} +function remove_empty_column!(ps::PresolveData{T}, j::Int) where{T} # Sanity check (ps.colflag[j] && (ps.nzcol[j] == 0)) || return nothing @@ -13,7 +13,7 @@ function remove_empty_column!(ps::PresolveData{Tv}, j::Int) where{Tv} cj = ps.obj[j] @debug "Removing empty column $j" cj lb ub - if cj > zero(Tv) + if cj > zero(T) if isfinite(lb) # Set variable to lower bound # Update objective constant @@ -28,24 +28,24 @@ function remove_empty_column!(ps::PresolveData{Tv}, j::Int) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Unbounded ray: xj = -1 ps.solution.primal_status = Sln_InfeasibilityCertificate ps.solution.dual_status = Sln_Unknown ps.solution.is_primal_ray = true ps.solution.is_dual_ray = false - ps.solution.z_primal = ps.solution.z_dual = -Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = -T(Inf) j_ = ps.new_var_idx[j] - ps.solution.x[j_] = -one(Tv) + ps.solution.x[j_] = -one(T) return nothing end - elseif cj < zero(Tv) + elseif cj < zero(T) if isfinite(ub) # Set variable to upper bound # Update objective constant @@ -60,32 +60,32 @@ function remove_empty_column!(ps::PresolveData{Tv}, j::Int) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Unbounded ray: xj = 1 ps.solution.primal_status = Sln_InfeasibilityCertificate ps.solution.dual_status = Sln_Unknown ps.solution.is_primal_ray = true ps.solution.is_dual_ray = false - ps.solution.z_primal = ps.solution.z_dual = -Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = -T(Inf) j_ = ps.new_var_idx[j] - ps.solution.x[j_] = one(Tv) + ps.solution.x[j_] = one(T) return end else # Any feasible value works if isfinite(lb) - push!(ps.ops, EmptyColumn(j, lb, zero(Tv))) + push!(ps.ops, EmptyColumn(j, lb, zero(T))) elseif isfinite(ub) - push!(ps.ops, EmptyColumn(j, ub, zero(Tv))) + push!(ps.ops, EmptyColumn(j, ub, zero(T))) else # Free variable with zero coefficient and empty column - push!(ps.ops, EmptyColumn(j, zero(Tv), zero(Tv))) + push!(ps.ops, EmptyColumn(j, zero(T), zero(T))) end end @@ -97,7 +97,7 @@ function remove_empty_column!(ps::PresolveData{Tv}, j::Int) where{Tv} return nothing end -function postsolve!(sol::Solution{Tv}, op::EmptyColumn{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::EmptyColumn{T}) where{T} sol.x[op.j] = op.x sol.s_lower[op.j] = pos_part(op.s) sol.s_upper[op.j] = neg_part(op.s) diff --git a/src/Presolve/empty_row.jl b/src/Presolve/empty_row.jl index 4f7f9c74..e9a9b5bd 100644 --- a/src/Presolve/empty_row.jl +++ b/src/Presolve/empty_row.jl @@ -1,19 +1,19 @@ # TODO: this is redundant with forcing constraint checks # => an empty row is automatically redundant or infeasible -struct EmptyRow{Tv} <: PresolveTransformation{Tv} +struct EmptyRow{T} <: PresolveTransformation{T} i::Int # row index - y::Tv # dual multiplier + y::T # dual multiplier end -function remove_empty_row!(ps::PresolveData{Tv}, i::Int) where{Tv} +function remove_empty_row!(ps::PresolveData{T}, i::Int) where{T} # Sanity checks (ps.rowflag[i] && ps.nzrow[i] == 0) || return nothing # Check bounds lb, ub = ps.lrow[i], ps.urow[i] - if ub < zero(Tv) + if ub < zero(T) # Infeasible @debug "Row $i is primal infeasible" ps.status = Trm_PrimalInfeasible @@ -22,22 +22,22 @@ function remove_empty_row!(ps::PresolveData{Tv}, i::Int) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Farkas ray: y⁺_i = 1 (any > 0 value works) ps.solution.primal_status = Sln_Unknown ps.solution.dual_status = Sln_InfeasibilityCertificate ps.solution.is_primal_ray = false ps.solution.is_dual_ray = true - ps.solution.z_primal = ps.solution.z_dual = Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = T(Inf) i_ = ps.new_con_idx[i] - ps.solution.y_upper[i] = one(Tv) + ps.solution.y_upper[i] = one(T) return - elseif lb > zero(Tv) + elseif lb > zero(T) @debug "Row $i is primal infeasible" ps.status = Trm_PrimalInfeasible ps.updated = true @@ -45,23 +45,23 @@ function remove_empty_row!(ps::PresolveData{Tv}, i::Int) where{Tv} # Resize problem compute_index_mapping!(ps) resize!(ps.solution, ps.nrow, ps.ncol) - ps.solution.x .= zero(Tv) - ps.solution.y_lower .= zero(Tv) - ps.solution.y_upper .= zero(Tv) - ps.solution.s_lower .= zero(Tv) - ps.solution.s_upper .= zero(Tv) + ps.solution.x .= zero(T) + ps.solution.y_lower .= zero(T) + ps.solution.y_upper .= zero(T) + ps.solution.s_lower .= zero(T) + ps.solution.s_upper .= zero(T) # Farkas ray: y⁺_i = 1 (any > 0 value works) ps.solution.primal_status = Sln_Unknown ps.solution.dual_status = Sln_InfeasibilityCertificate ps.solution.is_primal_ray = false ps.solution.is_dual_ray = true - ps.solution.z_primal = ps.solution.z_dual = Tv(Inf) + ps.solution.z_primal = ps.solution.z_dual = T(Inf) i_ = ps.new_con_idx[i] - ps.solution.y_lower[i_] = one(Tv) + ps.solution.y_lower[i_] = one(T) return else - push!(ps.ops, EmptyRow(i, zero(Tv))) + push!(ps.ops, EmptyRow(i, zero(T))) end # Book-keeping @@ -70,7 +70,7 @@ function remove_empty_row!(ps::PresolveData{Tv}, i::Int) where{Tv} ps.nrow -= 1 end -function postsolve!(sol::Solution{Tv}, op::EmptyRow{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::EmptyRow{T}) where{T} sol.y_lower[op.i] = pos_part(op.y) sol.y_upper[op.i] = neg_part(op.y) return nothing diff --git a/src/Presolve/fixed_variable.jl b/src/Presolve/fixed_variable.jl index 2a1c7046..655f263f 100644 --- a/src/Presolve/fixed_variable.jl +++ b/src/Presolve/fixed_variable.jl @@ -1,8 +1,8 @@ -struct FixedVariable{Tv} <: PresolveTransformation{Tv} +struct FixedVariable{T} <: PresolveTransformation{T} j::Int # variable index - x::Tv # primal value - c::Tv # current objective coeff - col::Col{Tv} # current column + x::T # primal value + c::T # current objective coeff + col::Col{T} # current column end function remove_fixed_variable!(ps::PresolveData, j::Int) @@ -56,9 +56,9 @@ function remove_fixed_variable!(ps::PresolveData, j::Int) return nothing end -function postsolve!(sol::Solution{Tv}, op::FixedVariable{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::FixedVariable{T}) where{T} sol.x[op.j] = op.x - s = sol.is_dual_ray ? zero(Tv) : op.c + s = sol.is_dual_ray ? zero(T) : op.c for (i, aij) in zip(op.col.nzind, op.col.nzval) s -= aij * (sol.y_lower[i] - sol.y_upper[i]) end diff --git a/src/Presolve/forcing_row.jl b/src/Presolve/forcing_row.jl index a77757aa..0d1dff6a 100644 --- a/src/Presolve/forcing_row.jl +++ b/src/Presolve/forcing_row.jl @@ -1,26 +1,26 @@ -struct ForcingRow{Tv} <: PresolveTransformation{Tv} +struct ForcingRow{T} <: PresolveTransformation{T} i::Int # Row index at_lower::Bool # Whether row is forced to its lower bound (false means upper) - row::Row{Tv} # Row - cols::Vector{Col{Tv}} # Columns of variables in forcing row - xs::Vector{Tv} # Primal values of variables in the row (upper or lower bound) - cs::Vector{Tv} # Objective coeffs of variables in forcing row + row::Row{T} # Row + cols::Vector{Col{T}} # Columns of variables in forcing row + xs::Vector{T} # Primal values of variables in the row (upper or lower bound) + cs::Vector{T} # Objective coeffs of variables in forcing row end -struct DominatedRow{Tv} <: PresolveTransformation{Tv} +struct DominatedRow{T} <: PresolveTransformation{T} i::Int # Row index end -function remove_forcing_row!(ps::PresolveData{Tv}, i::Int) where{Tv} +function remove_forcing_row!(ps::PresolveData{T}, i::Int) where{T} ps.rowflag[i] || return ps.nzrow[i] == 1 && return # skip row singletons # Implied row bounds row = ps.pb0.arows[i] - l_ = u_ = zero(Tv) + l_ = u_ = zero(T) for (j, aij) in zip(row.nzind, row.nzval) ps.colflag[j] || continue - if aij < zero(Tv) + if aij < zero(T) l_ += aij * ps.ucol[j] u_ += aij * ps.lcol[j] else @@ -39,7 +39,7 @@ function remove_forcing_row!(ps::PresolveData{Tv}, i::Int) where{Tv} ps.updated = true ps.nrow -= 1 - push!(ps.ops, DominatedRow{Tv}(i)) + push!(ps.ops, DominatedRow{T}(i)) # Update column non-zeros for (j, aij) in zip(row.nzind, row.nzval) ps.colflag[j] || continue @@ -54,15 +54,15 @@ function remove_forcing_row!(ps::PresolveData{Tv}, i::Int) where{Tv} [ j for (j, aij) in zip(row.nzind, row.nzval) if ps.colflag[j]], [aij for (j, aij) in zip(row.nzind, row.nzval) if ps.colflag[j]] ) - cols_ = Col{Tv}[] - xs = Tv[] - cs = Tv[] + cols_ = Col{T}[] + xs = T[] + cs = T[] for (j, aij) in zip(row.nzind, row.nzval) ps.colflag[j] || continue # Extract column j col = ps.pb0.acols[j] - col_ = Col{Tv}(Int[], Tv[]) + col_ = Col{T}(Int[], T[]) # Fix variable to its bound # TODO: put this in a function and mutualize with fixed variables @@ -116,15 +116,15 @@ function remove_forcing_row!(ps::PresolveData{Tv}, i::Int) where{Tv} [ j for (j, aij) in zip(row.nzind, row.nzval) if ps.colflag[j]], [aij for (j, aij) in zip(row.nzind, row.nzval) if ps.colflag[j]] ) - cols_ = Col{Tv}[] - xs = Tv[] - cs = Tv[] + cols_ = Col{T}[] + xs = T[] + cs = T[] for (j, aij) in zip(row.nzind, row.nzval) ps.colflag[j] || continue # Extract column j col = ps.pb0.acols[j] - col_ = Col{Tv}(Int[], Tv[]) + col_ = Col{T}(Int[], T[]) # Fix variable to its bound # TODO: put this in a function and mutualize with fixed variables @@ -175,14 +175,14 @@ function remove_forcing_row!(ps::PresolveData{Tv}, i::Int) where{Tv} return nothing end -function postsolve!(sol::Solution{Tv}, op::DominatedRow{Tv}) where{Tv} - sol.y_lower[op.i] = zero(Tv) - sol.y_upper[op.i] = zero(Tv) +function postsolve!(sol::Solution{T}, op::DominatedRow{T}) where{T} + sol.y_lower[op.i] = zero(T) + sol.y_upper[op.i] = zero(T) return nothing end # TODO: postsolve of forcing rows -function postsolve!(sol::Solution{Tv}, op::ForcingRow{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::ForcingRow{T}) where{T} # Primal for (j, xj) in zip(op.row.nzind, op.xs) diff --git a/src/Presolve/free_column_singleton.jl b/src/Presolve/free_column_singleton.jl index 98244ad3..21602b8b 100644 --- a/src/Presolve/free_column_singleton.jl +++ b/src/Presolve/free_column_singleton.jl @@ -1,14 +1,14 @@ -struct FreeColumnSingleton{Tv} <: PresolveTransformation{Tv} +struct FreeColumnSingleton{T} <: PresolveTransformation{T} i::Int # Row index j::Int # Column index - l::Tv # Row lower bound - u::Tv # Row upper bound - aij::Tv - y::Tv # Dual variable - row::Row{Tv} + l::T # Row lower bound + u::T # Row upper bound + aij::T + y::T # Dual variable + row::Row{T} end -function remove_free_column_singleton!(ps::PresolveData{Tv}, j::Int) where{Tv} +function remove_free_column_singleton!(ps::PresolveData{T}, j::Int) where{T} ps.colflag[j] && ps.nzcol[j] == 1 || return nothing # only column singletons @@ -17,7 +17,7 @@ function remove_free_column_singleton!(ps::PresolveData{Tv}, j::Int) where{Tv} # Find non-zero index # TODO: put this in a function nz = 0 - i, aij = 0, zero(Tv) + i, aij = 0, zero(T) for (i_, a_) in zip(col.nzind, col.nzval) if ps.rowflag[i_] nz += !iszero(a_); nz <= 1 || break @@ -41,7 +41,7 @@ function remove_free_column_singleton!(ps::PresolveData{Tv}, j::Int) where{Tv} l, u = ps.lcol[j], ps.ucol[j] if isfinite(l) || isfinite(u) # Not a free variable, compute implied bounds - if aij > zero(Tv) + if aij > zero(T) l_, u_ = lr, ur for (k, aik) in zip(row.nzind, row.nzval) (ps.colflag[k] && k != j) || continue @@ -84,8 +84,8 @@ function remove_free_column_singleton!(ps::PresolveData{Tv}, j::Int) where{Tv} y = ps.obj[j] / aij # dual of row i # Update objective - ps.obj0 += (y >= zero(Tv)) ? y * lr : y * ur - row_ = Row{Tv}(Int[], Tv[]) + ps.obj0 += (y >= zero(T)) ? y * lr : y * ur + row_ = Row{T}(Int[], T[]) for (j_, aij_) in zip(row.nzind, row.nzval) ps.colflag[j_] && (j_ != j) || continue @@ -108,16 +108,16 @@ function remove_free_column_singleton!(ps::PresolveData{Tv}, j::Int) where{Tv} return nothing end -function postsolve!(sol::Solution{Tv}, op::FreeColumnSingleton{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::FreeColumnSingleton{T}) where{T} # Dual y = op.y sol.y_lower[op.i] = pos_part(y) sol.y_upper[op.i] = neg_part(y) - sol.s_lower[op.j] = zero(Tv) - sol.s_upper[op.j] = zero(Tv) + sol.s_lower[op.j] = zero(T) + sol.s_upper[op.j] = zero(T) # Primal - sol.x[op.j] = sol.is_primal_ray ? zero(Tv) : (y >= zero(Tv) ? op.l : op.u) + sol.x[op.j] = sol.is_primal_ray ? zero(T) : (y >= zero(T) ? op.l : op.u) for (k, aik) in zip(op.row.nzind, op.row.nzval) sol.x[op.j] -= aik * sol.x[k] end diff --git a/src/Presolve/row_singleton.jl b/src/Presolve/row_singleton.jl index ad34a0a3..11b9ee2f 100644 --- a/src/Presolve/row_singleton.jl +++ b/src/Presolve/row_singleton.jl @@ -1,13 +1,13 @@ -struct RowSingleton{Tv} <: PresolveTransformation{Tv} +struct RowSingleton{T} <: PresolveTransformation{T} i::Int # Row index j::Int # Column index - aij::Tv # Row coefficient + aij::T # Row coefficient force_lower::Bool # Whether row was forcing the lower bound force_upper::Bool # Whether row was forcing the upper bound end -function remove_row_singleton!(ps::PresolveData{Tv}, i::Int) where{Tv} +function remove_row_singleton!(ps::PresolveData{T}, i::Int) where{T} # Sanity checks (ps.rowflag[i] && ps.nzrow[i] == 1) || return nothing @@ -15,7 +15,7 @@ function remove_row_singleton!(ps::PresolveData{Tv}, i::Int) where{Tv} # Find non-zero coefficient and its column index nz = 0 - j, aij = 0, zero(Tv) + j, aij = 0, zero(T) for (j_, aij_) in zip(row.nzind, row.nzval) if ps.colflag[j_] && !iszero(aij_) nz += 1; nz <= 1 || break # not a row singleton @@ -37,7 +37,7 @@ function remove_row_singleton!(ps::PresolveData{Tv}, i::Int) where{Tv} end # Compute implied bounds - if aij > zero(Tv) + if aij > zero(T) l = ps.lrow[i] / aij u = ps.urow[i] / aij else @@ -78,23 +78,23 @@ function remove_row_singleton!(ps::PresolveData{Tv}, i::Int) where{Tv} return nothing end -function postsolve!(sol::Solution{Tv}, op::RowSingleton{Tv}) where{Tv} +function postsolve!(sol::Solution{T}, op::RowSingleton{T}) where{T} if op.force_lower - if op.aij > zero(Tv) + if op.aij > zero(T) sol.y_lower[op.i] = sol.s_lower[op.j] / op.aij else sol.y_upper[op.i] = sol.s_lower[op.j] / abs(op.aij) end - sol.s_lower[op.j] = zero(Tv) + sol.s_lower[op.j] = zero(T) end if op.force_upper - if op.aij > zero(Tv) + if op.aij > zero(T) sol.y_upper[op.i] = sol.s_upper[op.j] / op.aij else sol.y_lower[op.i] = sol.s_upper[op.j] / abs(op.aij) end - sol.s_upper[op.j] = zero(Tv) + sol.s_upper[op.j] = zero(T) end return nothing diff --git a/src/model.jl b/src/model.jl index 4d6f7d86..05b57c5d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -1,7 +1,7 @@ -mutable struct Model{Tv} +mutable struct Model{T} # Parameters - params::Parameters{Tv} + params::Parameters{T} # TODO: model status #= Use an enum @@ -16,22 +16,22 @@ mutable struct Model{Tv} status::TerminationStatus # Problem data - pbdata::ProblemData{Tv} + pbdata::ProblemData{T} # Presolved problem # If presolved is disabled, this will point to m.pbdata - presolve_data::Union{Nothing, PresolveData{Tv}} + presolve_data::Union{Nothing, PresolveData{T}} # IPM solver # If required, the problem is transformed to standard form # when instantiating the IPMSolver object. - solver::Union{Nothing, AbstractIPMOptimizer{Tv}} + solver::Union{Nothing, AbstractIPMOptimizer{T}} # Problem solution (in original space) - solution::Union{Nothing, Solution{Tv}} + solution::Union{Nothing, Solution{T}} - Model{Tv}() where{Tv} = new{Tv}( - Parameters{Tv}(), Trm_NotCalled, ProblemData{Tv}(), + Model{T}() where{T} = new{T}( + Parameters{T}(), Trm_NotCalled, ProblemData{T}(), nothing, nothing, nothing ) end @@ -49,8 +49,8 @@ end import Base.empty! -function Base.empty!(m::Model{Tv}) where{Tv} - m.pbdata = ProblemData{Tv}() +function Base.empty!(m::Model{T}) where{T} + m.pbdata = ProblemData{T}() m.status = Trm_NotCalled m.presolve_data = nothing m.solver = nothing diff --git a/src/parameters.jl b/src/parameters.jl index 87bc8b3e..85ed0a3e 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -1,33 +1,33 @@ # TODO: enum for some parameters? # Or something easier to expand? """ - Parameters{Tv} + Parameters{T} """ -Base.@kwdef mutable struct Parameters{Tv} +Base.@kwdef mutable struct Parameters{T} # User limits BarrierIterationsLimit::Int = 100 TimeLimit::Float64 = Inf # Numerical tolerances - BarrierTolerancePFeas::Tv = sqrt(eps(Tv)) # primal feasibility - BarrierToleranceDFeas::Tv = sqrt(eps(Tv)) # dual feasibility - BarrierToleranceRGap::Tv = sqrt(eps(Tv)) # optimality - BarrierToleranceIFeas::Tv = sqrt(eps(Tv)) # infeasibility + BarrierTolerancePFeas::T = sqrt(eps(T)) # primal feasibility + BarrierToleranceDFeas::T = sqrt(eps(T)) # dual feasibility + BarrierToleranceRGap::T = sqrt(eps(T)) # optimality + BarrierToleranceIFeas::T = sqrt(eps(T)) # infeasibility # Algorithmic parameters BarrierAlgorithm::Int = 1 # TODO: docs BarrierCorrectionLimit::Int = 5 # Maximum number of centrality corrections - BarrierStepDampFactor::Tv = Tv(9_995 // 10_000) # Damp step size by this much - BarrierGammaMin::Tv = Tv(1 // 10) - BarrierCentralityOutlierThreshold::Tv = Tv(1 // 10) # Relative threshold for centrality outliers + BarrierStepDampFactor::T = T(9_995 // 10_000) # Damp step size by this much + BarrierGammaMin::T = T(1 // 10) + BarrierCentralityOutlierThreshold::T = T(1 // 10) # Relative threshold for centrality outliers # Linear algebra MatrixOptions::TLA.MatrixOptions = TLA.MatrixOptions(SparseMatrixCSC) - KKTOptions::KKT.SolverOptions = KKT.default_options(Tv) - BarrierPRegMin::Tv = sqrt(eps(Tv)) # primal - BarrierDRegMin::Tv = sqrt(eps(Tv)) # dual + KKTOptions::KKT.SolverOptions = KKT.default_options(T) + BarrierPRegMin::T = sqrt(eps(T)) # primal + BarrierDRegMin::T = sqrt(eps(T)) # dual # TODO: iterative refinement # I/O diff --git a/src/problemData.jl b/src/problemData.jl index c15891e2..abf2b201 100644 --- a/src/problemData.jl +++ b/src/problemData.jl @@ -1,17 +1,17 @@ using SparseArrays -mutable struct RowOrCol{Tv} +mutable struct RowOrCol{T} nzind::Vector{Int} - nzval::Vector{Tv} + nzval::Vector{T} end const Row = RowOrCol const Col = RowOrCol """ - ProblemData{Tv} + ProblemData{T} -Data structure for storing problem data in precision `Tv`. +Data structure for storing problem data in precision `T`. The LP is represented in canonical form @@ -23,7 +23,7 @@ The LP is represented in canonical form \\end{array} ``` """ -mutable struct ProblemData{Tv} +mutable struct ProblemData{T} name::String @@ -34,42 +34,42 @@ mutable struct ProblemData{Tv} # Objective # TODO: objective sense objsense::Bool # true is min, false is max - obj::Vector{Tv} - obj0::Tv # Constant objective offset + obj::Vector{T} + obj0::T # Constant objective offset # Constraint matrix # We store both rows and columns. It is redundant but simplifies access. # TODO: put this in its own data structure? (would allow more flexibility in modelling) - arows::Vector{Row{Tv}} - acols::Vector{Col{Tv}} + arows::Vector{Row{T}} + acols::Vector{Col{T}} # TODO: Data structures for QP # qrows # qcols # Bounds - lcon::Vector{Tv} - ucon::Vector{Tv} - lvar::Vector{Tv} - uvar::Vector{Tv} + lcon::Vector{T} + ucon::Vector{T} + lvar::Vector{T} + uvar::Vector{T} # Names con_names::Vector{String} var_names::Vector{String} # Only allow empty problems to be instantiated for now - ProblemData{Tv}(pbname::String="") where {Tv} = new{Tv}( + ProblemData{T}(pbname::String="") where {T} = new{T}( pbname, 0, 0, - true, Tv[], zero(Tv), - Row{Tv}[], Col{Tv}[], - Tv[], Tv[], Tv[], Tv[], + true, T[], zero(T), + Row{T}[], Col{T}[], + T[], T[], T[], T[], String[], String[] ) end import Base.empty! -function Base.empty!(pb::ProblemData{Tv}) where{Tv} +function Base.empty!(pb::ProblemData{T}) where{T} pb.name = "" @@ -77,16 +77,16 @@ function Base.empty!(pb::ProblemData{Tv}) where{Tv} pb.nvar = 0 pb.objsense = true - pb.obj = Tv[] - pb.obj0 = zero(Tv) + pb.obj = T[] + pb.obj0 = zero(T) - pb.arows = Row{Tv}[] - pb.acols = Col{Tv}[] + pb.arows = Row{T}[] + pb.acols = Col{T}[] - pb.lcon = Tv[] - pb.ucon = Tv[] - pb.lvar = Tv[] - pb.uvar = Tv[] + pb.lcon = T[] + pb.ucon = T[] + pb.lvar = T[] + pb.uvar = T[] pb.con_names = String[] pb.var_names = String[] @@ -126,20 +126,20 @@ end Add one linear constraint to the problem. # Arguments -* `pb::ProblemData{Tv}`: the problem to which the new row is added +* `pb::ProblemData{T}`: the problem to which the new row is added * `rind::Vector{Int}`: column indices in the new row -* `rval::Vector{Tv}`: non-zero values in the new row -* `l::Tv` -* `u::Tv` +* `rval::Vector{T}`: non-zero values in the new row +* `l::T` +* `u::T` * `name::String`: row name (defaults to `""`) * `issorted::Bool`: indicates whether the row indices are already issorted. """ -function add_constraint!(pb::ProblemData{Tv}, - rind::Vector{Int}, rval::Vector{Tv}, - l::Tv, u::Tv, +function add_constraint!(pb::ProblemData{T}, + rind::Vector{Int}, rval::Vector{T}, + l::T, u::T, name::String=""; issorted::Bool=false -)::Int where{Tv} +)::Int where{T} # Sanity checks nz = length(rind) nz == length(rval) || throw(DimensionMismatch( @@ -154,7 +154,7 @@ function add_constraint!(pb::ProblemData{Tv}, if nz == 0 # emtpy row - push!(pb.arows, Row{Tv}(Int[], Tv[])) + push!(pb.arows, Row{T}(Int[], T[])) return pb.ncon end @@ -166,11 +166,11 @@ function add_constraint!(pb::ProblemData{Tv}, # Create new row if issorted - row = Row{Tv}(copy(rind), copy(rval)) + row = Row{T}(copy(rind), copy(rval)) else # Sort indices first p = sortperm(rind) - row = Row{Tv}(rind[p], rval[p]) + row = Row{T}(rind[p], rval[p]) end push!(pb.arows, row) @@ -192,21 +192,21 @@ end Add one variable to the problem. # Arguments -* `pb::ProblemData{Tv}`: the problem to which the new column is added +* `pb::ProblemData{T}`: the problem to which the new column is added * `cind::Vector{Int}`: row indices in the new column -* `cval::Vector{Tv}`: non-zero values in the new column -* `obj::Tv`: objective coefficient -* `l::Tv`: column lower bound -* `u::Tv`: column upper bound +* `cval::Vector{T}`: non-zero values in the new column +* `obj::T`: objective coefficient +* `l::T`: column lower bound +* `u::T`: column upper bound * `name::String`: column name (defaults to `""`) * `issorted::Bool`: indicates whether the column indices are already issorted. """ -function add_variable!(pb::ProblemData{Tv}, - cind::Vector{Int}, cval::Vector{Tv}, - obj::Tv, l::Tv, u::Tv, +function add_variable!(pb::ProblemData{T}, + cind::Vector{Int}, cval::Vector{T}, + obj::T, l::T, u::T, name::String=""; issorted::Bool=false -)::Int where{Tv} +)::Int where{T} # Sanity checks nz = length(cind) nz == length(cval) || throw(DimensionMismatch( @@ -221,7 +221,7 @@ function add_variable!(pb::ProblemData{Tv}, push!(pb.obj, obj) if nz == 0 - push!(pb.acols, Col{Tv}(Int[], Tv[])) + push!(pb.acols, Col{T}(Int[], T[])) return pb.nvar # empty column end @@ -229,11 +229,11 @@ function add_variable!(pb::ProblemData{Tv}, # Create a new column if issorted - col = Col{Tv}(copy(cind), copy(cval)) + col = Col{T}(copy(cind), copy(cval)) else # Sort indices p = sortperm(cind) - col = Col{Tv}(cind[p], cind[p]) + col = Col{T}(cind[p], cind[p]) end push!(pb.acols, col) @@ -254,14 +254,14 @@ end Load entire problem. """ -function load_problem!(pb::ProblemData{Tv}, +function load_problem!(pb::ProblemData{T}, name::String, - objsense::Bool, obj::Vector{Tv}, obj0::Tv, + objsense::Bool, obj::Vector{T}, obj0::T, A::SparseMatrixCSC, - lcon::Vector{Tv}, ucon::Vector{Tv}, - lvar::Vector{Tv}, uvar::Vector{Tv}, + lcon::Vector{T}, ucon::Vector{T}, + lvar::Vector{T}, uvar::Vector{T}, con_names::Vector{String}, var_names::Vector{String} -) where{Tv} +) where{T} empty!(pb) # Sanity checks @@ -289,17 +289,17 @@ function load_problem!(pb::ProblemData{Tv}, pb.var_names = copy(var_names) # Load coefficients - pb.acols = Vector{Col{Tv}}(undef, nvar) - pb.arows = Vector{Row{Tv}}(undef, ncon) + pb.acols = Vector{Col{T}}(undef, nvar) + pb.arows = Vector{Row{T}}(undef, ncon) for j in 1:nvar col = A[:, j] - pb.acols[j] = Col{Tv}(col.nzind, col.nzval) + pb.acols[j] = Col{T}(col.nzind, col.nzval) end At = sparse(A') for i in 1:ncon row = At[:, i] - pb.arows[i] = Row{Tv}(row.nzind, row.nzval) + pb.arows[i] = Row{T}(row.nzind, row.nzval) end return pb @@ -314,7 +314,7 @@ end Delete a single constraint from problem `pb`. """ -function delete_constraint!(pb::ProblemData{Tv}, rind::Int) where{Tv} +function delete_constraint!(pb::ProblemData{T}, rind::Int) where{T} # Sanity checks 1 <= rind <= pb.ncon || error("Invalid row index $rind") @@ -359,7 +359,7 @@ Delete rows in collection `rind` from problem `pb`. * `pb::ProblemData` * `rinds`: collection of row indices to be removed """ -function delete_constraints!(pb::ProblemData{Tv}, rinds) where{Tv} +function delete_constraints!(pb::ProblemData{T}, rinds) where{T} # TODO: don't use fallback for i in rinds delete_constraint!(pb, i) @@ -372,7 +372,7 @@ end Delete a single column from problem `pb`. """ -function delete_variable!(pb::ProblemData{Tv}, cind::Int) where{Tv} +function delete_variable!(pb::ProblemData{T}, cind::Int) where{T} # Sanity checks 1 <= cind <= pb.nvar || error("Invalid column index $cind") @@ -418,7 +418,7 @@ Delete a collection of columns from problem `pb`. * `pb::ProblemData` * `cinds`: collection of row indices to be removed """ -function delete_variables!(pb::ProblemData{Tv}, cinds) where{Tv} +function delete_variables!(pb::ProblemData{T}, cinds) where{T} # TODO: don't use fallback for j in cinds delete_variable!(pb, j) @@ -432,12 +432,12 @@ end Set the coefficient `(i, j)` to value `v`. # Arguments -* `pb::ProblemData{Tv}`: the problem whose coefficient +* `pb::ProblemData{T}`: the problem whose coefficient * `i::Int`: row index * `j::Int`: column index -* `v::Tv`: coefficient value +* `v::T`: coefficient value """ -function set_coefficient!(pb::ProblemData{Tv}, i::Int, j::Int, v::Tv) where{Tv} +function set_coefficient!(pb::ProblemData{T}, i::Int, j::Int, v::T) where{T} # Sanity checks 1 <= i <= pb.ncon && 1 <= j <= pb.nvar || error( "Cannot access coeff $((i, j)) in a model of size ($(pb.ncon), $(pb.nvar))" @@ -451,11 +451,11 @@ function set_coefficient!(pb::ProblemData{Tv}, i::Int, j::Int, v::Tv) where{Tv} end """ - _set_coefficient!(roc::RowOrCol{Tv}, ind::Int, v::Tv) + _set_coefficient!(roc::RowOrCol{T}, ind::Int, v::T) Set coefficient to value `v`. """ -function _set_coefficient!(roc::RowOrCol{Tv}, ind::Int, v::Tv) where{Tv} +function _set_coefficient!(roc::RowOrCol{T}, ind::Int, v::T) where{T} # Check if index already exists k = searchsortedfirst(roc.nzind, ind) diff --git a/src/solution.jl b/src/solution.jl index 1266c612..bfd9b64c 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -1,4 +1,4 @@ -mutable struct Solution{Tv} +mutable struct Solution{T} m::Int n::Int @@ -7,22 +7,22 @@ mutable struct Solution{Tv} is_primal_ray::Bool is_dual_ray::Bool - z_primal::Tv - z_dual::Tv + z_primal::T + z_dual::T - x::Vector{Tv} - Ax::Vector{Tv} - y_lower::Vector{Tv} - y_upper::Vector{Tv} - s_lower::Vector{Tv} - s_upper::Vector{Tv} + x::Vector{T} + Ax::Vector{T} + y_lower::Vector{T} + y_upper::Vector{T} + s_lower::Vector{T} + s_upper::Vector{T} - Solution{Tv}(m, n) where{Tv} = new{Tv}( + Solution{T}(m, n) where{T} = new{T}( m, n, Sln_Unknown, Sln_Unknown, false, false, - zero(Tv), zero(Tv), - zeros(Tv, n), zeros(Tv, m), - zeros(Tv, m), zeros(Tv, m), - zeros(Tv, n), zeros(Tv, n) + zero(T), zero(T), + zeros(T, n), zeros(T, m), + zeros(T, m), zeros(T, m), + zeros(T, n), zeros(T, n) ) end diff --git a/test/KKT/krylov.jl b/test/KKT/krylov.jl index be643be0..7c038329 100644 --- a/test/KKT/krylov.jl +++ b/test/KKT/krylov.jl @@ -10,7 +10,7 @@ import Krylov 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.KrylovSPDSolver(f, A) + kkt = KKT.KrylovSPD(f, A) KKT.run_ls_tests(A, kkt) end end @@ -21,7 +21,7 @@ import Krylov 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.KrylovSIDSolver(f, A) + kkt = KKT.KrylovSID(f, A) KKT.run_ls_tests(A, kkt) end end @@ -32,7 +32,7 @@ import Krylov 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.KrylovSQDSolver(f, A) + kkt = KKT.KrylovSQD(f, A) KKT.run_ls_tests(A, kkt) end end diff --git a/test/KKT/lapack.jl b/test/KKT/lapack.jl index 16f6bcc7..d58b9a98 100644 --- a/test/KKT/lapack.jl +++ b/test/KKT/lapack.jl @@ -1,11 +1,11 @@ @testset "LAPACK" begin - for Tv in TvTYPES - @testset "$Tv" begin - A = Matrix{Tv}([ + for T in TvTYPES + @testset "$T" begin + A = Matrix{T}([ 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.Dense_SymPosDef(A) + kkt = KKT.DenseSPD(A) KKT.run_ls_tests(A, kkt) end end diff --git a/test/KKT/ldlfact.jl b/test/KKT/ldlfact.jl index 312a3e69..c5ddea23 100644 --- a/test/KKT/ldlfact.jl +++ b/test/KKT/ldlfact.jl @@ -1,11 +1,11 @@ @testset "LDLFact" begin - for Tv in TvTYPES - @testset "$Tv" begin - A = SparseMatrixCSC{Tv, Int}([ + for T in TvTYPES + @testset "$T" begin + A = SparseMatrixCSC{T, Int}([ 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.LDLFact_SymQuasDef(A) + kkt = KKT.LDLFactSQD(A) KKT.run_ls_tests(A, kkt) end end diff --git a/test/Presolve/empty_column.jl b/test/Presolve/empty_column.jl index 9f61e880..de581b50 100644 --- a/test/Presolve/empty_column.jl +++ b/test/Presolve/empty_column.jl @@ -1,4 +1,4 @@ -function emtpy_column_tests(Tv::Type) +function emtpy_column_tests(T::Type) # We test all the following combinations: #= @@ -17,18 +17,18 @@ function emtpy_column_tests(Tv::Type) ------------------------------ =# function build_problem(l, u, c) - pb = Tulip.ProblemData{Tv}() + pb = Tulip.ProblemData{T}() Tulip.load_problem!(pb, "Test", - true, [c], zero(Tv), - spzeros(Tv, 0, 1), Tv[], Tv[], [l], [u], String[], ["x"] + true, [c], zero(T), + spzeros(T, 0, 1), T[], T[], [l], [u], String[], ["x"] ) return pb end - L = Tv.([-Inf, -1]) - U = Tv.([ 1, Inf]) - C = Tv.([-1, 0, 1]) + L = T.([-Inf, -1]) + U = T.([ 1, Inf]) + C = T.([-1, 0, 1]) for l in L, u in U, c in C @testset "$((l, u, c))" begin pb = build_problem(l, u, c) @@ -75,7 +75,7 @@ function emtpy_column_tests(Tv::Type) end @testset "Empty column" begin - for Tv in TvTYPES - @testset "$Tv" begin emtpy_column_tests(Tv) end + for T in TvTYPES + @testset "$T" begin emtpy_column_tests(T) end end end \ No newline at end of file diff --git a/test/Presolve/empty_row.jl b/test/Presolve/empty_row.jl index c15cb863..8485341d 100644 --- a/test/Presolve/empty_row.jl +++ b/test/Presolve/empty_row.jl @@ -1,4 +1,4 @@ -function empty_row_tests(Tv::Type) +function empty_row_tests(T::Type) # Build the following model #= @@ -6,17 +6,17 @@ function empty_row_tests(Tv::Type) s.t. -1 ⩽ 0 * x + 0 * y + 0 * z ⩽ 1 1 ⩽ 0 * x + 0 * y + 0 * z ⩽ 2 =# - pb = Tulip.ProblemData{Tv}() + pb = Tulip.ProblemData{T}() m, n = 2, 3 - A = spzeros(Tv, m, n) + A = spzeros(T, m, n) - b = ones(Tv, m) - c = ones(Tv, n) + b = ones(T, m) + c = ones(T, n) Tulip.load_problem!(pb, "test", - true, c, zero(Tv), - A, Tv.([-1, 1]), Tv.([1, 2]), zeros(Tv, n), fill(Tv(Inf), n), + true, c, zero(T), + A, T.([-1, 1]), T.([1, 2]), zeros(T, n), fill(T(Inf), n), ["c1", "c2"], ["x", "y", "z"] ) @@ -35,7 +35,7 @@ function empty_row_tests(Tv::Type) @test length(ps.ops) == 1 op = ps.ops[1] - @test isa(op, Tulip.EmptyRow{Tv}) + @test isa(op, Tulip.EmptyRow{T}) @test op.i == 1 @test iszero(op.y) @@ -51,35 +51,35 @@ function empty_row_tests(Tv::Type) # Check solution status & objective value sol = ps.solution @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate - @test sol.z_primal == sol.z_dual == Tv(Inf) + @test sol.z_primal == sol.z_dual == T(Inf) # Check Farkas ray # (current problem only has 1 row) - @test sol.y_lower[1] > zero(Tv) + @test sol.y_lower[1] > zero(T) - test_empty_row_1(Tv) - test_empty_row_2(Tv) + test_empty_row_1(T) + test_empty_row_2(T) return end -function test_empty_row_1(Tv::Type) +function test_empty_row_1(T::Type) # Empty row with l > 0 #= min x s.t. 1 ⩽ 0 * x ⩽ 2 x >= 0 =# - pb = Tulip.ProblemData{Tv}() + pb = Tulip.ProblemData{T}() m, n = 1, 1 - A = spzeros(Tv, m, n) - c = ones(Tv, n) + A = spzeros(T, m, n) + c = ones(T, n) Tulip.load_problem!(pb, "test", - true, c, zero(Tv), - A, Tv.([1]), Tv.([2]), zeros(Tv, n), fill(Tv(Inf), n), + true, c, zero(T), + A, T.([1]), T.([2]), zeros(T, n), fill(T(Inf), n), ["c1"], ["x"] ) @@ -94,31 +94,31 @@ function test_empty_row_1(Tv::Type) # Check solution status & objective value sol = ps.solution @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate - @test sol.z_primal == sol.z_dual == Tv(Inf) + @test sol.z_primal == sol.z_dual == T(Inf) # Check Farkas ray # (current problem only has 1 row) - @test sol.y_lower[1] > zero(Tv) + @test sol.y_lower[1] > zero(T) return nothing end -function test_empty_row_2(Tv::Type) +function test_empty_row_2(T::Type) # Empty row with u < 0 #= min x s.t. -2 ⩽ 0 * x ⩽ -1 x >= 0 =# - pb = Tulip.ProblemData{Tv}() + pb = Tulip.ProblemData{T}() m, n = 1, 1 - A = spzeros(Tv, m, n) - c = ones(Tv, n) + A = spzeros(T, m, n) + c = ones(T, n) Tulip.load_problem!(pb, "test", - true, c, zero(Tv), - A, Tv.([-2]), Tv.([-1]), zeros(Tv, n), fill(Tv(Inf), n), + true, c, zero(T), + A, T.([-2]), T.([-1]), zeros(T, n), fill(T(Inf), n), ["c1"], ["x"] ) @@ -133,17 +133,17 @@ function test_empty_row_2(Tv::Type) # Check solution status & objective value sol = ps.solution @test sol.dual_status == Tulip.Sln_InfeasibilityCertificate - @test sol.z_primal == sol.z_dual == Tv(Inf) + @test sol.z_primal == sol.z_dual == T(Inf) # Check Farkas ray # (current problem only has 1 row) - @test sol.y_upper[1] > zero(Tv) + @test sol.y_upper[1] > zero(T) return nothing end @testset "Empty row" begin - for Tv in TvTYPES - @testset "$Tv" begin empty_row_tests(Tv) end + for T in TvTYPES + @testset "$T" begin empty_row_tests(T) end end end \ No newline at end of file diff --git a/test/examples.jl b/test/examples.jl index 4fd13806..f9c4ac66 100644 --- a/test/examples.jl +++ b/test/examples.jl @@ -2,19 +2,19 @@ const examples_dir = joinpath(@__FILE__, "../../examples") @testset "Optimal" begin include(joinpath(examples_dir, "optimal.jl")) - for Tv in TvTYPES - @testset "$Tv" begin ex_optimal(Tv) end + for T in TvTYPES + @testset "$T" begin ex_optimal(T) end end end @testset "PrimalInfeas" begin include(joinpath(examples_dir, "infeasible.jl")) - for Tv in TvTYPES - @testset "$Tv" begin ex_infeasible(Tv) end + for T in TvTYPES + @testset "$T" begin ex_infeasible(T) end end end @testset "DualInfeas" begin include(joinpath(examples_dir, "unbounded.jl")) - for Tv in TvTYPES - @testset "$Tv" begin ex_unbounded(Tv) end + for T in TvTYPES + @testset "$T" begin ex_unbounded(T) end end end