diff --git a/.gitignore b/.gitignore index 9db7b5ce..98fe4c05 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,9 @@ *.jl.*.cov *.jl.mem Manifest.toml + +/app +/dat +/exp +/.vscode +/*.jl \ No newline at end of file diff --git a/Project.toml b/Project.toml index 57cdd75b..51cad6e9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ name = "Tulip" uuid = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa" authors = ["Mathieu Tanneau "] -version = "0.7.5" +version = "0.8.0" [deps] CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd" CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" @@ -21,6 +22,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] CodecBzip2 = "0.7.2" CodecZlib = "0.7.0" +Krylov = "0.7.3" LDLFactorizations = "0.6, 0.7, 0.8" LinearOperators = "2.0" MathOptInterface = "0.9.5" diff --git a/docs/make.jl b/docs/make.jl index f7430d5c..2675c39a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,9 +1,12 @@ using Documenter, Tulip +const _FAST = findfirst(isequal("--fast"), ARGS) !== nothing + makedocs( sitename = "Tulip.jl", authors = "Mathieu Tanneau", format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"), + doctest = !_FAST, pages = [ "Home" => "index.md", "Tutorials" => Any[ @@ -11,7 +14,6 @@ makedocs( ], "User manual" => Any[ "Problem formulation" => "manual/formulation.md", - "Solving linear systems" => "manual/linear_systems.md", "Algorithms" => Any[ "Homogeneous Self-Dual" => "manual/IPM/HSD.md", "Predictor-Corrector" => "manual/IPM/MPC.md" @@ -20,7 +22,13 @@ makedocs( ], "Reference" => Any[ "Presolve" => "reference/Presolve/presolve.md", - "KKT solvers" => "reference/KKT/kkt_solvers.md", + "KKT" => [ + "reference/KKT/kkt.md", + "reference/KKT/tlp_cholmod.md", + "reference/KKT/tlp_dense.md", + "reference/KKT/tlp_krylov.md", + "reference/KKT/tlp_ldlfact.md", + ], "Julia API" => "reference/API.md", "Attributes" => "reference/attributes.md", "Options" => "reference/options.md" @@ -30,4 +38,4 @@ makedocs( deploydocs( repo = "github.com/ds4dm/Tulip.jl.git" -) \ No newline at end of file +) diff --git a/docs/src/manual/linear_systems.md b/docs/src/manual/linear_systems.md deleted file mode 100644 index 8b2b3905..00000000 --- a/docs/src/manual/linear_systems.md +++ /dev/null @@ -1,62 +0,0 @@ -```@meta -CurrentModule = Tulip.KKT -``` - -# Solving linear systems - -The interior-point algorithm in Tulip requires the solution, at each iteration, of the following symmetric quasi-definite _augmented system_ -```math -\begin{bmatrix} - -(\Theta^{-1} + R_{p}) & A^{T}\\ - A & R_{d} -\end{bmatrix} -\begin{bmatrix} - \Delta x\\ - \Delta y -\end{bmatrix} -= -\begin{bmatrix} - \xi_d\\ - \xi_p -\end{bmatrix} -``` -where -* ``\Delta x, \Delta y`` are primal and dual search directions, -* ``A`` is the problem's constraint matrix, -* ``\Theta``, ``R_p`` and ``R_d`` are positive diagonal matrices, -* ``\xi_p, \xi_d`` are right-hand side vectors. - -The augmented system above can be reduced to the positive-definite _normal equations system_ -```math -\begin{array}{rl} -\left( - A (\Theta^{-1} + R_{p})^{-1} A^{T} + R_{d} -\right) -\Delta y -& = -\xi_{p} + A (Θ^{-1} + R_{p})^{-1} \xi_{d}\\ -\Delta x &= (Θ^{-1} + R_{p})^{-1} (A^{T} \Delta y - \xi_{d}) -\end{array} -``` -If available and when selected, this reduction is transparent to the interior-point algorithm. - -To enable the use of fast external libraries and/or specialized routines, the resolution of linear systems is performed by an [`AbstractKKTSolver`] object. - - -## Supported linear solvers - -Here is a list of currently supported linear solvers: - -| Linear solver type | System | Backend | Method | -|:-------------------|:-------|:--------|:-------| -| [`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/kkt.md b/docs/src/reference/KKT/kkt.md new file mode 100644 index 00000000..b9151408 --- /dev/null +++ b/docs/src/reference/KKT/kkt.md @@ -0,0 +1,77 @@ +```@meta +CurrentModule = Tulip.KKT +``` + +# Overview + +The `KKT` module provides a modular, customizable interface for developing and selecting various approaches to solve the KKT systems. + +## KKT backends + +```@docs +AbstractKKTBackend +``` + +```@docs +DefaultKKTBackend +``` + +## KKT systems + +All formulations below refer to a linear program in primal-dual standard form +```math + \begin{array}{rl} + (P) \ \ \ + \displaystyle \min_{x} + & c^{\top} x\\ + s.t. + & A x = b\\ + & l \leq x \leq u + \end{array} + \quad \quad \quad + \begin{array}{rl} + (D) \ \ \ + \displaystyle \max_{y, z} + & b^{\top} y + l^{\top}z^{l} - u^{\top}z^{u}\\ + s.t. + & A^{\top}y + z^{l} - z^{u} = c\\ + & z^{l}, z^{u} \geq 0 + \end{array} +``` + +```@docs +AbstractKKTSystem +``` + +```@docs +DefaultKKTSystem +``` + +```@docs +K2 +``` + +```@docs +K1 +``` + +## KKT solvers + +```@docs +AbstractKKTSolver +``` + +Custom linear solvers should (preferably) inherit from the `AbstractKKTSolver` class, +and extend the following functions: + +```@docs +setup +``` + +```@docs +update! +``` + +```@docs +solve! +``` \ No newline at end of file diff --git a/docs/src/reference/KKT/kkt_solvers.md b/docs/src/reference/KKT/kkt_solvers.md deleted file mode 100644 index eb74b6ea..00000000 --- a/docs/src/reference/KKT/kkt_solvers.md +++ /dev/null @@ -1,75 +0,0 @@ -```@meta -CurrentModule = Tulip.KKT -``` - -## Overview - -### AbstractKKTSolver - -This is the base type from which all implementations should derive. - -```@docs -AbstractKKTSolver -``` - -Custom linear solvers should inherit from the `AbstractKKTSolver` class, -and extend the following two functions: - -```@docs -update! -``` - -```@docs -solve! -``` - -### Choosing between linear solvers - -```@docs -KKTOptions -``` - -## Dense/LAPACK - -```@docs -DenseSPD -``` - -## CHOLMOD - -```@docs -CholmodSolver -``` - -```@docs -CholmodSQD -``` - -```@docs -CholmodSPD -``` - -## LDLFactorizations - -```@docs -LDLFactSQD -``` - -## Krylov - -!!! warning - Iterative methods are still an experimental feature. - Some numerical and performance issues should be expected. - - -```@docs -KrylovSPD -``` - -```@docs -KrylovSID -``` - -```@docs -KrylovSQD -``` \ No newline at end of file diff --git a/docs/src/reference/KKT/tlp_cholmod.md b/docs/src/reference/KKT/tlp_cholmod.md new file mode 100644 index 00000000..1e78f6ad --- /dev/null +++ b/docs/src/reference/KKT/tlp_cholmod.md @@ -0,0 +1,13 @@ +```@meta +CurrentModule = Tulip.KKT +``` + +# TlpCholmod + +```@docs +TlpCholmod.Backend +``` + +```@docs +TlpCholmod.CholmodSolver +``` \ No newline at end of file diff --git a/docs/src/reference/KKT/tlp_dense.md b/docs/src/reference/KKT/tlp_dense.md new file mode 100644 index 00000000..05e18660 --- /dev/null +++ b/docs/src/reference/KKT/tlp_dense.md @@ -0,0 +1,13 @@ +```@meta +CurrentModule = Tulip.KKT.TlpDense +``` + +# TlpDense + +```@docs +TlpDense.Backend +``` + +```@docs +TlpDense.DenseSolver +``` \ No newline at end of file diff --git a/docs/src/reference/KKT/tlp_krylov.md b/docs/src/reference/KKT/tlp_krylov.md new file mode 100644 index 00000000..80c926f5 --- /dev/null +++ b/docs/src/reference/KKT/tlp_krylov.md @@ -0,0 +1,30 @@ +```@meta +CurrentModule = Tulip.KKT.TlpKrylov +``` + +# TlpKrylov + +!!! warning + Iterative methods are still an experimental feature. + Some numerical and performance issues should be expected. + + +```@docs +TlpKrylov.Backend +``` + +```@docs +TlpKrylov.AbstractKrylovSolver +``` + +```@docs +TlpKrylov.SPDSolver +``` + +```@docs +TlpKrylov.SIDSolver +``` + +```@docs +TlpKrylov.SQDSolver +``` \ No newline at end of file diff --git a/docs/src/reference/KKT/tlp_ldlfact.md b/docs/src/reference/KKT/tlp_ldlfact.md new file mode 100644 index 00000000..5e1b3e07 --- /dev/null +++ b/docs/src/reference/KKT/tlp_ldlfact.md @@ -0,0 +1,13 @@ +```@meta +CurrentModule = Tulip.KKT.TlpLDLFactorizations +``` + +# TlpLDLFactorizations + +```@docs +TlpLDLFactorizations.Backend +``` + +```@docs +TlpLDLFactorizations.LDLFactSolver +``` \ No newline at end of file diff --git a/docs/src/reference/options.md b/docs/src/reference/options.md index 7d0afc02..7047026a 100644 --- a/docs/src/reference/options.md +++ b/docs/src/reference/options.md @@ -13,28 +13,7 @@ For instance, in double-precision arithmetic, i.e., `Tv=Float64`, we have ``\eps Factory ``` -## Presolve - -These parameters control the execution of the presolve phase. -They should be called as `"Presolve_"`. - - -## Linear Algebra - -These parameters control the linear algebra implementation - -| Parameter | Description | Default | -|:----------|:------------|:--------| -| `MatrixFactory` | See [`Factory`](@ref) | `Factory(SparseMatrixCSC)` - - -## KKT solvers - -| Parameter | Description | Default | -|:----------|:------------|:--------| -| `Factory` | See [`Factory`](@ref) | [`KKT.CholmodSQD`](@ref) for `Float64`, [`KKT.LDLFactSQD`](@ref) otherwise | - -## Interior-Point +## IPM ### Tolerances @@ -66,6 +45,33 @@ Numerical tolerances for the interior-point algorithm. | `IterationsLimit` | Maximum number of barrier iterations | `100` | | `TimeLimit` | Time limit, in seconds | `Inf` | + + + + + + + +## KKT + +| Parameter | Description | Default | +|:----------|:------------|:--------| +| `Backend` | See [KKT backends](@ref) | [`KKT.DefaultKKTBackend`](@ref) | +| `System` | See [KKT systems](@ref) | [`KKT.DefaultKKTSystem`](@ref) | + +## Linear Algebra + +These parameters control the linear algebra implementation + +| Parameter | Description | Default | +|:----------|:------------|:--------| +| `MatrixFactory` | See [`Factory`](@ref) | `Factory(SparseMatrixCSC)` + +## Presolve + +These parameters control the execution of the presolve phase. +They should be called as `"Presolve_"`. + ## Others | Parameter | Description | Default | diff --git a/src/IPM/HSD/HSD.jl b/src/IPM/HSD/HSD.jl index 1c073932..b50f1149 100644 --- a/src/IPM/HSD/HSD.jl +++ b/src/IPM/HSD/HSD.jl @@ -34,7 +34,7 @@ mutable struct HSD{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T} function HSD( dat::IPMData{T, Tv, Tb, Ta}, kkt_options::KKTOptions{T} ) where{T, Tv<:AbstractVector{T}, Tb<:AbstractVector{Bool}, Ta<:AbstractMatrix{T}} - + m, n = dat.nrow, dat.ncol p = sum(dat.lflag) + sum(dat.uflag) @@ -51,7 +51,7 @@ mutable struct HSD{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T} regD = tones(Tv, m) regG = one(T) - kkt = KKT.setup(kkt_options.Factory.T, dat.A; kkt_options.Factory.options...) + kkt = KKT.setup(dat.A, kkt_options.System, kkt_options.Backend) Tk = typeof(kkt) return new{T, Tv, Tb, Ta, Tk}(dat, @@ -161,7 +161,7 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher else hsd.dual_status = Sln_Unknown end - + # Check for optimal solution if ρp <= ϵp && ρd <= ϵd && ρg <= ϵg hsd.primal_status = Sln_Optimal @@ -169,7 +169,7 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher hsd.solver_status = Trm_Optimal return nothing end - + # Check for infeasibility certificates if max( norm(dat.A * pt.x, Inf), @@ -244,7 +244,7 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T} hsd.pt.y .= zero(T) hsd.pt.zl .= one(T) .* dat.lflag hsd.pt.zu .= one(T) .* dat.uflag - + hsd.pt.τ = one(T) hsd.pt.κ = one(T) @@ -268,12 +268,12 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T} if params.OutputLevel > 0 # Display log @printf "%4d" hsd.niter - + # Objectives ϵ = dat.objsense ? one(T) : -one(T) @printf " %+14.7e" ϵ * hsd.primal_objective @printf " %+14.7e" ϵ * hsd.dual_objective - + # Residuals @printf " %8.2e" max(hsd.res.rp_nrm, hsd.res.ru_nrm) @printf " %8.2e" hsd.res.rd_nrm @@ -306,14 +306,14 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T} || hsd.solver_status == Trm_DualInfeasible ) break - elseif hsd.niter >= params.IterationsLimit + elseif hsd.niter >= params.IterationsLimit hsd.solver_status = Trm_IterationLimit break elseif ttot >= params.TimeLimit hsd.solver_status = Trm_TimeLimit break end - + # TODO: step # For now, include the factorization in the step function @@ -325,7 +325,7 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T} if isa(err, PosDefException) || isa(err, SingularException) # Numerical trouble while computing the factorization hsd.solver_status = Trm_NumericalProblem - + elseif isa(err, OutOfMemoryError) # Out of memory hsd.solver_status = Trm_MemoryLimit @@ -349,4 +349,4 @@ function ipm_optimize!(hsd::HSD{T}, params::IPMOptions{T}) where{T} return nothing -end \ No newline at end of file +end diff --git a/src/IPM/MPC/MPC.jl b/src/IPM/MPC/MPC.jl index eec253e6..76b37782 100644 --- a/src/IPM/MPC/MPC.jl +++ b/src/IPM/MPC/MPC.jl @@ -76,7 +76,7 @@ mutable struct MPC{T, Tv, Tb, Ta, Tk} <: AbstractIPMOptimizer{T} regP = tones(Tv, n) regD = tones(Tv, m) - kkt = KKT.setup(kkt_options.Factory.T, dat.A; kkt_options.Factory.options...) + kkt = KKT.setup(dat.A, kkt_options.System, kkt_options.Backend) Tk = typeof(kkt) return new{T, Tv, Tb, Ta, Tk}(dat, diff --git a/src/KKT/Cholmod/cholmod.jl b/src/KKT/Cholmod/cholmod.jl new file mode 100644 index 00000000..255db19a --- /dev/null +++ b/src/KKT/Cholmod/cholmod.jl @@ -0,0 +1,70 @@ +module TlpCholmod + +using LinearAlgebra +using SparseArrays +using SuiteSparse.CHOLMOD + +using ..KKT: AbstractKKTBackend, AbstractKKTSolver +using ..KKT: AbstractKKTSystem, K1, K2 +import ..KKT: setup, update!, solve!, backend, linear_system + +""" + Backend + +CHOLMOD backend for solving linear systems. + +See [`CholmodSolver`](@ref) for further details. +""" +struct Backend <: AbstractKKTBackend end + +""" + CholmodSolver{T,S<:AbstractKKTSystem} + +CHOLMOD-based KKT solver. + +# Supported arithmetics +* `Float64` + +# Supported systems +* [`K2`](@ref) via ``LDLᵀ`` factorization +* [`K1`](@ref) via Cholesky (``LLᵀ``) factorization + +# Examples + +* To solve the augmented system with CHOLMOD's ``LDL^{T}`` factorization: +```julia +set_parameter(tlp_model, "KKT_Backend", Tulip.KKT.TlpCholmod.Backend()) +set_parameter(tlp_model, "KKT_System", Tulip.KKT.K2()) +``` + +* To solve the normal equations system with CHOLMOD's Cholesky factorization: +```julia +set_parameter(tlp_model, "KKT_Backend", Tulip.KKT.TlpCholmod.Backend()) +set_parameter(tlp_model, "KKT_System", Tulip.KKT.K1()) +``` +""" +mutable struct CholmodSolver{T,S} <: AbstractKKTSolver{T} + # Problem data + m::Int + n::Int + A::SparseMatrixCSC{T,Int} + + # Workspace + # TODO: store K as CHOLMOD.Sparse instead of SparseMatrixCSC + θ::Vector{T} # Diagonal scaling + regP::Vector{T} # Primal regularization + regD::Vector{T} # Dual regularization + K::SparseMatrixCSC{T,Int} # KKT matrix + F::CHOLMOD.Factor{T} # Factorization + ξ::Vector{T} # RHS of KKT system +end + +backend(::CholmodSolver) = "CHOLMOD" + +# Convert to sparse matrix if other type is used +setup(A, system, backend::Backend) = setup(convert(SparseMatrixCSC, A), system, backend) + +include("spd.jl") # Normal equations +include("sqd.jl") # Augmented system + +end # module diff --git a/src/KKT/Cholmod/spd.jl b/src/KKT/Cholmod/spd.jl new file mode 100644 index 00000000..cf479aa4 --- /dev/null +++ b/src/KKT/Cholmod/spd.jl @@ -0,0 +1,70 @@ +const CholmodSPD = CholmodSolver{Float64,K1} + +linear_system(::CholmodSPD) = "Normal equations (K1)" + +function setup(A::SparseMatrixCSC{Float64}, ::K1, ::Backend) + m, n = size(A) + + θ = ones(Float64, n) + regP = ones(Float64, n) + regD = ones(Float64, m) + ξ = zeros(Float64, m) + + # TODO: analyze + in-place A*D*A' product + K = sparse(A * A') + spdiagm(0 => regD) + + # TODO: PSD-ness checks + F = cholesky(Symmetric(K)) + + return CholmodSolver{Float64,K1}(m, n, A, θ, regP, regD, K, F, ξ) +end + +function update!(kkt::CholmodSPD, θ, regP, regD) + m, n = kkt.m, kkt.n + + # Sanity checks + length(θ) == n || throw(DimensionMismatch( + "length(θ)=$(length(θ)) but KKT solver has n=$n." + )) + length(regP) == n || throw(DimensionMismatch( + "length(regP)=$(length(regP)) but KKT solver has n=$n" + )) + length(regD) == m || throw(DimensionMismatch( + "length(regD)=$(length(regD)) but KKT solver has m=$m" + )) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + # Form normal equations matrix + # TODO: use in-place update of S + D = inv(Diagonal(kkt.θ .+ kkt.regP)) + kkt.K = (kkt.A * D * kkt.A') + spdiagm(0 => kkt.regD) + + # Update factorization + cholesky!(kkt.F, Symmetric(kkt.K), check=false) + issuccess(kkt.F) || throw(PosDefException(0)) + + return nothing +end + +function solve!(dx, dy, kkt::CholmodSPD, ξp, ξd) + m, n = kkt.m, kkt.n + + D = inv(Diagonal(kkt.θ .+ kkt.regP)) + copyto!(kkt.ξ, ξp) + mul!(kkt.ξ, kkt.A, D * ξd, true, true) + + # Solve normal equations + # CHOLMOD doesn't have in-place solve, so this line will allocate + dy .= kkt.F \ kkt.ξ + + # Recover dx + copyto!(dx, ξd) + mul!(dx, kkt.A', dy, 1.0, -1.0) + lmul!(D, dx) + + # TODO: iterative refinement + return nothing +end diff --git a/src/KKT/Cholmod/sqd.jl b/src/KKT/Cholmod/sqd.jl new file mode 100644 index 00000000..4b1e609e --- /dev/null +++ b/src/KKT/Cholmod/sqd.jl @@ -0,0 +1,74 @@ +const CholmodSQD = CholmodSolver{Float64,K2} + +linear_system(::CholmodSQD) = "Augmented system (K2)" + +function setup(A::SparseMatrixCSC{Float64,Int}, ::K2, ::Backend) + m, n = size(A) + + θ = ones(Float64, n) + regP = ones(Float64, n) + regD = ones(Float64, m) + ξ = zeros(Float64, m+n) + + K = [ + spdiagm(0 => -θ) A'; + spzeros(Float64, m, n) spdiagm(0 => ones(m)) + ] + + # TODO: Symbolic factorization only + F = ldlt(Symmetric(K)) + + return CholmodSolver{Float64,K2}(m, n, A, θ, regP, regD, K, F, ξ) +end + +function update!(kkt::CholmodSQD, θ, regP, regD) + m, n = kkt.m, kkt.n + + # Sanity checks + length(θ) == n || throw(DimensionMismatch( + "length(θ)=$(length(θ)) but KKT solver has n=$n." + )) + length(regP) == n || throw(DimensionMismatch( + "length(regP)=$(length(regP)) but KKT solver has n=$n" + )) + length(regD) == m || throw(DimensionMismatch( + "length(regD)=$(length(regD)) but KKT solver has m=$m" + )) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + # Update KKT matrix + # K is stored as upper-triangular, and only its diagonal is changed + @inbounds for j in 1:kkt.n + k = kkt.K.colptr[1+j] - 1 + kkt.K.nzval[k] = -kkt.θ[j] - regP[j] + end + @inbounds for i in 1:kkt.m + k = kkt.K.colptr[1+kkt.n+i] - 1 + kkt.K.nzval[k] = regD[i] + end + + ldlt!(kkt.F, Symmetric(kkt.K)) + return nothing +end + +function solve!(dx, dy, kkt::CholmodSQD, ξp, ξd) + m, n = kkt.m, kkt.n + + # Setup right-hand side + @views copyto!(kkt.ξ[1:n], ξd) + @views copyto!(kkt.ξ[(n+1):end], ξp) + + # Solve augmented system + # CHOLMOD doesn't have in-place solve, so this line will allocate + δ = kkt.F \ kkt.ξ + + # Recover dx, dy + @views copyto!(dx, δ[1:n]) + @views copyto!(dy, δ[(n+1):end]) + + # TODO: iterative refinement + return nothing +end diff --git a/src/KKT/Dense/lapack.jl b/src/KKT/Dense/lapack.jl new file mode 100644 index 00000000..fd581288 --- /dev/null +++ b/src/KKT/Dense/lapack.jl @@ -0,0 +1,121 @@ +module TlpDense + +using LinearAlgebra +using LinearAlgebra:BlasReal + +using ..KKT: AbstractKKTBackend, AbstractKKTSolver +using ..KKT: AbstractKKTSystem, K1, K2 +import ..KKT: setup, update!, solve!, backend, linear_system + +""" + Backend + +Dense linear algebra backend for solving linear systems. + +See [`DenseSolver`](@ref) for further details. +""" +struct Backend <: AbstractKKTBackend end + +""" + DenseSolver{T} + +Dense linear algebra-based KKT solver. + +# Supported arithmetics + +All arithmetics are supported. + +BLAS/LAPACK routines are used automatically with `Float32` and `Float64` arithmetic. + +# Supported systems +* [`K1`](@ref) via Cholesky factorization +""" +mutable struct DenseSolver{T} <: AbstractKKTSolver{T} + # Problem data + m::Int + n::Int + A::Matrix{T} + + # Workspace + _A::Matrix{T} # Place-holder for scaled copy of A + θ::Vector{T} # Diagonal scaling + regP::Vector{T} # Primal regularization + regD::Vector{T} # Dual regularization + K::Matrix{T} # KKT matrix + ξ::Vector{T} # RHS of KKT system +end + +backend(::DenseSolver) = "Julia.LinearAlgebra" +backend(::DenseSolver{<:BlasReal}) = "LAPACK $(LinearAlgebra.BLAS.vendor())" +linear_system(::DenseSolver) = "Normal equations (K1)" + +function setup(A::Matrix{T}, ::K1, ::Backend) where{T} + m, n = size(A) + + _A = Matrix{T}(undef, m, n) + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) + K = Matrix{T}(undef, m, m) + ξ = zeros(T, m) + + return DenseSolver{T}(m, n, A, _A, θ, regP, regD, K, ξ) +end + +function update!(kkt::DenseSolver, θ, regP, regD) + m, n = kkt.m, kkt.n + + # Sanity checks + length(θ) == n || throw(DimensionMismatch( + "length(θ)=$(length(θ)) but KKT solver has n=$n." + )) + length(regP) == n || throw(DimensionMismatch( + "length(regP)=$(length(regP)) but KKT solver has n=$n" + )) + length(regD) == m || throw(DimensionMismatch( + "length(regD)=$(length(regD)) but KKT solver has m=$m" + )) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + # 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 = sqrt(inv(Diagonal(kkt.θ .+ kkt.regP))) + rmul!(kkt._A, D) # B = A * √D + mul!(kkt.K, kkt._A, transpose(kkt._A), true, false) # Now K = A*D*A' + # Finally, add dual regularizations to the diagonal + @inbounds for i in 1:kkt.m + kkt.K[i, i] += kkt.regD[i] + end + + # In-place Cholesky factorization + cholesky!(Symmetric(kkt.K)) + + return nothing +end + +function solve!(dx, dy, kkt::DenseSolver{T}, ξp, ξd) where{T} + m, n = kkt.m, kkt.n + + # Set-up right-hand side + D = inv(Diagonal(kkt.θ .+ kkt.regP)) + copyto!(dy, ξp) + mul!(dy, kkt.A, D * ξd, true, true) + + # Solve normal equations + ldiv!(UpperTriangular(kkt.K)', dy) + ldiv!(UpperTriangular(kkt.K) , dy) + + # Recover dx + copyto!(dx, ξd) + mul!(dx, kkt.A', dy, one(T), -one(T)) + lmul!(D, dx) + + # TODO: Iterative refinement + return nothing +end + +end # module diff --git a/src/KKT/KKT.jl b/src/KKT/KKT.jl index f8ea8044..9196c813 100644 --- a/src/KKT/KKT.jl +++ b/src/KKT/KKT.jl @@ -1,20 +1,43 @@ module KKT using LinearAlgebra +using SparseArrays -import ..Tulip.Factory +using LinearOperators -export AbstractKKTSolver, KKTOptions +export AbstractKKTSystem, AbstractKKTBackend, AbstractKKTSolver +export KKTOptions + +""" + AbstractKKTSystem + +Abstract type for KKT systems +""" +abstract type AbstractKKTSystem end + +include("systems.jl") + +""" + AbstractKKTBackend + +Abstract type for KKT backend, i.e., the actual linear solver. +""" +abstract type AbstractKKTBackend end + +""" + DefaultKKTBackend + +Default setting for KKT backend. + +Currently defaults to [`TlpCholmod.Backend`](@ref) for `Float64` arithmetic, + and [`TlpLDLFact.Backend`](@ref) otherwise. +""" +struct DefaultKKTBackend <: AbstractKKTBackend end """ AbstractKKTSolver{T} -Abstract container for solving an augmented system -``` - [-(Θ⁻¹ + Rp) Aᵀ] [dx] = [ξd] - [ A Rd] [dy] [ξp] -``` -where `ξd` and `ξp` are given right-hand side. +Abstract container for solving KKT systems in arithmetic `T`. """ abstract type AbstractKKTSolver{T} end @@ -24,22 +47,20 @@ abstract type AbstractKKTSolver{T} end KKT solver options. """ Base.@kwdef mutable struct KKTOptions{T} - Factory::Factory{<:AbstractKKTSolver} = default_factory(T) + Backend::AbstractKKTBackend = DefaultKKTBackend() + System::AbstractKKTSystem = DefaultKKTSystem() end - """ - setup(T::Type{<:AbstractKKTSolver}, args...; kwargs...) + setup(A, system, backend; kwargs...) Instantiate a KKT solver object. """ -function setup(Ts::Type{<:AbstractKKTSolver}, args...; kwargs...) - return Ts(args...; kwargs...) -end +function setup end -# +# # Specialized implementations should extend the functions below -# +# """ update!(kkt, θinv, regP, regD) @@ -61,7 +82,6 @@ for given right-hand sides `ξd` and `ξp`. """ function update! end - """ solve!(dx, dy, kkt, ξp, ξd) @@ -84,7 +104,7 @@ function solve! end Return the arithmetic used by the solver. """ -arithmetic(kkt::AbstractKKTSolver{T}) where T = T +arithmetic(::AbstractKKTSolver{T}) where T = T """ backend(kkt) @@ -101,25 +121,23 @@ Return which system is solved by the kkt solver. linear_system(::AbstractKKTSolver) = "Unkown" # Generic tests -include("test.jl") +include("Test/test.jl") # Custom linear solvers -include("lapack.jl") -include("cholmod.jl") -include("ldlfact.jl") -include("krylov.jl") - -""" - default_factory(T) - -Use CHOLMOD for `Float64` and LDLFactorizations otherwise. -""" -function default_factory(::Type{T}) where{T} +include("Dense/lapack.jl") +include("Cholmod/cholmod.jl") +include("LDLFactorizations/ldlfact.jl") +const TlpLDLFact = TlpLDLFactorizations +include("Krylov/krylov.jl") + +# Default backend and system choices +function setup(A, ::DefaultKKTSystem, ::DefaultKKTBackend) + T = eltype(A) if T == Float64 - return Factory(CholmodSolver; normal_equations=false) + return setup(A, K2(), TlpCholmod.Backend()) else - return Factory(LDLFactSQD) + return setup(A, K2(), TlpLDLFact.Backend()) end end -end # module \ No newline at end of file +end # module diff --git a/src/KKT/Krylov/defs.jl b/src/KKT/Krylov/defs.jl new file mode 100644 index 00000000..d8e37345 --- /dev/null +++ b/src/KKT/Krylov/defs.jl @@ -0,0 +1,48 @@ +const _KRYLOV_SPD = Union{ + Krylov.CgSolver, + Krylov.CrSolver, +} + +const _KRYLOV_SID = Union{ + Krylov.MinresSolver, + Krylov.MinresQlpSolver, + Krylov.SymmlqSolver +} + +const _KRYLOV_SQD = Union{ + Krylov.TricgSolver, + Krylov.TrimrSolver, +} + +const _KRYLOV_LN = Union{ + Krylov.LnlqSolver, + Krylov.CraigSolver, + Krylov.CraigmrSolver, +} + +const _KRYLOV_LS = Union{ + Krylov.LslqSolver, + Krylov.LsqrSolver, + Krylov.LsmrSolver, +} + +# Helper functions +for (KS, fun) in [ + (Krylov.CgSolver,Krylov.cg!) + (Krylov.CrSolver,Krylov.cr!) + (Krylov.MinresSolver,Krylov.minres!) + (Krylov.MinresQlpSolver,Krylov.minres_qlp!) + (Krylov.SymmlqSolver,Krylov.symmlq!) + (Krylov.TricgSolver,Krylov.tricg!) + (Krylov.TrimrSolver,Krylov.trimr!) + (Krylov.LnlqSolver,Krylov.lnlq!) + (Krylov.CraigSolver,Krylov.craig!) + (Krylov.CraigmrSolver,Krylov.craigmr!) + (Krylov.LslqSolver,Krylov.lslq!) + (Krylov.LsqrSolver,Krylov.lsqr!) + (Krylov.LsmrSolver,Krylov.lsmr!) +] + @eval begin + @inline _krylov!(solver::$KS, args...; kwargs...) = $(fun)(solver, args...; kwargs...) + end +end diff --git a/src/KKT/Krylov/krylov.jl b/src/KKT/Krylov/krylov.jl new file mode 100644 index 00000000..7353bb31 --- /dev/null +++ b/src/KKT/Krylov/krylov.jl @@ -0,0 +1,57 @@ +module TlpKrylov + +using LinearAlgebra + +using Krylov +using LinearOperators +const LO = LinearOperators + +using ..KKT: AbstractKKTBackend, AbstractKKTSolver +using ..KKT: AbstractKKTSystem, K1, K2 +import ..KKT: arithmetic, backend, linear_system +import ..KKT: setup, update!, solve! + +include("defs.jl") + +""" + Backend{KS<:Krylov.KrylovSolver,V<:AbstractVector} + +[Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl)-based backend for solving linear systems. + +The type is parametrized by: +* `KS<:Krylov.KrylovSolver`: workspace type for the Krylov method. + Also defines the Krylov method to be used. +* `V<:AbstractVector`: the vector storage type used within the Krylov method. + This should be set to `Vector{T}` (for arithmetic `T`) unless, e.g., one uses a GPU. + +See the [Krylov.jl documentation](https://juliasmoothoptimizers.github.io/Krylov.jl/dev/inplace/) for further details. + +# Example usage + +All the following examples assume everything runs on a CPU in `Float64` arithmetic. +* To use the conjugate gradient: +```julia +backend = KKT.TlpKrylov.Backend(Krylov.CgSolver, Vector{Float64}) +``` +* To use MINRES: +```julia +backend = KKT.TlpKrylov.Backend(Krylov.MinresSolver, Vector{Float64}) +``` +""" +struct Backend{KS,V} <: AbstractKKTBackend + krylov_solver::Type{KS} + vector_storage::Type{V} +end + +""" + AbstractKrylovSolver{T} + +Abstract type for Kyrlov-based linear solvers. +""" +abstract type AbstractKrylovSolver{T} <: AbstractKKTSolver{T} end + +include("spd.jl") +include("sid.jl") +include("sqd.jl") + +end # module diff --git a/src/KKT/Krylov/sid.jl b/src/KKT/Krylov/sid.jl new file mode 100644 index 00000000..fae6d49c --- /dev/null +++ b/src/KKT/Krylov/sid.jl @@ -0,0 +1,108 @@ +""" + SIDSolver +""" +mutable struct SIDSolver{T,V,Ta,KL,KS} <: AbstractKrylovSolver{T} + # Problem data + m::Int + n::Int + A::Ta + + # IPM-related workspace + θ::Vector{T} + regP::Vector{T} + regD::Vector{T} + + # Krylov-related workspace + Θp::Diagonal{T,V} + Θd::Diagonal{T,V} + ξ::V + opK::KL + + # Krylov solver & related options + atol::T + rtol::T + krylov_solver::KS + + # TODO: preconditioner +end + +backend(kkt::SIDSolver) = "$(typeof(kkt.krylov_solver))" +linear_system(kkt::SIDSolver) = "K2" + +function setup(A, ::K2, backend::Backend{KS,V}) where{KS<:_KRYLOV_SID,V} + Ta = typeof(A) + T = eltype(A) + T == eltype(V) || error("eltype(A)=$T incompatible with eltype of Krylov vector storage $V.") + + m, n = size(A) + + # Workspace + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) + + Θp = Diagonal(V(undef, n)) + Θd = Diagonal(V(undef, m)) + ξ = V(undef, m+n) + + # Define linear operator for the augmented system + # This linear operator is symmetric indefinite + opK = LO.LinearOperator(T, m+n, m+n, true, false, + (u, w, α, β) -> begin + @views u1 = u[1:n] + @views u2 = u[(n+1):(m+n)] + @views w1 = w[1:n] + @views w2 = w[n+1:n+m] + + mul!(u1, Θp, w1, α, β) + mul!(u1, A', w2, α, one(T)) + + mul!(u2, A , w1, α, β) + mul!(u2, Θd, w2, α, one(T)) + u + end + ) + + # Allocate Krylov solver's workspace + atol = sqrt(eps(T)) + rtol = sqrt(eps(T)) + krylov_solver = KS(m+n, m+n, V) + + return SIDSolver{T,V,Ta,typeof(opK),typeof(krylov_solver)}( + m, n, A, + θ, regP, regD, + Θp, Θd, ξ, + opK, + atol, rtol, + krylov_solver + ) +end + +function update!(kkt::SIDSolver, θ, regP, regD) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + copyto!(kkt.Θd.diag, regD) + copyto!(kkt.Θp.diag, -(kkt.θ .+ kkt.regP)) + + return nothing +end + +function solve!(dx, dy, kkt::SIDSolver{T}, ξp, ξd) where{T} + m, n = kkt.m, kkt.n + + @views copyto!(kkt.ξ[1:m], ξp) + @views copyto!(kkt.ξ[(m+1):(m+n)], ξd) + + # Solve the augmented system + _krylov!(kkt.krylov_solver, kkt.opK, kkt.ξ; atol=kkt.atol, rtol=kkt.rtol) + + # Recover dx, dy + copyto!(dx, kkt.krylov_solver.x[1:n]) + copyto!(dy, kkt.krylov_solver.x[(n+1):(m+n)]) + + # TODO: iterative refinement (?) + return nothing +end diff --git a/src/KKT/Krylov/spd.jl b/src/KKT/Krylov/spd.jl new file mode 100644 index 00000000..b3d5607d --- /dev/null +++ b/src/KKT/Krylov/spd.jl @@ -0,0 +1,110 @@ +""" + SPDSolver +""" +mutable struct SPDSolver{T,V,Ta,KL,KS} <: AbstractKrylovSolver{T} + # Problem data + m::Int + n::Int + A::Ta + + # IPM-related workspace + θ::Vector{T} + regP::Vector{T} + regD::Vector{T} + + # Krylov-related workspace + D::Diagonal{T,V} + Rd::Diagonal{T,V} + ξ::V + opK::KL + + # Krylov solver & related options + atol::T + rtol::T + krylov_solver::KS + + # TODO: preconditioner +end + +backend(kkt::SPDSolver) = "$(typeof(kkt.krylov_solver))" +linear_system(kkt::SPDSolver) = "K1" + +function setup(A, ::K1, backend::Backend{KS,V}) where{KS<:Union{_KRYLOV_SPD,_KRYLOV_SID},V} + Ta = typeof(A) + T = eltype(A) + T == eltype(V) || error("eltype(A)=$T incompatible with eltype of Krylov vector storage $V.") + + m, n = size(A) + + # Workspace + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) + + D = Diagonal(V(undef, n)) + Rd = Diagonal(V(undef, m)) + ξ = V(undef, m) + + # Define linear operator for normal equations system + # We need to allocate one temporary vector + # This linear operator is symmetric definite positive, + # so we only need to define + # u ⟵ α (A×D×A'+Rd) × w + β u + # i.e., u = α(ADA')×w + (αRd)×w + βu + vtmp = V(undef, n) + opK = LO.LinearOperator(T, m, m, true, true, + (u, w, α, β) -> begin + mul!(vtmp, A', w) + lmul!(D, vtmp) + mul!(u, A, vtmp, α, β) + mul!(u, Rd, w, α, one(T)) + u + end + ) + + # Allocate Krylov solver's workspace + atol = sqrt(eps(T)) + rtol = sqrt(eps(T)) + krylov_solver = KS(m, m, V) + + return SPDSolver{T,V,Ta,typeof(opK),typeof(krylov_solver)}( + m, n, A, + θ, regP, regD, + D, Rd, ξ, + opK, + atol, rtol, + krylov_solver + ) +end + +function update!(kkt::SPDSolver, θ, regP, regD) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + copyto!(kkt.Rd.diag, regD) + copyto!(kkt.D.diag, inv.(kkt.θ .+ kkt.regP)) + + return nothing +end + +function solve!(dx, dy, kkt::SPDSolver{T}, ξp, ξd) where{T} + m, n = kkt.m, kkt.n + + copyto!(kkt.ξ, ξp) + + mul!(kkt.ξ, kkt.A, kkt.D * ξd, true, true) + + # Solve the normal equations + _krylov!(kkt.krylov_solver, kkt.opK, kkt.ξ; atol=kkt.atol, rtol=kkt.rtol) + copyto!(dy, kkt.krylov_solver.x) + + # Recover dx + copyto!(dx, ξd) + mul!(dx, kkt.A', dy, one(T), -one(T)) + lmul!(kkt.D, dx) + + # TODO: iterative refinement (?) + return nothing +end diff --git a/src/KKT/Krylov/sqd.jl b/src/KKT/Krylov/sqd.jl new file mode 100644 index 00000000..1a2b2c48 --- /dev/null +++ b/src/KKT/Krylov/sqd.jl @@ -0,0 +1,100 @@ +""" + SQDSolver +""" +mutable struct SQDSolver{T,V,Ta,KS} <: AbstractKrylovSolver{T} + # Problem data + m::Int + n::Int + A::Ta + + # IPM-related workspace + θ::Vector{T} + regP::Vector{T} + regD::Vector{T} + + # Krylov-related workspace + Θp::Diagonal{T,V} + Θp⁻¹::Diagonal{T,V} + Θd::Diagonal{T,V} + Θd⁻¹::Diagonal{T,V} + ξp::V + ξd::V + + # Krylov solver & related options + atol::T + rtol::T + krylov_solver::KS + + # TODO: preconditioner +end + +backend(kkt::SQDSolver) = "$(typeof(kkt.krylov_solver))" +linear_system(kkt::SQDSolver) = "K2" + +function setup(A, ::K2, backend::Backend{KS,V}) where{KS<:_KRYLOV_SQD,V} + Ta = typeof(A) + T = eltype(A) + T == eltype(V) || error("eltype(A)=$T incompatible with eltype of Krylov vector storage $V.") + + m, n = size(A) + + # Workspace + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) + + Θp = Diagonal(V(undef, n)) + Θp⁻¹ = Diagonal(V(undef, n)) + Θd = Diagonal(V(undef, m)) + Θd⁻¹ = Diagonal(V(undef, m)) + ξp = V(undef, m) + ξd = V(undef, n) + + # Allocate Krylov solver's workspace + atol = sqrt(eps(T)) + rtol = sqrt(eps(T)) + krylov_solver = KS(m, n, V) + + return SQDSolver{T,V,Ta,typeof(krylov_solver)}( + m, n, A, + θ, regP, regD, + Θp, Θp⁻¹, Θd, Θd⁻¹, + ξp, ξd, + atol, rtol, + krylov_solver + ) +end + +function update!(kkt::SQDSolver, θ, regP, regD) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + copyto!(kkt.Θp.diag, -(kkt.θ .+ kkt.regP)) + copyto!(kkt.Θp⁻¹.diag, inv.(kkt.θ .+ kkt.regP)) # Θp⁻¹ will be negated by tricg/trimr + copyto!(kkt.Θd.diag, kkt.regD) + copyto!(kkt.Θd⁻¹.diag, inv.(kkt.regD)) + + return nothing +end + +function solve!(dx, dy, kkt::SQDSolver{T}, ξp, ξd) where{T} + copyto!(kkt.ξp, ξp) + copyto!(kkt.ξd, ξd) + + # Solve the augmented system + _krylov!(kkt.krylov_solver, kkt.A, kkt.ξp, kkt.ξd; + M=kkt.Θd⁻¹, + N=kkt.Θp⁻¹, + atol=kkt.atol, + rtol=kkt.rtol + ) + + # Recover dx, dy + copyto!(dx, kkt.krylov_solver.yₖ) + copyto!(dy, kkt.krylov_solver.xₖ) + + # TODO: iterative refinement (?) + return nothing +end diff --git a/src/KKT/LDLFactorizations/ldlfact.jl b/src/KKT/LDLFactorizations/ldlfact.jl new file mode 100644 index 00000000..dfd88b14 --- /dev/null +++ b/src/KKT/LDLFactorizations/ldlfact.jl @@ -0,0 +1,141 @@ +module TlpLDLFactorizations + +using LinearAlgebra +using SparseArrays + +using LDLFactorizations +const LDLF = LDLFactorizations + +using ..KKT: AbstractKKTBackend, AbstractKKTSolver +using ..KKT: AbstractKKTSystem, K1, K2 +import ..KKT: setup, update!, solve!, backend, linear_system + + +""" + Backend + +LDLFactorizations backend for solving linear systems. + +See [`LDLFactSolver`](@ref) for further details. +""" +struct Backend <: AbstractKKTBackend end + +""" + LDLFactSolver{T,S<:AbstractKKTSystem} + +[`LDLFactorizations.jl`](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl)-based KKT solver. + +# Supported arithmetics +* All arithmetics are supported + +# Supported systems +* [`K2`](@ref) via ``LDLᵀ`` factorization + +# Examples + +To solve the augmented system with LDLFactorizations' ``LDL^{T}`` factorization: +```julia +set_parameter(tlp_model, "KKT_Backend", Tulip.KKT.TlpLDLFact.Backend()) +set_parameter(tlp_model, "KKT_System", Tulip.KKT.K2()) +``` +""" +mutable struct LDLFactSolver{T,S} <: AbstractKKTSolver{T} + # Problem data + m::Int + n::Int + A::SparseMatrixCSC{T,Int} + + # Workspace + θ::Vector{T} # Diagonal scaling + regP::Vector{T} # Primal regularization + regD::Vector{T} # Dual regularization + K::SparseMatrixCSC{T,Int} # KKT matrix + F::LDLF.LDLFactorization{T} # Factorization + ξ::Vector{T} # RHS of KKT system +end + +backend(::LDLFactSolver) = "LDLFactorizations" +linear_system(::LDLFactSolver) = "Augmented system (K2)" + +# Convert A to sparse matrix if needed +setup(A, system, backend::Backend) = setup(convert(SparseMatrixCSC, A), system, backend) + +function setup(A::SparseMatrixCSC{T,Int}, ::K2, ::Backend) where{T} + m, n = size(A) + + θ = ones(T, n) + regP = ones(T, n) + regD = ones(T, m) + ξ = zeros(T, m+n) + + K = [ + spdiagm(0 => -θ) A'; + spzeros(T, m, n) spdiagm(0 => ones(T, m)) + ] + + # TODO: Symbolic factorization only + F = LDLF.ldl_analyze(Symmetric(K)) + + return LDLFactSolver{T,K2}(m, n, A, θ, regP, regD, K, F, ξ) +end + +function update!(kkt::LDLFactSolver{T,K2}, θ, regP, regD) where{T} + m, n = kkt.m, kkt.n + + # Sanity checks + length(θ) == n || throw(DimensionMismatch( + "length(θ)=$(length(θ)) but KKT solver has n=$n." + )) + length(regP) == n || throw(DimensionMismatch( + "length(regP)=$(length(regP)) but KKT solver has n=$n" + )) + length(regD) == m || throw(DimensionMismatch( + "length(regD)=$(length(regD)) but KKT solver has m=$m" + )) + + copyto!(kkt.θ, θ) + copyto!(kkt.regP, regP) + copyto!(kkt.regD, regD) + + # Update KKT matrix + # K is stored as upper-triangular, and only its diagonal is changed + @inbounds for j in 1:kkt.n + k = kkt.K.colptr[1+j] - 1 + kkt.K.nzval[k] = -kkt.θ[j] - regP[j] + end + @inbounds for i in 1:kkt.m + k = kkt.K.colptr[1+kkt.n+i] - 1 + kkt.K.nzval[k] = regD[i] + end + + # Update factorization + try + LDLF.ldl_factorize!(Symmetric(kkt.K), kkt.F) + catch err + isa(err, LDLF.SQDException) && throw(PosDefException(-1)) + rethrow(err) + end + + return nothing +end + +function solve!(dx, dy, kkt::LDLFactSolver{T,K2}, ξp, ξd) where{T} + m, n = kkt.m, kkt.n + + # Setup right-hand side + @views copyto!(kkt.ξ[1:n], ξd) + @views copyto!(kkt.ξ[(n+1):end], ξp) + + # Solve augmented system + # CHOLMOD doesn't have in-place solve, so this line will allocate + LDLF.ldiv!(kkt.F, kkt.ξ) + + # Recover dx, dy + @views copyto!(dx, kkt.ξ[1:n]) + @views copyto!(dy, kkt.ξ[(n+1):end]) + + # TODO: iterative refinement + return nothing +end + +end # module diff --git a/src/KKT/test.jl b/src/KKT/Test/test.jl similarity index 100% rename from src/KKT/test.jl rename to src/KKT/Test/test.jl diff --git a/src/KKT/cholmod.jl b/src/KKT/cholmod.jl deleted file mode 100644 index 6cdf7f7a..00000000 --- a/src/KKT/cholmod.jl +++ /dev/null @@ -1,323 +0,0 @@ -using SparseArrays -using SuiteSparse.CHOLMOD - -""" - CholmodSolver - -Uses CHOLMOD's factorization routines to solve the augmented system. - -To use an LDLᵀ factorization of the augmented system -(see [`CholmodSQD`](@ref)) -```julia -model.params.KKT.Factory = Tulip.Factory(CholmodSolver, normal_equations=false) -``` - -To use a Cholesky factorization of the normal equations system -(see [`CholmodSPD`](@ref)) -```julia -model.params.KKT.Factory = Tulip.Factory(CholmodSolver, normal_equations=true) -``` - -!!! warning - CHOLMOD can only be used with `Float64` arithmetic. -""" -abstract type CholmodSolver <: AbstractKKTSolver{Float64} end - -function CholmodSolver( - A::AbstractMatrix{Float64}; - normal_equations::Bool=false, - kwargs... -) - if normal_equations - return CholmodSPD(A; kwargs...) - else - return CholmodSQD(A; kwargs...) - end -end - - -# ============================================================================== -# CholmodSQD -# ============================================================================== - -""" - CholmodSQD - -Linear solver for the 2x2 augmented system with ``A`` sparse. - -Uses an LDLᵀ factorization of the quasi-definite augmented system. -""" -mutable struct CholmodSQD <: CholmodSolver - m::Int # Number of rows - n::Int # Number of columns - - # TODO: allow for user-provided ordering, - # and add flag (default to false) to know whether user ordering should be used - - # Problem data - A::SparseMatrixCSC{Float64} - θ::Vector{Float64} # diagonal scaling - regP::Vector{Float64} # primal-dual regularization - regD::Vector{Float64} # dual regularization - ξ::Vector{Float64} # right-hand side of the augmented system - - # Left-hand side matrix - S::SparseMatrixCSC{Float64, Int} # TODO: use `CHOLMOD.Sparse` instead - - # Factorization - F::CHOLMOD.Factor{Float64} - - # TODO: constructor with initial memory allocation - function CholmodSQD(A::AbstractMatrix{Float64}) - m, n = size(A) - θ = ones(Float64, n) - regP = ones(Float64, n) - regD = ones(Float64, m) - ξ = zeros(Float64, n+m) - - S = [ - spdiagm(0 => -θ) A'; - spzeros(Float64, m, n) spdiagm(0 => ones(m)) - ] - - # TODO: Symbolic factorization only - F = ldlt(Symmetric(S)) - - return new(m, n, A, θ, regP, regD, ξ, S, F) - end - -end - -backend(::CholmodSQD) = "CHOLMOD" -linear_system(::CholmodSQD) = "Augmented system" - -""" - update!(kkt, θ, regP, regD) - -Update LDLᵀ factorization of the augmented system. - -Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute - the factorization. -Throws a `PosDefException` if matrix is not quasi-definite. -""" -function update!( - kkt::CholmodSQD, - θ::AbstractVector{Float64}, - regP::AbstractVector{Float64}, - regD::AbstractVector{Float64} -) - # Sanity checks - length(θ) == kkt.n || throw(DimensionMismatch( - "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - # Update diagonal scaling - kkt.θ .= θ - # Update regularizers - kkt.regP .= regP - kkt.regD .= regD - - # Update S. S is stored as upper-triangular and only its diagonal changes. - @inbounds for j in 1:kkt.n - k = kkt.S.colptr[1+j] - 1 - kkt.S.nzval[k] = -kkt.θ[j] - regP[j] - end - @inbounds for i in 1:kkt.m - k = kkt.S.colptr[1+kkt.n+i] - 1 - kkt.S.nzval[k] = regD[i] - end - - # `S` is first converted into CHOLMOD's internal data structure, - # so this line allocates. - ldlt!(kkt.F, Symmetric(kkt.S)) - - return nothing -end - -""" - solve!(dx, dy, kkt, ξp, ξd) - -Solve the augmented system, overwriting `dx, dy` with the result. -""" -function solve!( - dx::Vector{Float64}, dy::Vector{Float64}, - kkt::CholmodSQD, - ξp::Vector{Float64}, ξd::Vector{Float64} -) - m, n = kkt.m, kkt.n - - # Set-up the right-hand side - @inbounds for i in 1:n - kkt.ξ[i] = ξd[i] - end - @inbounds for i in 1:m - kkt.ξ[i+n] = ξp[i] - end - - # Solve augmented system - d = kkt.F \ kkt.ξ - - # Recover dx, dy - @inbounds for i in 1:n - dx[i] = d[i] - end - @inbounds for i in 1:m - dy[i] = d[i+n] - end - - # Iterative refinement - # TODO: - # * Max number of refine steps - # * Check for residuals before refining - # * Check whether residuals did improve - # resP = kkt.A*dx + kkt.regD .* dy - ξp - # resD = - dx .* (kkt.θ + kkt.regP) + kkt.A' * dy - ξd - # println("\n|resP| = $(norm(resP, Inf))\n|resD| = $(norm(resD, Inf))") - - # ξ1 = [resD; resP] - # d1 = kkt.F \ ξ1 - - # # Update search direction - # @views dx .-= d1[1:n] - # @views dy .-= d1[(n+1):(m+n)] - - return nothing -end - -# ============================================================================== -# CholmodSPD -# ============================================================================== - -""" - CholmodSPD - -Linear solver for the 2x2 augmented system -``` - [-(Θ⁻¹ + Rp) Aᵀ] [dx] = [ξd] - [ A Rd] [dy] [ξp] -``` -with ``A`` sparse. - -Uses a Cholesky factorization of the positive definite normal equations system -``` -(A * ((Θ⁻¹ + Rp)⁻¹ * Aᵀ + Rd) dy = ξp + A * (θ⁻¹ + Rp)⁻¹ * ξd - dx = (Θ⁻¹ + Rp)⁻¹ * (Aᵀ * dy - ξd) -``` -""" -mutable struct CholmodSPD <: CholmodSolver - m::Int # Number of rows - n::Int # Number of columns - - # TODO: allow for user-provided ordering, - # and add flag (default to false) to know whether user ordering should be used - - # Problem data - A::SparseMatrixCSC{Float64, Int} - θ::Vector{Float64} # Diagonal scaling - # Regularization - regP::Vector{Float64} # primal - regD::Vector{Float64} # dual - - # Factorization - F::CHOLMOD.Factor{Float64} - - # Constructor and initial memory allocation - # TODO: symbolic only + allocation - function CholmodSPD(A::AbstractMatrix{Float64}) - m, n = size(A) - θ = ones(Float64, n) - - S = sparse(A * A') + spdiagm(0 => ones(m)) - - # TODO: PSD-ness checks - F = cholesky(Symmetric(S)) - - return new(m, n, A, θ, zeros(Float64, n), ones(Float64, m), F) - end - -end - -backend(::CholmodSPD) = "CHOLMOD - Cholesky" -linear_system(::CholmodSPD) = "Normal equations" - -""" - update!(kkt, θ, regP, regD) - -Compute normal equation system matrix, and update the factorization. -""" -function update!( - kkt::CholmodSPD, - θ::AbstractVector{Float64}, - regP::AbstractVector{Float64}, - regD::AbstractVector{Float64} -) - - # Sanity checks - length(θ) == kkt.n || throw(DimensionMismatch( - "θ has length $(length(θ)) but linear solver is for n=$(kkt.n)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - kkt.θ .= θ - kkt.regP .= regP # Primal regularization is disabled for normal equations - kkt.regD .= regD - - # Re-compute factorization - # D = (Θ^{-1} + Rp)^{-1} - D = Diagonal(one(Float64) ./ (kkt.θ .+ kkt.regP)) - Rd = spdiagm(0 => kkt.regD) - S = kkt.A * D * kkt.A' + Rd - - # Update factorization - cholesky!(kkt.F, Symmetric(S), check=false) - issuccess(kkt.F) || throw(PosDefException(0)) - - return nothing -end - -""" - solve!(dx, dy, kkt, ξp, ξd) - -Solve the augmented system, overwriting `dx, dy` with the result. -""" -function solve!( - dx::Vector{Float64}, dy::Vector{Float64}, - kkt::CholmodSPD, - ξp::Vector{Float64}, ξd::Vector{Float64} -) - m, n = kkt.m, kkt.n - - d = one(Float64) ./ (kkt.θ .+ kkt.regP) - D = Diagonal(d) - - # Set-up right-hand side - ξ_ = ξp .+ kkt.A * (D * ξd) - - # Solve augmented system - dy .= (kkt.F \ ξ_) - - # Recover dx - dx .= D * (kkt.A' * dy - ξd) - - # TODO: Iterative refinement - # * Max number of refine steps - # * Check for residuals before refining - # * Check whether residuals did improve - # resP = kkt.A * dx + kkt.regD .* dy - ξp - # resD = - D \ dx + kkt.A' * dy - ξd - # println("\n|resP| = $(norm(resP, Inf))\n|resD| = $(norm(resD, Inf))") - - - return nothing -end diff --git a/src/KKT/krylov.jl b/src/KKT/krylov.jl deleted file mode 100644 index 62019c0f..00000000 --- a/src/KKT/krylov.jl +++ /dev/null @@ -1,423 +0,0 @@ -import LinearOperators -const LO = LinearOperators - -""" - AbstractKrylovSolver{T} - -Abstract type for Kyrlov-based linear solvers. -""" -abstract type AbstractKrylovSolver{T} <: AbstractKKTSolver{T} end - -# ============================================================================== -# KrylovSPD: -# ============================================================================== - -""" - KrylovSPD{T, F, Tv, Ta} - -Krylov-based **S**ymmetric **P**ositive-**D**efinite (SPD) linear solver. - -Applies a Krylov method to the normal equations systems, then recovers a solution -to the augmented system. -The selected Krylov method must therefore handle symmetric positive-definite systems. -Suitable methods are CG or CR. - -A `KrylovSPD` is selected as follows -```julia -model.params.KKT.Factory = Tulip.Factory( - KrylovSPD, method=Krylov.cg, - atol=1e-12, rtol=1e-12 -) -``` - -The `method` argument is a function `f::F` whose signature must match -```julia -dy, _ = f(S, ξ; atol, rtol) -``` -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 KrylovSPD{T, F, Tv, Ta} <: AbstractKrylovSolver{T} - m::Int - n::Int - - f::F # Krylov function - atol::T # absolute tolerance - rtol::T # relative tolerance - - # Problem data - A::Ta # Constraint matrix - θ::Tv # scaling - regP::Tv # primal regularization - regD::Tv # dual regularization - - function KrylovSPD(f::Function, A::Ta; - atol::T=eps(T)^(3 // 4), - rtol::T=eps(T)^(3 // 4) - ) where{T, Ta <: AbstractMatrix{T}} - F = typeof(f) - m, n = size(A) - - return new{T, F, Vector{T}, Ta}( - m, n, - f, atol, rtol, - A, ones(T, n), ones(T, n), ones(T, m) - ) - end -end - -setup(::Type{KrylovSPD}, A; method, kwargs...) = KrylovSPD(method, A; kwargs...) - -backend(kkt::KrylovSPD) = "Krylov ($(kkt.f))" -linear_system(::KrylovSPD) = "Normal equations" - - -""" - update!(kkt::KrylovSPD, θ, regP, regD) - -Update diagonal scaling ``θ`` and primal-dual regularizations. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - # Update diagonal scaling - kkt.θ .= θ - # Update regularizers - kkt.regP .= regP - kkt.regD .= regD - - return nothing -end - -""" - solve!(dx, dy, kkt::KrylovSPD, ξp, ξd) - -Solve the normal equations system using the selected Krylov method. -""" -function solve!( - 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(T) ./ (kkt.θ .+ kkt.regP) - D = Diagonal(d) - - # Set-up right-hand side - ξ_ = ξp .+ kkt.A * (D * ξd) - - # Form linear operator S = A * D * A' + Rd - # Here we pre-allocate the intermediary vector - v1 = zeros(T, n) - opS = LO.LinearOperator(T, m, m, true, true, - (v2, w, α, β) -> begin - mul!(v1, kkt.A', w) - v1 .*= d - mul!(v2, kkt.A, v1, α, β) - v2 .+= α * kkt.regD .* w - v2 - end - ) - - # Solve normal equations - _dy, stats = kkt.f(opS, ξ_, atol=kkt.atol, rtol=kkt.rtol) - dy .= _dy - - # Recover dx - dx .= D * (kkt.A' * dy - ξd) - - return nothing -end - - -# ============================================================================== -# KrylovSID: -# ============================================================================== - -""" - KrylovSID{T, F, Tv, Ta} - -Krylov-based **S**ymmetric **I**n**D**efinite (SID) linear solver. - -Applies a Krylov method to solve the augmented system, -without exploiting its quasi-definiteness properties. -The selected Krylov method must therefore handle symmetric indefinite systems. -Suitable methods are MINRES or MINRES-QLP. - -A `KrylovSID` is selected as follows -```julia -model.params.KKT.Factory = Tulip.Factory( - KrylovSID, method=Krylov.minres, - atol=1e-12, rtol=1e-12 -) -``` - -The `method` argument is a function `f::F` whose signature must match -```julia -Δ, _ = f(K, ξ; atol, rtol) -``` -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 KrylovSID{T, F, Tv, Ta} <: AbstractKrylovSolver{T} - m::Int - n::Int - - f::F # Krylov function - atol::T # absolute tolerance - rtol::T # relative tolerance - - # Problem data - A::Ta # Constraint matrix - θ::Tv # scaling - regP::Tv # primal regularization - regD::Tv # dual regularization - - function KrylovSID(f::Function, A::Ta; - atol::T=eps(T)^(3 // 4), - rtol::T=eps(T)^(3 // 4) - ) where{T, Ta <: AbstractMatrix{T}} - F = typeof(f) - m, n = size(A) - - return new{T, F, Vector{T}, Ta}( - m, n, - f, atol, rtol, - A, ones(T, n), ones(T, n), ones(T, m) - ) - end -end - -setup(::Type{KrylovSID}, A; method, kwargs...) = KrylovSID(method, A; kwargs...) - -backend(kkt::KrylovSID) = "Krylov ($(kkt.f))" -linear_system(::KrylovSID) = "Augmented system" - - -""" - update!(kkt::KrylovSID, θ, regP, regD) - -Update diagonal scaling ``θ`` and primal-dual regularizations. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - # Update diagonal scaling - kkt.θ .= θ - # Update regularizers - kkt.regP .= regP - kkt.regD .= regD - - return nothing -end - -""" - solve!(dx, dy, kkt::KrylovSID, ξp, ξd) - -Solve the augmented system using the selected Krylov method. -""" -function solve!( - dx::Vector{T}, dy::Vector{T}, - kkt::KrylovSID{T}, - ξp::Vector{T}, ξd::Vector{T} -) where{T} - m, n = kkt.m, kkt.n - - # Setup - A = kkt.A - D = Diagonal(-(kkt.θ .+ kkt.regP)) # D = Θ⁻¹ + Rp - - # Set-up right-hand side - ξ = [ξd; ξp] - - # Form linear operator K = [-(Θ⁻¹ + Rp) Aᵀ; A Rd] - opK = LO.LinearOperator(T, n+m, n+m, true, false, - (z, w, α, β) -> begin - @views z1 = z[1:n] - @views z2 = z[n+1:n+m] - - @views w1 = w[1:n] - @views w2 = w[n+1:n+m] - - # z1 = -α(Θ⁻¹ + Rp) w1 + αA'w2 + βz1 - mul!(z1, D, w1, α, β) # z1 = -α(Θ⁻¹ + Rp) w1 + β z1 - mul!(z1, A', w2, α, one(T)) # z1 += α A' w2 - - # z2 = A w1 + Rd w2 - mul!(z2, A, w1, α, β) # z2 = α A w1 + β z2 - mul!(z2, Diagonal(kkt.regD), w2, α, one(T)) # z2 += α * Rd * w2 - z - end - ) - - # Solve augmented system - Δ, stats = kkt.f(opK, ξ, atol=kkt.atol, rtol=kkt.rtol) - - # Recover dx, dy - dx .= Δ[1:n] - dy .= Δ[n+1:n+m] - - return nothing -end - - -# ============================================================================== -# KrylovSQD: -# ============================================================================== - -""" - KrylovSQD{T, F, Tv, Ta} - -Krylov-based **S**ymmetric **Q**uasi-**D**efinite (SQD) linear solver. - -Applies a Krylov method to solve the augmented system, taking advantage of -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 `KrylovSID` is selected as follows -```julia -model.params.KKT.Factory = Tulip.Factory( - KrylovSQD, method=Krylov.trimr, - atol=1e-12, rtol=1e-12 -) -``` - -The `method` argument is a function `f::F` whose signature must match -``` -dy, dx, _ = f(A, ξp, ξd; M, N, atol, rtol) -``` -where the augmented system is of the form -``` -[ -N⁻¹ Aᵀ ] [dx] = [ξd] -[ A M⁻¹] [dy] = [ξp] -``` -i.e., ``N = (Θ^{-1} + R_{p})^{-1}`` and ``M = R_{d}^{-1}``. -""" -mutable struct KrylovSQD{T, F, Tv, Ta} <: AbstractKrylovSolver{T} - m::Int - n::Int - - f::F # Krylov function - atol::T # absolute tolerance - rtol::T # relative tolerance - - # Problem data - A::Ta # Constraint matrix - θ::Tv # scaling - regP::Tv # primal regularization - regD::Tv # dual regularization - - function KrylovSQD(f::Function, A::Ta; - atol::T=eps(T)^(3 // 4), - rtol::T=eps(T)^(3 // 4) - ) where{T, Ta <: AbstractMatrix{T}} - F = typeof(f) - m, n = size(A) - - return new{T, F, Vector{T}, Ta}( - m, n, - f, atol, rtol, - A, ones(T, n), ones(T, n), ones(T, m) - ) - end -end - -setup(::Type{KrylovSQD}, A; method, kwargs...) = KrylovSQD(method, A; kwargs...) - -backend(kkt::KrylovSQD) = "Krylov ($(kkt.f))" -linear_system(::KrylovSQD) = "Augmented system" - - -""" - update!(kkt::KrylovSQD, θ, regP, regD) - -Update diagonal scaling ``θ`` and primal-dual regularizations. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - # Update diagonal scaling - kkt.θ .= θ - # Update regularizers - kkt.regP .= regP - kkt.regD .= regD - - return nothing -end - -""" - solve!(dx, dy, kkt::KrylovSQD, ξp, ξd) - -Solve the augmented system using the selected Kyrlov method. -""" -function solve!( - 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(T) ./ (kkt.θ .+ kkt.regP) - N = Diagonal(d) - - M = Diagonal(one(T) ./ kkt.regD) - - # Solve augmented system - _dy, _dx, stats = kkt.f(kkt.A, ξp, ξd, - M=M, N=N, - atol=kkt.atol, rtol=kkt.rtol - ) - - dx .= _dx - dy .= _dy - - return nothing -end diff --git a/src/KKT/lapack.jl b/src/KKT/lapack.jl deleted file mode 100644 index b1090315..00000000 --- a/src/KKT/lapack.jl +++ /dev/null @@ -1,202 +0,0 @@ -using LinearAlgebra.LAPACK - -const BlasReal = LinearAlgebra.BlasReal - -""" - DenseSPD{T} - -KKT solver for dense linear systems. - -Uses a Cholesky factorization of the normal equations system. Not recommended - for systems of size larger than a few thousands. - -```julia -model = Tulip.Model{Float64}() -model.params.KKT.Factory = Tulip.Factory(Tulip.KKT.DenseSPD) -``` - -!!! info - Dispatches to BLAS/LAPACK in `Float32` and `Float64` arithmetic. -""" -mutable struct DenseSPD{T} <: AbstractKKTSolver{T} - m::Int # Number of rows in A - n::Int # Number of columns in A - - # Problem data - A::Matrix{T} - θ::Vector{T} - - # Regularizations - regP::Vector{T} # primal - regD::Vector{T} # dual - - _A::Matrix{T} # place-holder to avoid allocating memory afterwards - - # Factorization - S::Matrix{T} # Normal equations matrix, that also holds the factorization - - # Constructor - function DenseSPD(A::AbstractMatrix{T}) where{T} - m, n = size(A) - θ = ones(T, n) - - # We just need to allocate memory here, - # so no factorization yet - S = zeros(T, m, m) - - return new{T}(m, n, A, θ, zeros(T, n), zeros(T, m), Matrix(A), S) - end -end - -backend(::DenseSPD) = "LAPACK" -linear_system(::DenseSPD) = "Normal equations" - -""" - update!(kkt::DenseSPD{T}, θ, regP, regD) - -Compute normal equations system matrix and update Cholesky factorization. - -Uses Julia's generic linear algebra. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - kkt.θ .= θ - kkt.regP .= regP - kkt.regD .= regD - - # 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(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 - @inbounds for i in 1:kkt.m - kkt.S[i, i] += kkt.regD[i] - end - - # Cholesky factorization - cholesky!(Symmetric(kkt.S)) - - return nothing -end - -""" - update!(kkt::DenseSPD{<:BlasReal}, θ, regP, regD) - -Compute normal equations system matrix and update Cholesky factorization. - -Uses BLAS and LAPACK routines. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - kkt.θ .= θ - kkt.regP .= regP - kkt.regD .= regD - - # 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(T) ./ (kkt.θ .+ kkt.regP))) - rmul!(kkt._A, D) # B = A * √D - BLAS.syrk!('U', 'N', one(T), kkt._A, zero(T), kkt.S) - - # Regularization - # TODO: only update diagonal of S if only dual regularization changes - @inbounds for i in 1:kkt.m - kkt.S[i, i] += kkt.regD[i] - end - - # Cholesky factorization - # TODO: regularize if needed - LAPACK.potrf!('U', kkt.S) - - return nothing -end - -""" - solve!(dx, dy, kkt, ξp, ξd) - -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{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 - dy .= ξp .+ kkt.A * (ξd ./ (kkt.θ .+ kkt.regP)) - - # Solve normal equations - ldiv!(UpperTriangular(kkt.S)', dy) - ldiv!(UpperTriangular(kkt.S) , dy) - - # Recover dx - # TODO: use more efficient mul! syntax - dx .= (kkt.A' * dy - ξd) ./ (kkt.θ .+ kkt.regP) - - # TODO: Iterative refinement - return nothing -end - -""" - solve!(dx, dy, kkt, ξp, ξd) - -Solve the augmented system, overwriting `dx, dy` with the result. - -Uses one LAPACK call for solving the normal equations system. -""" -function solve!( - 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 - dy .= ξp .+ kkt.A * (ξd ./ (kkt.θ .+ kkt.regP)) - - # Solve augmented system - LAPACK.potrs!('U', kkt.S, dy) # potrs! over-writes the right-hand side - - # Recover dx - # TODO: use more efficient mul! syntax - dx .= (kkt.A' * dy - ξd) ./ (kkt.θ .+ kkt.regP) - - # TODO: Iterative refinement - return nothing -end \ No newline at end of file diff --git a/src/KKT/ldlfact.jl b/src/KKT/ldlfact.jl deleted file mode 100644 index 7c2d9902..00000000 --- a/src/KKT/ldlfact.jl +++ /dev/null @@ -1,153 +0,0 @@ -import LDLFactorizations -const LDLF = LDLFactorizations - -# ============================================================================== -# LDLFactSQD -# ============================================================================== - -""" - 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.KKT.Factory = Tulip.Factory(LDLFactSQD) -``` -""" -mutable struct LDLFactSQD{T<:Real} <: AbstractKKTSolver{T} - m::Int # Number of rows - n::Int # Number of columns - - # TODO: allow for user-provided ordering, - # and add flag (default to false) to know whether user ordering should be used - - # Problem data - A::SparseMatrixCSC{T, Int} - θ::Vector{T} # diagonal scaling - regP::Vector{T} # primal regularization - regD::Vector{T} # dual regularization - ξ::Vector{T} # right-hand side of the augmented system - - # Left-hand side matrix - S::SparseMatrixCSC{T, Int} - - # Factorization - F::LDLF.LDLFactorization{T} - - # TODO: constructor with initial memory allocation - function LDLFactSQD(A::AbstractMatrix{T}) where{T<:Real} - m, n = size(A) - θ = ones(T, n) - regP = ones(T, n) - regD = ones(T, m) - ξ = zeros(T, n+m) - - S = [ - spdiagm(0 => -θ) A'; - spzeros(T, m, n) spdiagm(0 => ones(T, m)) - ] - - # TODO: PSD-ness checks - # TODO: symbolic factorization only - F = LDLF.ldl_analyze(Symmetric(S)) - - return new{T}(m, n, A, θ, regP, regD, ξ, S, F) - end - -end - -setup(::Type{LDLFactSQD}, A) = LDLFactSQD(A) - -backend(::LDLFactSQD) = "LDLFactorizations.jl" -linear_system(::LDLFactSQD) = "Augmented system" - -""" - update!(kkt, θ, regP, regD) - -Update LDLᵀ factorization of the augmented system. - -Update diagonal scaling ``\\theta``, primal-dual regularizations, and re-compute - the factorization. -Throws a `PosDefException` if matrix is not quasi-definite. -""" -function update!( - 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)." - )) - length(regP) == kkt.n || throw(DimensionMismatch( - "regP has length $(length(regP)) but linear solver has n=$(kkt.n)" - )) - length(regD) == kkt.m || throw(DimensionMismatch( - "regD has length $(length(regD)) but linear solver has m=$(kkt.m)" - )) - - # Update diagonal scaling - kkt.θ .= θ - # Update regularizers - kkt.regP .= regP - kkt.regD .= regD - - # Update S. - # S is stored as upper-triangular and only its diagonal changes. - @inbounds for j in 1:kkt.n - k = kkt.S.colptr[1+j] - 1 - kkt.S.nzval[k] = -kkt.θ[j] - regP[j] - end - @inbounds for i in 1:kkt.m - k = kkt.S.colptr[1+kkt.n+i] - 1 - kkt.S.nzval[k] = regD[i] - end - - # TODO: PSD-ness checks - try - LDLF.ldl_factorize!(Symmetric(kkt.S), kkt.F) - catch err - isa(err, LDLF.SQDException) && throw(PosDefException(-1)) - rethrow(err) - end - - return nothing -end - -""" - solve!(dx, dy, kkt, ξp, ξd) - -Solve the augmented system, overwriting `dx, dy` with the result. -""" -function solve!( - 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 the right-hand side - @inbounds for i in 1:n - kkt.ξ[i] = ξd[i] - end - @inbounds for i in 1:m - kkt.ξ[i+n] = ξp[i] - end - - # Solve the augmented system - LDLF.ldiv!(kkt.F, kkt.ξ) - - # Recover dx, dy - @inbounds for i in 1:n - dx[i] = kkt.ξ[i] - end - @inbounds for i in 1:m - dy[i] = kkt.ξ[i+n] - end - - # TODO: Iterative refinement - return nothing -end diff --git a/src/KKT/systems.jl b/src/KKT/systems.jl new file mode 100644 index 00000000..dd3ef6f3 --- /dev/null +++ b/src/KKT/systems.jl @@ -0,0 +1,54 @@ +""" + DefaultKKTSystem + +Default KKT system setting. Currently equivalent to [`K2`](@ref) +""" +struct DefaultKKTSystem <: AbstractKKTSystem end + +@doc raw""" + K2 <: AbstractKKTSystem + +Augmented system +```math + \begin{bmatrix} + -(\Theta^{-1} + R_{p}) & A^{\top}\\ + A & R_{d} + \end{bmatrix} + \begin{bmatrix} + \Delta x\\ + \Delta y + \end{bmatrix} + = + \begin{bmatrix} + \xi_d\\ + \xi_p + \end{bmatrix} +``` +where +* ``\Theta^{-1} = X^{-l}Z^{l} + X^{-u} Z^{u}`` +* ``R_{p}, R_{d}`` are current primal and dual regularizations +* ``\xi_{d}, \xi_{p}`` are given right-hand sides +""" +struct K2 <: AbstractKKTSystem end + +@doc raw""" + K1 <: AbstractKKTSystem + +Normal equations system +```math + \begin{array}{rl} + \left( + A (\Theta^{-1} + R_{p})^{-1} A^{\top} + R_{d} + \right) + \Delta y + & = + \xi_{p} + A (Θ^{-1} + R_{p})^{-1} \xi_{d}\\ + \Delta x &= (Θ^{-1} + R_{p})^{-1} (A^{\top} \Delta y - \xi_{d}) + \end{array} +``` +where +* ``\Theta^{-1} = X^{-l}Z^{l} + X^{-u} Z^{u}`` +* ``R_{p}, R_{d}`` are current primal and dual regularizations +* ``\xi_{d}, \xi_{p}`` are the augmented system's right-hand side +""" +struct K1 <: AbstractKKTSystem end diff --git a/src/Tulip.jl b/src/Tulip.jl index c0a52ea8..bcfe2f60 100644 --- a/src/Tulip.jl +++ b/src/Tulip.jl @@ -7,7 +7,7 @@ using SparseArrays using TimerOutputs -version() = v"0.7.5" +version() = v"0.8.0" include("utils.jl") diff --git a/test/KKT/cholmod.jl b/test/KKT/Cholmod/cholmod.jl similarity index 65% rename from test/KKT/cholmod.jl rename to test/KKT/Cholmod/cholmod.jl index 633fcfa3..787df38a 100644 --- a/test/KKT/cholmod.jl +++ b/test/KKT/Cholmod/cholmod.jl @@ -6,13 +6,12 @@ ]) @testset "LDL" begin - kkt = KKT.CholmodSolver(A, normal_equations=false) + kkt = KKT.setup(A, KKT.K2(), KKT.TlpCholmod.Backend()) KKT.run_ls_tests(A, kkt) end - + @testset "Cholesky" begin - kkt = KKT.CholmodSolver(A, normal_equations=true) + kkt = KKT.setup(A, KKT.K1(), KKT.TlpCholmod.Backend()) KKT.run_ls_tests(A, kkt) end - -end \ No newline at end of file +end diff --git a/test/KKT/lapack.jl b/test/KKT/Dense/lapack.jl similarity index 76% rename from test/KKT/lapack.jl rename to test/KKT/Dense/lapack.jl index d58b9a98..157af29a 100644 --- a/test/KKT/lapack.jl +++ b/test/KKT/Dense/lapack.jl @@ -5,8 +5,8 @@ 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.DenseSPD(A) + kkt = KKT.setup(A, KKT.K1(), KKT.TlpDense.Backend()) KKT.run_ls_tests(A, kkt) end end -end \ No newline at end of file +end diff --git a/test/KKT/KKT.jl b/test/KKT/KKT.jl index 8f3b8137..db884bd8 100644 --- a/test/KKT/KKT.jl +++ b/test/KKT/KKT.jl @@ -1,6 +1,6 @@ const KKT = Tulip.KKT -include("lapack.jl") -include("cholmod.jl") -include("ldlfact.jl") -include("krylov.jl") \ No newline at end of file +include("Dense/lapack.jl") +include("Cholmod/cholmod.jl") +include("LDLFactorizations/ldlfact.jl") +include("Krylov/krylov.jl") diff --git a/test/KKT/Krylov/krylov.jl b/test/KKT/Krylov/krylov.jl new file mode 100644 index 00000000..56eceab0 --- /dev/null +++ b/test/KKT/Krylov/krylov.jl @@ -0,0 +1,7 @@ +using Krylov + +@testset "Krylov" begin + include("spd.jl") + include("sid.jl") + include("sqd.jl") +end diff --git a/test/KKT/Krylov/sid.jl b/test/KKT/Krylov/sid.jl new file mode 100644 index 00000000..d594f5ff --- /dev/null +++ b/test/KKT/Krylov/sid.jl @@ -0,0 +1,19 @@ +function test_krylov_sid(T, ksolver) + A = SparseMatrixCSC{T, Int}([ + 1 0 1 0; + 0 1 0 1 + ]) + + kkt = KKT.setup(A, KKT.K2(), KKT.TlpKrylov.Backend(ksolver, Vector{T})) + KKT.run_ls_tests(A, kkt) + + return nothing +end + +@testset "SID" begin + for T in TvTYPES, ksolver in [MinresSolver, MinresQlpSolver, SymmlqSolver] + @testset "$ksolver ($T)" begin + test_krylov_sid(T, ksolver) + end + end +end diff --git a/test/KKT/Krylov/spd.jl b/test/KKT/Krylov/spd.jl new file mode 100644 index 00000000..ce13b9f9 --- /dev/null +++ b/test/KKT/Krylov/spd.jl @@ -0,0 +1,19 @@ +function test_krylov_spd(T, ksolver) + A = SparseMatrixCSC{T, Int}([ + 1 0 1 0; + 0 1 0 1 + ]) + + kkt = KKT.setup(A, KKT.K1(), KKT.TlpKrylov.Backend(ksolver, Vector{T})) + KKT.run_ls_tests(A, kkt) + + return nothing +end + +@testset "SPD" begin + for T in TvTYPES, ksolver in [CgSolver, MinresSolver] + @testset "$ksolver ($T)" begin + test_krylov_spd(T, ksolver) + end + end +end diff --git a/test/KKT/Krylov/sqd.jl b/test/KKT/Krylov/sqd.jl new file mode 100644 index 00000000..edfe401f --- /dev/null +++ b/test/KKT/Krylov/sqd.jl @@ -0,0 +1,19 @@ +function test_krylov_sqd(T, ksolver) + A = SparseMatrixCSC{T, Int}([ + 1 0 1 0; + 0 1 0 1 + ]) + + kkt = KKT.setup(A, KKT.K2(), KKT.TlpKrylov.Backend(ksolver, Vector{T})) + KKT.run_ls_tests(A, kkt) + + return nothing +end + +@testset "SQD" begin + for T in TvTYPES, ksolver in [TricgSolver, TrimrSolver] + @testset "$ksolver ($T)" begin + test_krylov_sqd(T, ksolver) + end + end +end diff --git a/test/KKT/ldlfact.jl b/test/KKT/LDLFactorizations/ldlfact.jl similarity index 76% rename from test/KKT/ldlfact.jl rename to test/KKT/LDLFactorizations/ldlfact.jl index c5ddea23..6d0c28c1 100644 --- a/test/KKT/ldlfact.jl +++ b/test/KKT/LDLFactorizations/ldlfact.jl @@ -5,8 +5,8 @@ 1 0 1 0; 0 1 0 1 ]) - kkt = KKT.LDLFactSQD(A) + kkt = KKT.setup(A, KKT.K2(), KKT.TlpLDLFact.Backend()) KKT.run_ls_tests(A, kkt) end end -end \ No newline at end of file +end diff --git a/test/KKT/krylov.jl b/test/KKT/krylov.jl deleted file mode 100644 index bf94fae6..00000000 --- a/test/KKT/krylov.jl +++ /dev/null @@ -1,41 +0,0 @@ -import Krylov - -@testset "Krylov" begin - for T in TvTYPES - @testset "$T" begin - # Kyrlov solvers for symmetric positive-definite systems - for f in [Krylov.cg, Krylov.minres, Krylov.minres_qlp] - @testset "PD: $f" begin - A = SparseMatrixCSC{T, Int}([ - 1 0 1 0; - 0 1 0 1 - ]) - kkt = KKT.KrylovSPD(f, A) - KKT.run_ls_tests(A, kkt) - end - end - # Kyrlov solvers for symmetric indefinite systems - for f in [Krylov.minres, Krylov.minres_qlp] - @testset "SID: $f" begin - A = SparseMatrixCSC{T, Int}([ - 1 0 1 0; - 0 1 0 1 - ]) - kkt = KKT.KrylovSID(f, A) - KKT.run_ls_tests(A, kkt) - end - end - # Kyrlov solvers for symmetric quasi-definite systems - for f in [Krylov.tricg, Krylov.trimr] - @testset "SQD: $f" begin - A = SparseMatrixCSC{T, Int}([ - 1 0 1 0; - 0 1 0 1 - ]) - kkt = KKT.KrylovSQD(f, A) - KKT.run_ls_tests(A, kkt) - end - end - end - end -end diff --git a/test/Project.toml b/test/Project.toml index f9ac870d..681614ca 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -9,3 +9,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Convex = "0.14" +Krylov = "0.7.3"