From 01459163ed27fb5fcda1adbb9cf924486d55862c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 26 May 2025 18:07:14 +0200 Subject: [PATCH 1/8] Refactor model --- Project.toml | 8 + docs/src/low-rank_data.md | 7 +- examples/k.jl | 3 +- examples/simple_linear.jl | 27 ++ examples/solve_sdpa.jl | 13 + src/Loraine.jl | 50 +++- src/Solvers.jl | 542 +++++++++++++++++++------------------ src/initial_point.jl | 78 +++--- src/kron_etc.jl | 11 +- src/makeBBBB.jl | 228 ---------------- src/model.jl | 287 -------------------- src/predictor_corrector.jl | 210 +++++++------- src/prepare_W.jl | 45 +-- test/MOI_wrapper.jl | 15 +- 14 files changed, 540 insertions(+), 984 deletions(-) create mode 100644 examples/simple_linear.jl delete mode 100644 src/makeBBBB.jl delete mode 100644 src/model.jl diff --git a/Project.toml b/Project.toml index 80fbf81..d836f16 100644 --- a/Project.toml +++ b/Project.toml @@ -10,11 +10,15 @@ FameSVD = "9ba2d756-9ce3-11e9-1a71-0ffcb019784d" GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" @@ -25,8 +29,12 @@ DocumenterTools = "0.1.21" FameSVD = "0.1" GenericLinearAlgebra = "0.3" JuMP = "1" +LowRankOpt = "0.2.1" MathOptInterface = "1" MultiFloats = "2" +MutableArithmetics = "1.6.4" +NLPModels = "0.21.5" Revise = "3" +SolverCore = "0.3.8" TimerOutputs = "0.5" julia = "1" diff --git a/docs/src/low-rank_data.md b/docs/src/low-rank_data.md index 3302c35..0a1916c 100644 --- a/docs/src/low-rank_data.md +++ b/docs/src/low-rank_data.md @@ -9,9 +9,4 @@ with data matrices ``B_i\in{\mathbb R}^{m\times k}`` and with ``k \ll n``. The complexity of ``H`` assembly is then reduced from ``nm^3`` to ``knm^2``. *In particular:* -If we know that ``A_i`` have rank one, the decomposition ``A_i = b_i b_i^\top`` is performed by Loraine automatically (`datarank = -1`). - -The model is represented internally via the following `struct`: -```@docs -Loraine.MyModel -``` +If we know that ``A_i`` have rank one, the decomposition ``A_i = b_i b_i^\top`` is performed by [LowRankOpt](https://github.com/blegat/LowRankOpt.jl/) automatically (`detect_rank_1_tol = 1e-6`). diff --git a/examples/k.jl b/examples/k.jl index ced0253..1a6e186 100644 --- a/examples/k.jl +++ b/examples/k.jl @@ -5,7 +5,7 @@ import Loraine using MultiFloats # model = Model(Loraine.Optimizer) -model = Model(Loraine.Optimizer{Float64x2}) +model = GenericModel{Float64x2}(Loraine.Optimizer{Float64x2}) # model = JuMP.GenericModel{Float64x2}(Loraine.Optimizer{Float64x2}) # @variable(model, x >= 0) @@ -30,6 +30,7 @@ using Test @test primal_status(model) == MOI.FEASIBLE_POINT @test dual_status(model) == MOI.FEASIBLE_POINT @test objective_value(model) ≈ 4 rtol = 1e-6 +@test dual_objective_value(model) ≈ 4 rtol = 1e-6 @test value(x) ≈ 2 rtol = 1e-6 diff --git a/examples/simple_linear.jl b/examples/simple_linear.jl new file mode 100644 index 0000000..00bbff1 --- /dev/null +++ b/examples/simple_linear.jl @@ -0,0 +1,27 @@ +model = Model(Loraine.Optimizer) +@variable(model, x) +@constraint(model, con_ref, x <= 2) +@objective(model, Max, 3x) +optimize!(model) +@test value(x) ≈ 2 rtol = 1e-6 +@test objective_value(model) ≈ 6 rtol = 1e-6 +@test termination_status(model) == MOI.OPTIMAL +@test primal_status(model) == MOI.FEASIBLE_POINT +@test dual_status(model) == MOI.FEASIBLE_POINT +@test dual(con_ref) ≈ -3 rtol = 1e-6 + +model = Model(Loraine.Optimizer) +@variable(model, x) +@variable(model, y) +@constraint(model, cx, x >= 0) +@constraint(model, cy, y <= 0) +@objective(model, Min, x - y) +optimize!(model) +@test value(x) ≈ 0 atol = 1e-6 +@test value(y) ≈ 0 atol = 1e-6 +@test objective_value(model) ≈ 0 atol = 1e-6 +@test termination_status(model) == MOI.OPTIMAL +@test primal_status(model) == MOI.FEASIBLE_POINT +@test dual_status(model) == MOI.FEASIBLE_POINT +@test dual(cx) ≈ 1 rtol = 1e-6 +@test dual(cy) ≈ -1 rtol = 1e-6 diff --git a/examples/solve_sdpa.jl b/examples/solve_sdpa.jl index 19695f1..fc191f6 100644 --- a/examples/solve_sdpa.jl +++ b/examples/solve_sdpa.jl @@ -62,6 +62,19 @@ using Test @test objective_value(model) ≈ 23 rtol = 1e-6 # # value.(X) +# With CG now +set_attribute(model, "kit", 1) +optimize!(model) +@test objective_value(model) ≈ 23 rtol = 1e-6 + +set_attribute(model, "preconditioner", 0) +optimize!(model) +@test objective_value(model) ≈ 23 rtol = 1e-6 + +set_attribute(model, "preconditioner", 2) +optimize!(model) +@test objective_value(model) ≈ 23 rtol = 1e-6 + # Mosek (CSDP, etc) for a comparison # Mosek must solve the dualized problem to be efficient diff --git a/src/Loraine.jl b/src/Loraine.jl index ca6e89d..5da368f 100644 --- a/src/Loraine.jl +++ b/src/Loraine.jl @@ -16,17 +16,61 @@ using FameSVD # using MKLSparse # using MKL +import MathOptInterface as MOI +import LowRankOpt as LRO +struct Optimizer{T} + dummy::T +end +function Optimizer{T}() where {T} + model = LRO.Optimizer{T}() + MOI.set( + model, + MOI.RawOptimizerAttribute("solver"), + Solvers.Solver{T}, + ) + return model +end +Optimizer() = Optimizer{Float64}() + #modules include("Solvers.jl") using .Solvers include("kron_etc.jl") -include("makeBBBB.jl") include("initial_point.jl") include("predictor_corrector.jl") include("prepare_W.jl") -include("MOI_wrapper.jl") -include("multithreading.jl") + +function prepare_model_data(d,drank) + msizes = Vector{Int64} + n = Int64(get(d, "nvar", 1)); + msizesa = get(d, "msizes", 1) + if length(msizesa) == 1 + msizes = [convert.(Int64,msizesa)] + else + msizes = convert.(Int64,msizesa[:]) + end + nlin = Int64(get(d, "nlin", 1)) + nlmi = Int64(get(d, "nlmi", 1)) + A = get(d, "A", 1); + @assert size(A, 1) == nlmi + b = -get(d, "c", 1); + @assert length(b) == n + b_const = -get(d, "b_const", 1); + + if nlin > 0 + d_lin = -get(d, "d", 1) + d_lin = d_lin[:] + C_lin = -get(d, "C", 1) + else + d_lin = sparse([0.; 0.]) + C_lin = sparse([0. 0.;0. 0.]) + end + + model = LRO.Model(A[:,2:end], _prepare_A(A,drank,κ)..., b, b_const, d_lin, C_lin, msizes) + + return model +end function loraine(d, options::Dict) diff --git a/src/Solvers.jl b/src/Solvers.jl index f89b2ec..4029e89 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -8,14 +8,55 @@ import Statistics: mean using Printf using TimerOutputs using MultiFloats +import NLPModels +import SolverCore +import LowRankOpt as LRO # using MKLSparse # using MKL include("kron_etc.jl") -include("makeBBBB.jl") -include("model.jl") -mutable struct MySolver{T} +struct FactoredMatrix{T} <: AbstractMatrix{T} + factor::Matrix{T} + factor_inv::Matrix{T} # inv(factor) + matrix::Matrix{T} # factor * factor_inv' +end +Base.size(A::FactoredMatrix) = size(A.matrix) +Base.getindex(A::FactoredMatrix, i, j) = Base.getindex(A.matrix, i, j) +# Multiply through the dense `matrix` field so we hit BLAS instead of the +# generic `AbstractMatrix` (getindex-based) fallback. +Base.:*(A::FactoredMatrix, B::AbstractMatrix) = A.matrix * B +Base.:*(A::AbstractMatrix, B::FactoredMatrix) = A * B.matrix +Base.:*(A::FactoredMatrix, B::FactoredMatrix) = A.matrix * B.matrix +function LinearAlgebra.mul!( + C::AbstractMatrix, + A::FactoredMatrix, + B::AbstractVecOrMat, + α::Number, + β::Number, +) + return LinearAlgebra.mul!(C, A.matrix, B, α, β) +end +function LinearAlgebra.mul!( + C::AbstractMatrix, + A::AbstractVecOrMat, + B::FactoredMatrix, + α::Number, + β::Number, +) + return LinearAlgebra.mul!(C, A, B.matrix, α, β) +end +function LinearAlgebra.mul!( + C::AbstractMatrix, + A::FactoredMatrix, + B::FactoredMatrix, + α::Number, + β::Number, +) + return LinearAlgebra.mul!(C, A.matrix, B.matrix, α, β) +end + +mutable struct MySolver{T,M} # main options kit::Int64 tol_cg::T @@ -36,7 +77,7 @@ mutable struct MySolver{T} to::Any # model and preprocessed model data - model::MyModel + model::M predict::Bool @@ -47,22 +88,18 @@ mutable struct MySolver{T} expon::T iter::Int64 DIMACS_error::T + BBBB::Matrix{T} cholBBBB status::Int regcount::Int - err1::T - err2::T - err3::T - err4::T - err5::T - err6::T + err::NTuple{6,T} - X - S - y + X::LRO.VectorizedSolution{T} + S::LRO.VectorizedSolution{T} + y::Vector{T} yold delX delS @@ -70,26 +107,21 @@ mutable struct MySolver{T} Xn Sn - X_lin - S_lin - Si_lin - S_lin_inv - delX_lin - delS_lin - Xn_lin - Sn_lin + Si_lin::Vector{T} + S_lin_inv::Vector{T} + delX_lin::Vector{T} + delS_lin::Vector{T} + Xn_lin::Vector{T} + Sn_lin::Vector{T} D - G - Gi - W + W::LRO.ShapedSolution{T,FactoredMatrix{T}} Si DDsi - Rp - Rd + Rp::Vector{T} + Rd::LRO.VectorizedSolution{T} Rc - Rd_lin cg_iter_pre cg_iter_cor @@ -106,6 +138,8 @@ mutable struct MySolver{T} RNT RNT_lin + y_buffer::Vector{T} # Bufer of the same size as `y` + function MySolver{T}( kit::Int64, tol_cg::Float64 , @@ -122,10 +156,11 @@ mutable struct MySolver{T} timing::Int64, maxit::Int64, datasparsity::Int64, - model::MyModel - ) where {T} + model::LRO.Model{T,A} + ) where {T,A} - solver = new{T}() + buffered_model = LRO.BufferedModelForSchur(model, datasparsity) + solver = new{T,typeof(buffered_model)}() solver.kit = kit solver.tol_cg = tol_cg solver.tol_cg_up = tol_cg_up @@ -141,7 +176,7 @@ mutable struct MySolver{T} solver.timing = timing solver.maxit = maxit solver.datasparsity = datasparsity - solver.model = model + solver.model = buffered_model return solver end end @@ -161,6 +196,53 @@ mutable struct Halpha end end +struct Solver{T,M} <: SolverCore.AbstractOptimizationSolver + solver::MySolver{T,M} + halpha::Halpha + stats::SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any} +end + +function Solver{T}(model::LRO.Model; kws...) where {T} + options = Dict(kw[1] => kw[2] for kw in kws) + stats = SolverCore.GenericExecutionStats(model) + solver, halpha = load(model, options; T) + solver.X = LRO.VectorizedSolution{T}(stats.solution, model.dim) + solver.y = stats.multipliers + Solver(solver, halpha, stats) +end + +function LRO.MOI.get(solver::Solver, ::LRO.Solution) + return LRO.VectorizedSolution(solver.stats.solution, solver.solver.model.model.dim) +end + +const STATUS_MAP = [ + :unknown, + :first_order, + :unbounded, # the primal `max ⟨C,X⟩` has an unbounded ray so the primal is unbounded + :infeasible, # the dual `max ⟨b,y⟩` has an unbounded ray so the primal `min ⟨C,X⟩` is infeasible + :max_iter, + :exception, +] + +function SolverCore.solve!( + solver::Solver, + model::NLPModels.AbstractNLPModel; # Same as `solver.model`, we can ignore + kws..., +) + for kw in kws + field = Symbol(kw[1]) + if field == :verbose + field = :verb # TODO rename to :verbose to follow NLPModels convention + end + setproperty!(solver.solver, field, kw[2]) + end + solver.solver.to = TimerOutput() + solve(solver.solver, solver.halpha) + solver.stats.status = STATUS_MAP[solver.solver.status + 1] + solver.stats.objective = NLPModels.obj(solver.solver.model, solver.solver.X) + return +end + include("initial_point.jl") include("prepare_W.jl") @@ -200,7 +282,7 @@ function load(model, options::Dict; T = Float64) initpoint = Int64(get(options, "initpoint", 0)) timing = Int64(get(options, "timing", 1)) maxit = Int64(get(options, "maxit", 100)) - datasparsity = Int64(get(options, "maxit", 8)) + datasparsity = Int64(get(options, "datasparsity", 8)) solver = MySolver{T}(kit, tol_cg, @@ -217,49 +299,12 @@ function load(model, options::Dict; T = Float64) timing, maxit, datasparsity, - MyModel(model.A, - model.AA, - model.B, - model.C, - model.nzA, - model.sigmaA, - model.qA, - model.b, - model.b_const, - model.d_lin, - model.C_lin, - model.n, - model.msizes, - model.nlin, - model.nlmi - ) + model, ) halpha = Halpha(kit) solver.cg_iter_tot = 0 - if verb > 0 - t1 = time() - @printf("\n *** Loraine.jl v0.2.6 ***\n") - @printf(" *** Initialisation STARTS\n") - end - - if verb > 0 - @printf(" Number of variables: %5d\n",model.n) - @printf(" LMI constraints : %5d\n",model.nlmi) - if model.nlmi>0 - @printf(" Matrix size(s) :") - Printf.format.(Ref(stdout), Ref(Printf.Format("%6d")), model.msizes); - @printf("\n") - end - @printf(" Linear constraints : %5d\n",model.nlin) - if solver.kit>0 - @printf(" Preconditioner : %5d\n",preconditioner) - else - @printf(" Preconditioner : none, using direct solver\n") - end - end - # Input parameters check if kit < 0 || kit > 1 solver.kit = 0 @@ -304,6 +349,23 @@ end function solve(solver::MySolver,halpha::Halpha) t1 = time() if solver.verb > 0 + @printf("\n *** Loraine.jl v0.2.5 ***\n") + + @printf(" Number of variables: %5d\n",solver.model.meta.ncon) + @printf(" LMI constraints : %5d\n",LRO.num_matrices(solver.model)) + if LRO.num_matrices(solver.model) > 0 + @printf(" Matrix size(s) :") + msizes = LRO.side_dimension.(Ref(solver.model), LRO.matrix_indices(solver.model)) + Printf.format.(Ref(stdout), Ref(Printf.Format("%6d")), msizes); + @printf("\n") + end + @printf(" Linear constraints : %5d\n",LRO.num_scalars(solver.model)) + if solver.kit>0 + @printf(" Preconditioner : %5d\n",solver.preconditioner) + else + @printf(" Preconditioner : none, using direct solver\n") + end + @printf(" *** IP STARTS\n") if solver.verb < 2 if solver.kit == 0 @@ -338,7 +400,7 @@ function solve(solver::MySolver,halpha::Halpha) if solver.preconditioner == 4 # if (cg_iter2>erank*nlmi*sqrt(n)/1 && iter>sqrt(n)/60)||cg_iter2>100 %for SNL problems - if (solver.cg_iter_cor / 2 > solver.erank * solver.model.nlmi * sqrt(solver.model.n)/20 && solver.iter > sqrt(solver.model.n) / 60) || solver.cg_iter_cor > 100 + if (solver.cg_iter_cor / 2 > solver.erank * LRO.num_matrices(solver.model) * sqrt(solver.model.meta.ncon)/20 && solver.iter > sqrt(solver.model.meta.ncon) / 60) || solver.cg_iter_cor > 100 solver.preconditioner = 1; solver.aamat = 2; if solver.verb > 0 println("Switching to preconditioner 1") @@ -362,25 +424,29 @@ end function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} - solver.X = Matrix{T}[] - solver.S = Matrix{T}[] - solver.y = Vector{T}[] + solver.S = similar(solver.X) + solver.Rd = similar(solver.X) + solver.Rp = zeros(T, solver.model.meta.ncon) + solver.y_buffer = similar(solver.Rp) solver.delX = Matrix{T}[] solver.delS = Matrix{T}[] solver.D = Vector{T}[] - solver.G = Matrix{T}[] - solver.Gi = Matrix{T}[] - solver.W = Matrix{T}[] + solver.W = LRO.ShapedSolution{T,FactoredMatrix{T}}( + zeros(T, LRO.num_scalars(solver.model)), + map(LRO.matrix_indices(solver.model)) do i + dim = LRO.side_dimension(solver.model, i) + FactoredMatrix(zeros(T, dim, dim), zeros(T, dim, dim), zeros(T, dim, dim)) + end, + ) solver.Si = Matrix{T}[] solver.DDsi = Vector{T}[] - solver.Rd = Matrix{T}[] solver.Rc = Matrix{T}[] - solver.alpha = zeros(solver.model.nlmi) - solver.beta = zeros(solver.model.nlmi) + solver.alpha = zeros(LRO.num_matrices(solver.model)) + solver.beta = zeros(LRO.num_matrices(solver.model)) solver.Xn = Matrix{T}[] solver.Sn = Matrix{T}[] @@ -388,70 +454,60 @@ function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} solver.regcount = 0 - for i = 1:solver.model.nlmi - push!(solver.X,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.S,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.delX,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.delS,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.D, zeros(solver.model.msizes[i])) - push!(solver.G,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.Gi,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.W,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.Si,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.DDsi,zeros(solver.model.msizes[i])) - push!(solver.Rd,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.Rc,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.Xn,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.Sn,zeros(solver.model.msizes[i], solver.model.msizes[i])) - push!(solver.RNT,zeros(solver.model.msizes[i], solver.model.msizes[i])) + solver.delS_lin = zeros(LRO.num_scalars(solver.model)) + for mat_idx in LRO.matrix_indices(solver.model) + dim = LRO.side_dimension(solver.model, mat_idx) + push!(solver.delX,zeros(dim, dim)) + push!(solver.delS,zeros(dim, dim)) + push!(solver.D, zeros(dim)) + push!(solver.Si,zeros(dim, dim)) + push!(solver.DDsi,zeros(dim)) + push!(solver.Rc,zeros(dim, dim)) + push!(solver.Xn,zeros(dim, dim)) + push!(solver.Sn,zeros(dim, dim)) + push!(solver.RNT,zeros(dim, dim)) end halpha.Umat = Matrix{T}[] halpha.Z = Matrix{T}[] halpha.AAAATtau = SparseMatrixCSC{T}[] + ncon = solver.model.meta.ncon - for i = 1:solver.model.nlmi - push!(halpha.Umat,zeros(solver.model.msizes[i], solver.erank)) - push!(halpha.Z,zeros(solver.model.msizes[i], solver.model.msizes[i])) - # tmp = Matrix(I(solver.model.msizes[i])) + for mat_idx in LRO.matrix_indices(solver.model) + dim = LRO.side_dimension(solver.model, mat_idx) + push!(halpha.Umat,zeros(dim, solver.erank)) + push!(halpha.Z,zeros(dim, dim)) + # tmp = Matrix(I(dim)) # push!(halpha.cholS,cholesky(tmp)) - push!(halpha.AAAATtau,spzeros(solver.model.n, solver.model.n)) + push!(halpha.AAAATtau, spzeros(ncon, ncon)) end if solver.kit == 1 - if solver.model.nlmi == 0 + if LRO.num_matrices(solver.model) == 0 if solver.verb > 0 - println("WARNING: Switching to a direct solver, no LMIs") + @warn("Switching to a direct solver, no LMIs") end solver.kit = 0 - elseif solver.model.nlmi > 0 && solver.erank >= maximum(solver.model.msizes) - 1 + elseif LRO.num_matrices(solver.model) > 0 && solver.erank >= maximum(Base.Fix1(LRO.side_dimension, solver.model), LRO.matrix_indices(solver.model)) - 1 if solver.verb > 0 - println("WARNING: Switching to a direct solver, erank bigger than matrix size") + @warn("Switching to a direct solver, erank bigger than matrix size") end solver.kit = 0 end end - # when datarank was set to -1 and conversion failed, we switch to datarank = 0 - if ~isempty(solver.model.B) - if solver.model.nlmi > 0 - for ilmi = 1:solver.model.nlmi - if nnz(solver.model.B[ilmi]) == 0 - solver.datarank = 0 - end - end - end + if solver.kit == 0 # if direct solver; compute the Hessian matrix + solver.BBBB = zeros(T, ncon, ncon) end end function myIPstep(solver::MySolver{T},halpha::Halpha) where {T} - mmm = Matrix{T}(undef, solver.model.n, solver.model.n) solver.iter += 1 if solver.iter > solver.maxit solver.status = 4 if solver.verb > 0 - println("WARNING: Stopped by iteration limit (stopping status = 4)") + @warn("Stopped by iteration limit (stopping status = 4)") end end solver.cg_iter_pre = 0 @@ -478,90 +534,68 @@ function myIPstep(solver::MySolver{T},halpha::Halpha) where {T} end function find_mu(solver) - trXS = 0 - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - trXS = trXS + sum(sum(solver.X[i] .* solver.S[i])) - end - end - mu = trXS - - if solver.model.nlin > 0 - mu = mu + tr(solver.X_lin' * solver.S_lin) - end - solver.mu = mu / (sum(solver.model.msizes) + solver.model.nlin) + solver.mu = dot(solver.X, solver.S) / (LRO.num_scalars(solver.model) + sum(Base.Fix1(LRO.side_dimension, solver.model), LRO.matrix_indices(solver.model), init = 0)) return solver.mu end function check_convergence(solver) # DIMACS error evaluation - solver.err1 = norm(solver.Rp) / (1 + norm(solver.model.b)) - (solver.err2,solver.err3,solver.err4,solver.err5,solver.err6) = [0.,0.,0.,0.,0.] - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - solver.err2 = solver.err2 + max(0, -eigmin(solver.X[i]) / (1 + norm(solver.model.b))) - solver.err3 = solver.err3 + norm(solver.Rd[i], 2) / (1 + norm(solver.model.C[i])) - solver.err4 = solver.err4 + max(0, -eigmin(solver.S[i]) / (1 + norm(solver.model.C[i]))) - # err5 = err5 + (vecC[i]"*vec(X[i])-b'*y)/(1+abs(vecC[i]'*vec(X[i]))+abs(b"*y)) - solver.err6 = solver.err6 + (vec(solver.S[i]))' * vec(solver.X[i]) / (1 + abs(vec(solver.model.C[i])' * vec(solver.X[i])) + abs(dot(solver.model.b', solver.y))) - end - end - - solver.err5 = (btrace(solver.model.nlmi, solver.model.C, solver.X) - dot(solver.model.b', solver.y)) / (1 + abs(btrace(solver.model.nlmi, solver.model.C, solver.X)) + abs(dot(solver.model.b', solver.y))) - if solver.model.nlin > 0 - solver.err2 = solver.err2 + max(0, -minimum(solver.X_lin) / (1 + norm(solver.model.b))) - solver.err3 = solver.err3 + norm(solver.Rd_lin) / (1 + norm(solver.model.d_lin)) - solver.err4 = solver.err4 + max(0, -minimum(solver.S_lin) / (1 + norm(solver.model.d_lin))) - solver.err5 = (btrace(solver.model.nlmi, solver.model.C, solver.X) + dot(solver.model.d_lin', solver.X_lin) - dot(solver.model.b',solver.y)) / (1 + abs(btrace(solver.model.nlmi, solver.model.C, solver.X)) + abs(dot(solver.model.b', solver.y))) - solver.err6 = solver.err6 + dot(solver.S_lin' , solver.X_lin) / (1 + abs(dot(solver.model.d_lin', solver.X_lin)) + abs(dot(solver.model.b', solver.y))) - end + pobj = NLPModels.obj(solver.model, solver.X) + dobj = LRO.dual_obj(solver.model, solver.y) + solver.err = LRO.errors( + solver.model, + solver.X; + y = solver.y, + primal_err = solver.Rp, + dual_slack = solver.S, + dual_err = solver.Rd, + pobj, + dobj, + ) - if solver.model.nlmi > 0 - DIMACS_error = solver.err1 + solver.err2 + solver.err3 + solver.err4 + abs(solver.err5) + solver.err6 - else - DIMACS_error = solver.err2 + solver.err3 + solver.err4 + abs(solver.err5) + solver.err6 - end + # `err[5]` may be negative so we need `abs` + DIMACS_error = sum(abs, solver.err) if solver.verb > 0 && solver.status == 0 - #@sprintf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.0d %9.0d %8.1e %6.0d %8.2f\n', iter, y[1:ddnvar]"*ddc[:], DIMACS_error, err1, err2, err3, err4, err5, err6, cg_iter1, cg_iter2, eq_norm, arank, titi) - # @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.0d %9.0d %6.0d\n", iter, dot(y, ctmp'), DIMACS_error, err1, err2, err3, err4, err5, err6, cg_iter1, cg_iter2, cg_iter2) + #@sprintf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.0d %9.0d %8.1e %6.0d %8.2f\n', iter, y[1:ddnvar]"*ddc[:], DIMACS_error, err[1], err[2], err[3], err[4], err[5], err[6], cg_iter1, cg_iter2, eq_norm, arank, titi) + # @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.0d %9.0d %6.0d\n", iter, dot(y, ctmp'), DIMACS_error, err[1], err[2], err[3], err[4], err[5], err[6], cg_iter1, cg_iter2, cg_iter2) if solver.verb > 1 if solver.kit == 0 - @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.2f\n", solver.iter, -dot(solver.y, solver.model.b') + solver.model.b_const, DIMACS_error, solver.err1, solver.err2, solver.err3, solver.err4, solver.err5, solver.err6,solver.itertime) + @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %8.2f\n", solver.iter, dobj, DIMACS_error, solver.err[1], solver.err[2], solver.err[3], solver.err[4], solver.err[5], solver.err[6],solver.itertime) else - @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %7.0d %7.0d %8.2f\n", solver.iter, -dot(solver.y, solver.model.b') + solver.model.b_const, DIMACS_error, solver.err1, solver.err2, solver.err3, solver.err4, solver.err5, solver.err6, solver.cg_iter_pre, solver.cg_iter_cor,solver.itertime) + @printf("%3.0d %16.8e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %7.0d %7.0d %8.2f\n", solver.iter, dobj, DIMACS_error, solver.err[1], solver.err[2], solver.err[3], solver.err[4], solver.err[5], solver.err[6], solver.cg_iter_pre, solver.cg_iter_cor,solver.itertime) end else if solver.kit == 0 - @printf("%3.0d %16.8e %9.2e %8.2f\n", solver.iter, -dot(solver.y, solver.model.b') + solver.model.b_const, DIMACS_error, solver.itertime) + @printf("%3.0d %16.8e %9.2e %8.2f\n", solver.iter, dobj, DIMACS_error, solver.itertime) else - @printf("%3.0d %16.8e %9.2e %9.0d %8.2f\n", solver.iter, -dot(solver.y, solver.model.b') + solver.model.b_const, DIMACS_error, solver.cg_iter_pre + solver.cg_iter_cor, solver.itertime) + @printf("%3.0d %16.8e %9.2e %9.0d %8.2f\n", solver.iter, dobj, DIMACS_error, solver.cg_iter_pre + solver.cg_iter_cor, solver.itertime) end end end if DIMACS_error < solver.eDIMACS solver.status = 1 - solver.y = solver.y if solver.verb > 0 - println("Primal objective: ", -dot(solver.y, solver.model.b') + solver.model.b_const) - if solver.model.nlin > 0 - println("Dual objective: ", -btrace(solver.model.nlmi, solver.model.C, solver.X) - dot(solver.model.d_lin', solver.X_lin)) - else - println("Dual objective: ", -btrace(solver.model.nlmi, solver.model.C, solver.X) ) - end - end + println("Primal objective: ", dobj) + println("Dual objective: ", pobj) + end end - if DIMACS_error > 1e55 + if pobj < -1e55 solver.status = 2 if solver.verb > 0 - println("WARNING: Problem probably infeasible (stopping status = 2)") + @warn("Problem probably infeasible (stopping status = 2)") end - elseif DIMACS_error > 1e55 || abs(dot(solver.y, solver.model.b')) > 1e55 + elseif dobj > 1e55 solver.status = 3 if solver.verb > 0 - println("WARNING: Problem probably unbounded or infeasible (stopping status = 3)") + @warn("Problem probably unbounded (stopping status = 3)") + end + elseif any(isnan, solver.err) + solver.status = 5 + if solver.verb > 0 + @warn("Got NaN (stopping status = 5)") end end @@ -569,47 +603,15 @@ end ```Functions for the iterative solver follow``` -struct MyA{T} - W::Vector{Matrix{T}} - AA::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}} - nlin::Int64 - C_lin::SparseArrays.SparseMatrixCSC{Float64,Int64} - X_lin - S_lin_inv +struct MyA{T,M} + W::LRO.ShapedSolution{T,FactoredMatrix{T}} + model::M to::TimerOutputs.TimerOutput end -function (t::MyA)(Ax::Vector{T}, x::Vector{T}) where {T} +function (t::MyA)(Ax::Vector, x::Vector) @timeit t.to "Ax" begin - nlmi = length(t.AA) - m = size(t.AA[1],1) - ax1 = zeros(m,1) - if nlmi > 0 - for ilmi = 1:nlmi - waxwtmp = Matrix{T}(undef,size(t.W[ilmi])) - waxw = Matrix{T}(undef,size(t.W[ilmi])) - # @timeit t.to "Ax1" begin - ax = Vector{T}(undef,size(t.AA[ilmi],2)) - # end - # @timeit t.to "Ax2" begin - mul!(ax, transpose(t.AA[ilmi]), x) - # ax = transpose(t.AA[ilmi]) * x - # end - # @timeit t.to "Ax3" begin - # waxw .= t.W[ilmi] * mat(ax) * t.W[ilmi] - mul!(waxwtmp,t.W[ilmi], mat(ax)) - mul!(waxw, waxwtmp, t.W[ilmi]) - # end - # @timeit t.to "Ax4" begin - ax1 .+= t.AA[ilmi] * waxw[:] - # end - end - end - if t.nlin>0 - ax1 .+= t.C_lin * ((t.X_lin .* t.S_lin_inv) .* (t.C_lin' * x)) - end - - mul!(Ax,I(m),ax1[:]) + LRO.eval_schur_complement!(t.model, t.W, x, Ax) end end @@ -623,23 +625,24 @@ end function Prec_for_CG_beta(solver,halpha) - nlmi = solver.model.nlmi + nlmi = LRO.num_matrices(solver.model) kk = solver.erank .* ones(Int64,nlmi,1) - nvar = solver.model.n + nvar = solver.model.meta.ncon ntot=0 if nlmi > 0 - for ilmi=1:nlmi - ntot = ntot + size(solver.W[ilmi],1) + for i in LRO.matrix_indices(solver.model) + ntot = ntot + size(solver.W[i],1) end end halpha.AAAATtau = zeros(nvar) if nlmi > 0 - for ilmi = 1:nlmi - n = size(solver.W[ilmi],1); + for i in LRO.matrix_indices(solver.model) + ilmi = i.value + n = size(solver.W[i],1); k = kk[ilmi]; - F = eigen(solver.W[ilmi]); + F = eigen(solver.W[i]); lambdaf = F.values lambda_s = lambdaf[1:n-k] @@ -656,14 +659,18 @@ function Prec_for_CG_beta(solver,halpha) halpha.AAAATtau += ttau^2 .* ZZZ end - if solver.model.nlin > 0 - halpha.AAAATtau .+= diag(solver.model.C_lin * spdiagm((solver.X_lin .* solver.S_lin_inv)[:]) * solver.model.C_lin') + if LRO.num_scalars(solver.model) > 0 + for i in 1:solver.model.meta.ncon + solver.BBBB[i, i] = 0 + LRO.add_schur_complement!(solver.model, solver.X[LRO.ScalarIndex] .* solver.S_lin_inv, ScalarIndex, solver.BBBB) + end + halpha.AAAATtau .+= diag(solver.BBBB) end end end -struct MyM_beta - AA +struct MyM_beta{M} + model::M AAAATtau end @@ -674,36 +681,37 @@ end function Prec_for_CG_tilS_prep(solver::MySolver{T},halpha) where {T} @timeit solver.to "prec" begin - nlmi = solver.model.nlmi - kk = solver.erank .* ones(Int64,nlmi,1) + nlmi = LRO.num_matrices(solver.model) + kk = solver.erank .* ones(Int64,nlmi) # kk[2] = 3 # halpha.Z = SparseMatrixCSC{T}[] halpha.Z = Matrix{T}[] - nvar = solver.model.n + nvar = solver.model.meta.ncon halpha.AAAATtau = spzeros(nvar,nvar) ntot=0 if nlmi > 0 - for ilmi=1:nlmi - ntot = ntot + size(solver.W[ilmi],1) + for i in LRO.matrix_indices(solver.model) + ntot = ntot + size(solver.W[i],1) end end sizeS=0 if nlmi > 0 - for ilmi=1:nlmi - sizeS += kk[ilmi] * size(solver.W[ilmi],1) + for i in LRO.matrix_indices(solver.model) + sizeS += kk[i.value] * size(solver.W[i],1) end end lbt = 1; lbs=1; if nlmi > 0 - for ilmi = 1:nlmi - n = size(solver.W[ilmi],1); - k = kk[ilmi]; + for i in LRO.matrix_indices(solver.model) + ilmi = i.value + n = size(solver.W[i],1) + k = kk[ilmi] - F = eigen(Float64.(solver.W[ilmi])); + F = eigen(Float64.(solver.W[i])) vectf = F.vectors lambdaf = F.values @@ -740,24 +748,21 @@ function Prec_for_CG_tilS_prep(solver::MySolver{T},halpha) where {T} end end - if solver.model.nlin > 0 - halpha.AAAATtau .+= solver.model.C_lin * spdiagm((solver.X_lin .* solver.S_lin_inv)[:]) * solver.model.C_lin' + if LRO.num_scalars(solver.model) > 0 + LRO.add_schur_complement!(solver.model, solver.X[LRO.ScalarIndex] .* solver.S_lin_inv, ScalarIndex, halpha.AAAATtau) end - didi = 0 - for ilmi = 1:nlmi - didi += size(solver.W[ilmi],1) - end k = kk[1] if k > 1 #slow formula # @timeit solver.to "prec3" begin - t = zeros(nvar, k*didi) + t = zeros(nvar, k*ntot) if nlmi > 0 - for ilmi = 1:nlmi - n = size(solver.W[ilmi],1) + for mat_idx in LRO.matrix_indices(solver.model) + ilmi = mat_idx.value + n = size(solver.W[mat_idx],1) k = kk[ilmi] TT = kron(halpha.Umat[ilmi],halpha.Z[ilmi]) - t[1:nvar,lbt:lbt+k*n-1] .= solver.model.AA[ilmi] * TT + t[1:nvar,lbt:lbt+k*n-1] .= jac(solver.model, mat_idx)' * TT lbt = lbt + k*n end end @@ -770,7 +775,7 @@ function Prec_for_CG_tilS_prep(solver::MySolver{T},halpha) where {T} AAAATtau_d = spdiagm(sqrt.(1 ./ diag(halpha.AAAATtau))); # @timeit solver.to "prec3" begin - # t = zeros(nvar, k*didi) + # t = zeros(nvar, k*ntot) # if nlmi > 0 # for ilmi = 1:nlmi # if kk[ilmi] == 0 @@ -797,7 +802,7 @@ function Prec_for_CG_tilS_prep(solver::MySolver{T},halpha) where {T} # # mul!(S,t',t) # end - S, lbt = prec_alpha_S!(solver,halpha,AAAATtau_d,kk,didi,lbt,sizeS) + S, lbt = prec_alpha_S!(solver,halpha,AAAATtau_d,kk,ntot,lbt,sizeS) end # Schur complement for the SMW formula @@ -808,8 +813,8 @@ function Prec_for_CG_tilS_prep(solver::MySolver{T},halpha) where {T} end -struct MyM - AA +struct MyM{M} + model::M AAAATtau Umat Z @@ -819,30 +824,33 @@ end function prec_alpha_S!(solver::MySolver{T},halpha,AAAATtau_d,kk,didi,lbt,sizeS) where {T} @timeit solver.to "prec3" begin S = Matrix{T}(undef,sizeS,sizeS) - nvar = solver.model.n + nvar = solver.model.meta.ncon t = Matrix{T}(undef,nvar,kk[1]*didi) - if solver.model.nlmi > 0 - for ilmi = 1:solver.model.nlmi + if LRO.num_matrices(solver.model) > 0 + for mat_idx in LRO.matrix_indices(solver.model) + ilmi = mat_idx.value if kk[ilmi] == 0 continue end - n = size(solver.W[ilmi],1) + n = size(solver.W[mat_idx],1) k = kk[ilmi] @timeit solver.to "prec30" begin - AAs = AAAATtau_d * solver.model.AA[ilmi] + # We can reuse the buffer for different `i` + # because we directly apply the multiplication + # with `Umat`. + AU = reduce(vcat, [ + (LRO.unsafe_jtprod( + solver.model, + -AAAATtau_d[i,:], + mat_idx, + ) * halpha.Umat[ilmi])' + for i in axes(AAAATtau_d, 1) + ]) end - # @timeit solver.to "prec31" begin - ii_, jj_, aa_ = findnz(AAs) - qq_ = floor.(Int64,(jj_ .- 1) ./ n) .+ 1 - pp_ = mod.(jj_ .- 1, n) .+ 1 - aau = Vector{T}(undef,length(aa_)) - aau .= aa_ .* halpha.Umat[ilmi][qq_] - AU = sparse(ii_,pp_,aau,nvar,n) - # end - if solver.model.nlmi>1 + if LRO.num_matrices(solver.model) > 1 # @timeit solver.to "prec32" begin - didi1 = size(solver.W[ilmi],1) + didi1 = size(solver.W[mat_idx],1) ttmp = Matrix{T}(undef,nvar,kk[ilmi]*didi1) mul!(ttmp, AU, halpha.Z[ilmi]) t[1:nvar,lbt:lbt+k*n-1] = ttmp @@ -866,17 +874,18 @@ end function (t::MyM)(Mx::Vector{T}, x::Vector{T}) where {T} nvar = size(x,1) - nlmi = length(t.AA) + nlmi = LRO.num_matrices(t.model) - yy2 = zeros(nvar,1) + yy2 = zeros(nvar) y33 = zeros(T,0) AAAAinvx = t.AAAATtau\x if nlmi > 0 - for ilmi = 1:nlmi - y22 = t.AA[ilmi]' * AAAAinvx - y33 = [y33; vec(t.Z[ilmi]' * mat(y22) * t.Umat[ilmi])] + for mat_idx in LRO.matrix_indices(t.model) + ilmi = mat_idx.value + y22 = LRO.unsafe_jtprod(t.model, AAAAinvx, mat_idx) + y33 = [y33; vec(t.Z[ilmi]' * y22 * t.Umat[ilmi])] end end @@ -884,7 +893,8 @@ function (t::MyM)(Mx::Vector{T}, x::Vector{T}) where {T} ii = 0 if nlmi > 0 - for ilmi = 1:nlmi + for mat_idx in LRO.matrix_indices(t.model) + ilmi = mat_idx.value n = size(t.Umat[ilmi],1) k = size(t.Umat[ilmi],2) yy = zeros(n*n) @@ -893,7 +903,7 @@ function (t::MyM)(Mx::Vector{T}, x::Vector{T}) where {T} yy .+= kron(t.Umat[ilmi][:,i],xx) ii += n end - yy2 .+= t.AA[ilmi] * yy + LRO.add_jprod!(t.model, reshape(yy, n, n), yy2, mat_idx) end end diff --git a/src/initial_point.jl b/src/initial_point.jl index 2634163..2704e15 100644 --- a/src/initial_point.jl +++ b/src/initial_point.jl @@ -16,66 +16,50 @@ end function find_initial!(solver) - C_lin = solver.model.C_lin' + b2 = 1 .+ abs.(LRO.cons_constant(solver.model)') + n = length(b2) + solver.y .= 0 - n = length(solver.model.b) - solver.y = zeros(n,1) - - b2 = 1 .+ abs.(solver.model.b') - f = zeros(1,n) - for i=1:solver.model.nlmi + f = zeros(n) + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value + dim = LRO.side_dimension(solver.model, mat_idx) if solver.initpoint == 0 Eps = 1.0 else - f = norm(b2)/(1+norm(solver.model.AA[i])) - Eps = sqrt.(solver.model.msizes[i]).* max(1,sqrt.(solver.model.msizes[i]).* f) + f = norm(b2)/(1+LRO.norm_jac(solver.model, mat_idx)) + Eps = sqrt(dim) * max(1, sqrt.(dim) * f) end - solver.X[i] = Eps * Matrix(1.0I, Int64(solver.model.msizes[i]), Int64(solver.model.msizes[i])) + solver.X[mat_idx] .= Eps * Matrix(1.0I, dim, dim) if solver.initpoint == 0 - Eta = solver.model.n + Eta = n else - mf = max(f,norm(solver.model.C[i],2)) - mf = (1 + mf)./ sqrt(solver.model.msizes[i]) - Eta = sqrt(solver.model.msizes[i]).* max(1,mf) + mf = max(f, norm(LRO.grad(solver.model, mat_idx), 2)) + mf = (1 + mf) / dim + Eta = sqrt(dim).* max(1, mf) end - solver.S[i] = Eta * Matrix(1.0I, Int64(solver.model.msizes[i]), Int64(solver.model.msizes[i])) + solver.S[mat_idx] .= Eta * Matrix(1.0I, dim, dim) end - p = zeros(1,n) - pp = zeros(1,n) - dd = size(solver.model.d_lin,1) - if solver.model.nlin>0 - if solver.initpoint == 0 - Epss = 1.0 - else - for j=1:n - normClin = 1+norm(solver.model.C_lin[j,:]) - p[j] = b2[j] ./ normClin; - end - Epss = max(1, maximum(p)) - end - solver.X_lin = 1 .* Epss * ones(dd,1) - - if solver.initpoint == 0 - Etaa = 1.0 - else - for j=1:n - pp[j]=norm(solver.model.C_lin[j,:]) - end - mf = max(maximum(pp),norm(solver.model.d_lin)) - mf = (0 + mf) ./ sqrt(dd) - Etaa = max(1,mf) - end - solver.S_lin = 1 .* Etaa * ones(dd,1) - solver.S_lin_inv = 1 ./ solver.S_lin + if solver.initpoint == 0 + pp = zeros(n) + p = zeros(n) else - solver.X_lin = Float64[]; solver.S_lin = Float64[] + pp = [norm(LRO.jac(solver.model, j, LRO.ScalarIndex)) for j in 1:n] + p = b2 ./ (1 .+ pp) end - if solver.model.nlin==0 - solver.X_lin=Float64[] - solver.S_lin=Float64[] - solver.S_lin_inv=Float64[] + Epss = max(1.0, maximum(p, init = 0.0)) + solver.X[LRO.ScalarIndex] .= Epss + + if solver.initpoint == 0 + Etaa = 1.0 + else + mf = max(maximum(pp, init = 0.0), norm(LRO.grad(solver.model, LRO.ScalarIndex))) + mf = (0 + mf) ./ sqrt(LRO.num_scalars(solver.model)) + Etaa = max(1, mf) end + solver.S[LRO.ScalarIndex] .= Etaa + solver.S_lin_inv = inv.(solver.S[LRO.ScalarIndex]) end diff --git a/src/kron_etc.jl b/src/kron_etc.jl index d3c6b9f..25ad9ad 100644 --- a/src/kron_etc.jl +++ b/src/kron_etc.jl @@ -7,7 +7,7 @@ function my_kron(A::Matrix{T}, B, C) where {T} TMP = Matrix{T}(undef,size(B,1),size(C,1)) mul!(TMP1,C,A') mul!(TMP,B,TMP1) - return vec(TMP) + return TMP end ########################################################################### function mat(vecA) @@ -18,15 +18,6 @@ function mat(vecA) end ########################################################################### -function btrace(nlmi, X, S) - # compute sum of traces of products of block matrices - trXS = 0 - @inbounds for i = 1:nlmi - trXS += dot(X[i], S[i]) - end - return trXS -end - # How to create a vector of sparse matrices: diff --git a/src/makeBBBB.jl b/src/makeBBBB.jl deleted file mode 100644 index 85b4e4e..0000000 --- a/src/makeBBBB.jl +++ /dev/null @@ -1,228 +0,0 @@ -function makeBBBB_rank1(n,nlmi,B,G,to) - @timeit to "BBBB_rank1" begin - tmp = zeros(Float64, n, n) - BBBB = zeros(Float64, n, n) - for ilmi = 1:nlmi - # @timeit to "BBBB_rank1_a" begin - BB = transpose(B[ilmi] * G[ilmi]) - # end - # @timeit to "BBBB_rank1_b" begin - mul!(tmp,BB',BB) - if ilmi == 1 - BBBB = tmp .^ 2 - else - BBBB += tmp .^ 2 - end - # end - end - end - return BBBB -end - -######################### - -function makeBBBBs(n,nlmi,A,AA,W,to,qA,sigmaA) - BBBB = zeros(Float64, n, n) - @inbounds for ilmi = 1:nlmi - Wilmi = W[ilmi] - AAilmi = AA[ilmi] - Ailmi = A[ilmi,:] - @timeit to "BBBBs" begin - BBBB += makeBBBBsi(ilmi,Ailmi,AAilmi,Wilmi,n,to,qA,sigmaA) - end - end - - return BBBB -end - -# Computes `⟨A * W, W * B⟩` for symmetric sparse matrices `A` and `B` -function _dot(A::SparseMatrixCSC, B::SparseMatrixCSC, W::Matrix) - @assert LinearAlgebra.checksquare(W) == LinearAlgebra.checksquare(A) == LinearAlgebra.checksquare(B) - # After these asserts, we know that `A`, `B` and `W` are square and - # have the same sizes so we can safely use `@inbounds` - result = zero(eltype(A)) - @inbounds for i in axes(A, 2) - nzA = nzrange(A, i) - if !isempty(nzA) - for j in axes(B, 2) - nzB = nzrange(B, j) - if !isempty(nzB) - AW = zero(result) - for k in nzA - AW += nonzeros(A)[k] * W[rowvals(A)[k], j] - end - WB = zero(result) - for k in nzB - WB += W[i, rowvals(B)[k]] * nonzeros(B)[k] - end - result += AW * WB - end - end - end - end - return result -end - -##### -function makeBBBBsi(ilmi,Ailmi,AAilmi,Wilmi::Matrix{T},n,to,qA,sigmaA) where {T} - BBBB = zeros(T, n, n) - tmp1 = Matrix{T}(undef,size(Wilmi, 2), size(Ailmi[1], 1)) - # tmp = Matrix{Float64}(undef,size(Wilmi, 1), size(Wilmi, 1)) - tmp2 = Matrix{T}(undef,size(AAilmi, 1), 1) - tmp3 = Vector{Float64}(undef,size(Wilmi, 1)) - ilmi1 = (ilmi-1)*n - - # @timeit to "BBBBsi" begin - - @inbounds for ii = 1:n - # tmp1 = zeros(Float64,size(Wilmi, 2), size(Ailmi[1], 1)) - i = sigmaA[ii,ilmi] - if nnz(Ailmi[i+1]) > 0 - if ii <= qA[1,ilmi] - tmp = zeros(T,size(Wilmi, 2), size(Ailmi[1], 1)) - # if 1==1 - # @show "one" - # @show ii - @timeit to "BBBBone" begin - @timeit to "BBBBone1" begin - mul!(tmp1,Wilmi,Ailmi[i+1]) - end - @timeit to "BBBBone2" begin - # mul!(tmp,tmp1,Wilmi) - tmp = tmp1 * Wilmi - end - @timeit to "BBBBone3" begin - tmp2 = AAilmi * vec(tmp) - # mul!(tmp2,AAilmi,vec(tmp)) - end - @timeit to "BBBBone4" begin - indi = sigmaA[ii:end,ilmi] - BBBB[indi,i] .= -tmp2[indi] - BBBB[i,indi] .= -tmp2[indi] - # @show BBBB[1:2,1:2] - end - end - # elseif ii <= qA[2,ilmi] - elseif 1==0 - # @show "two" - @timeit to "BBBBtwo" begin - mul!(tmp1,Ailmi[i+1],Wilmi) - @inbounds for jj = ii:n - j = sigmaA[jj,ilmi] - Ajjj = Ailmi[j+1] - if !iszero(nnz(Ajjj)) - ttt = 0.0 - # @timeit to "BBBBtwo_i" begin - @inbounds for jjjjAj in axes(Ajjj, 2) - for k in nzrange(Ajjj, jjjjAj) - # @timeit to "BBBBtwo_ii_A" begin - iiijAj = rowvals(Ajjj)[k] - # end - # vvvj = -vvv_j[iAj] - # @timeit to "BBBBtwo_i_B" begin - ttt1 = dot(tmp1[:,iiijAj],Wilmi[:,jjjjAj]) - # end - # @timeit to "BBBBtwo_i_C" begin - ttt += ttt1 * nonzeros(Ajjj)[k] - end - # end - end - # end - BBBB[i,j] = ttt - if !=(i,j) - BBBB[j,i] = ttt - end - end - end - end - # end - else - @timeit to "BBBBthree" begin - # @show "three" - # @show ilmi - # @show ii - if !iszero(nnz(Ailmi[i+1])) - if nnz(Ailmi[i+1]) > 1 - # @timeit to "BBBBthree>1" begin - # @show iii_i,myAiii.jind - # iii_is = iii_i[1:Int64(sqrt(length(iii_i)))] - # jjj_i = myAiii.jind - # vvv_i = myAiii.nzval - @inbounds for jj = ii:n - j = sigmaA[jj,ilmi] - if !iszero(nnz(Ailmi[j+1])) - # @timeit to "BBBBthree_1>1" begin - # iii_js = iii_j[1:Int64(sqrt(length(iii_j)))] - # jjj_j = myAjjj.jind - # vvv_j = myAjjj.nzval - # ttt = 0.0 - # end - # @timeit to "BBBBthree_2>1" begin - ttt = _dot(Ailmi[i+1], Ailmi[j+1], Wilmi) - # end - # @inbounds for iAj in eachindex(iii_j) - # ttt1 = 0.0 - # iiijAj = iii_j[iAj] - # jjjjAj = jjj_j[iAj] - # vvvj = vvv_j[iAj] - # @inbounds for iAi in eachindex(iii_i) - # iiiiAi = iii_i[iAi] - # jjjiAi = jjj_i[iAi] - # vvvi = vvv_i[iAi] - # ttt1 += vvvi * Wilmi[iiiiAi,iiijAj] * Wilmi[jjjiAi,jjjjAj] - # # ttt1 -= vvv_i[iAi] * Wilmi[iii_i[iAi],iiijAj] * Wilmi[jjj_i[iAi],jjjjAj] - # end - # ttt += ttt1 * vvvj - # end - # @timeit to "BBBBthree_3>1" begin - if i >= j - BBBB[i,j] = ttt - else - BBBB[j,i] = ttt - end - # end - end - # end - end - else - @timeit to "BBBBthree=1" begin - # A is symmetric - iiiiAi = jjjiAi = only(rowvals(Ailmi[i+1])) - vvvi = only(nonzeros(Ailmi[i+1])) - @inbounds for jj = ii:n - j = sigmaA[jj,ilmi] - Ajjj = Ailmi[j+1] - # As we sort the matrices in decreasing `nnz` order, - # the rest of matrices is either zero or have only - # one entry - if !iszero(nnz(Ajjj)) - iiijAj = jjjjAj = only(rowvals(Ajjj)) - vvvj = only(nonzeros(Ajjj)) - ttt = vvvi * Wilmi[iiiiAi,iiijAj] * Wilmi[jjjiAi,jjjjAj] * vvvj - if i >= j - BBBB[i,j] = ttt - else - BBBB[j,i] = ttt - end - end - end - end - end - end - end - end - end - end -# end - return BBBB -end - - -function makeRHS(nlmi,AA,W,S,Rp,Rd) - h = Rp # RHS for the Hessian equation - for i = 1:nlmi - # h = h + AA[i] * my_kron(G[i], G[i], (G[i]' * Rd[i] * G[i] + diagm(D[i]))) - h = h + AA[i]*vec(W[i]*(Rd[i]+S[i])*W[i]); #equivalent - end -return h -end diff --git a/src/model.jl b/src/model.jl deleted file mode 100644 index f2ad152..0000000 --- a/src/model.jl +++ /dev/null @@ -1,287 +0,0 @@ -export prepare_model_data, MyModel - -using SparseArrays -using Printf -using TimerOutputs -using LinearAlgebra - -""" - MyModel - -Model representing the problem: -```math -\\begin{aligned} -\\max {} & b^\\top y - b_\\text{const} -\\\\ -& \\sum_{j=1}^n y_j A_{i,j} \\preceq C_i -\\qquad -\\forall i \\in \\{1,\\ldots,\\text{nlmi}\\} -\\\\ -& C_\\text{lin} y \\le d_\\text{lin} -\\end{aligned} -``` -The fields of the `struct` as related to the arrays of the above formulation as follows: - -* The ``i``th PSD constraint is of size `msize[i] × msisze[i]` -* The matrix ``C_i`` is given by `C[i]` which should be equal to `-A[i,1]`. -* The matrix ``A_{i,j}`` is given by `-A[i,j+1]` as well as `myA[(i-1)*n + j]`. -* The vectorization `vec(A[i,j+1])` is also given by `-AA[i][:,j]` -* If `datarank == -1`, ``A_{i,j}`` is also equal to `-B[i][j,:] * B[i][j,:]'`. -* The matrix ``A_{i,j}`` has `nzA[j,i]` nonzero entries -* The index `j = sigmaA[k,i]` is the `k`th matrix ``A_{i,j}`` of the largest number of nonzeros. -* The first `qA[1,i] = qA[2,i]` matrices are considered as dense in the computation. -""" -mutable struct MyModel - A::Matrix{SparseArrays.SparseMatrixCSC{Float64,Int}} - AA::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}} - B::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}} - C::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}} - nzA::Matrix{Int64} - sigmaA::Matrix{Int64} - qA::Matrix{Int64} - b::Vector{Float64} - b_const::Float64 - d_lin::SparseArrays.SparseVector{Float64, Int64} - C_lin::SparseArrays.SparseMatrixCSC{Float64, Int64} - n::Int64 - msizes::Vector{Int64} - nlin::Int64 - nlmi::Int64 - - function MyModel( - A::Matrix{SparseArrays.SparseMatrixCSC{Float64,Int}}, - AA::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}}, - B::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}}, - C::Vector{SparseArrays.SparseMatrixCSC{Float64,Int}}, - nzA::Matrix{Int64}, - sigmaA::Matrix{Int64}, - qA::Matrix{Int64}, - b::Vector{Float64}, - b_const::Float64, - d_lin::SparseArrays.SparseVector{Float64, Int64}, - C_lin::SparseArrays.SparseMatrixCSC{Float64, Int64}, - n::Int64, - msizes::Vector{Int64}, - nlin::Int64, - nlmi::Int64 - ) - - model = new() - model.A = A - model.AA = AA - model.B = B - model.C = C - model.nzA = nzA - model.sigmaA = sigmaA - model.qA = qA - model.b = b - model.b_const = b_const - model.d_lin = d_lin - model.C_lin = C_lin - model.n = n - model.msizes = msizes - model.nlin = nlin - model.nlmi = nlmi - return model - end -end - - -function prepare_model_data(d,drank) - -msizes = Vector{Int64} -n = Int64(get(d, "nvar", 1)); -msizesa = get(d, "msizes", 1) -if length(msizesa) == 1 - msizes = [convert.(Int64,msizesa)] -else - msizes = convert.(Int64,msizesa[:]) -end -nlin = Int64(get(d, "nlin", 1)) -nlmi = Int64(get(d, "nlmi", 1)) -A = get(d, "A", 1); -b = -get(d, "c", 1); -b_const = -get(d, "b_const", 1); - -if nlin > 0 - d_lin = -get(d, "d", 1) - d_lin = d_lin[:] - C_lin = -get(d, "C", 1) -else - d_lin = sparse([0.; 0.]) - C_lin = sparse([0. 0.;0. 0.]) -end - -model = MyModel(A, _prepare_A(A,drank,κ)..., b, b_const, d_lin, C_lin, n, msizes, nlin, nlmi) - -return model -end - -function _prepare_A(A, datarank, κ) - - nlmi = size(A, 1) - n = size(A, 2) - 1 - AA = SparseMatrixCSC{Float64,Int}[] - B = SparseMatrixCSC{Float64,Int}[] - C = SparseMatrixCSC{Float64,Int}[] - nzA = zeros(Int64,n,nlmi) - sigmaA = zeros(Int64,n,nlmi) - qA = zeros(Int64,2,nlmi) - - for i = 1:nlmi - - push!(C, copy(-A[i, 1])) - - Ai = A[i,:] - m = size(Ai,1) - AAA = prep_AA!(Ai,n) - push!(AA, copy(AAA')) - - if datarank == -1 - Btmp = prep_B(A,n,i) - push!(B, Btmp) - end - - prep_sparse!(A,n,m,i,nzA,sigmaA,qA,κ) - - end - - return AA, B, C, nzA, sigmaA, qA -end - - -function prep_sparse!(A,n,m,i,nzA,sigmaA,qA,κ) - # Simplified data sparsity handling - - for j = 1:n - nzA[j,i] = nnz(A[i,j+1]) - end - sigmaA[:,i] = sortperm(nzA[:,i], rev = true) - sisi = nzA[sigmaA[:,i],i] - # @show sisi - - qA[1,i] = n - kappa = κ - for j = 1:n - if sisi[j] <= kappa - qA[1,i] = j-1 - break - end - end - qA[2,i] = qA[1,i] - - # @show qA -end - -function prep_B(A,n,i) - m = size(A[i, 1],1) - Btmp = spzeros(n,m) - - for k = 1:n - ii = rowvals(A[i, k + 1]) - bidx = unique(ii) - if !isempty(bidx) - tmp = Matrix(A[i, k + 1][bidx, bidx]) - # utmp, vtmp = eigen(Hermitian(tmp)) - utmp, vtmp = eigen((tmp + tmp') ./ 2) - bbb = sign.(vtmp[:, end]) .* sqrt.(diag(tmp)) - tmp2 = bbb * bbb' - if norm(tmp - tmp2) > 5.0e-6 - error("Obtained an error of `$(norm(tmp - tmp2)) > 5e-6` when converting matrix into rank `1`, use `datarank = 0` to disable the rank-1 conversion.") - end - Btmp[k, bidx] = bbb - end - end - - return Btmp -end - -function prep_AA!(Ai,n) - - @inbounds Threads.@threads for j = 1:n - if isempty(Ai[j+1]) - Ai[j+1][1, 1] = 0 - end - end - - ntmp = size(Ai[1], 1) * size(Ai[1], 2) - - nnz = 0 - @inbounds for j = 1:n - nnz += SparseArrays.nnz(Ai[j+1]) - end - - iii = zeros(Int64, nnz) - jjj = zeros(Int64, nnz) - vvv = zeros(Float64, nnz) - lb = 1 - @inbounds for j = 1:n - ii,vv = findnz(-(Ai[j+1])[:]) - lf = lb+length(ii)-1 - iii[lb:lf] = ii - jjj[lb:lf] = j .* ones(Int64,length(ii)) - vvv[lb:lf] = float(vv) - lb = lf+1 - end - AAA = sparse(iii,jjj,vvv,ntmp,n) - - return AAA -end - -# end #module - - -# function prep_sparse!(A,n,m,i,nzA,sigmaA,qA) -# This is the Kojima et al data sparsity handling -# d1 = zeros(Float64,n) -# d2 = zeros(Float64,n) -# d3 = zeros(Float64,n) - -# kappa = 100 -# kappa = 500000/m -# for j = 1:n -# nzA[j,i] = nnz(A[i,j+1]) -# end -# sigmaA[:,i] = sortperm(nzA[:,i], rev = true) -# @show nzA[:,i] -# # @show nzA[sigmaA[:,i],i] -# sisi = sort(nzA[sigmaA[:,i],i], rev = true) -# # @show sigmaA[:,i] -# # @show sisi -# cs = cumsum(sisi[end:-1:1]) -# cs = cs[n:-1:1] -# # @show cs - -# for j = 1:n -# d1[j] = kappa * m * nzA[sigmaA[j,i],i] + m^3 + kappa * cs[j] -# d2[j] = kappa * m * nzA[sigmaA[j,i],i] + kappa * (n+1) * cs[j] -# d3[j] = kappa * (2 * kappa * nzA[sigmaA[j,i],i] + 1) * cs[j] -# end - - -# qA[1,i] = 0 -# ktmp = 0 -# for j = 1:n -# if d1[j] > min(d2[j],d3[j]) -# qA[1,i] = j-1 -# ktmp = 1 -# break -# end -# end -# if ktmp == 0 -# qA[1,i] = n -# qA[2,i] = n -# else -# qA[2,i] = 0 -# for j = max(1,qA[1,i]):n -# if d2[j] >= d1[j] || d2[j] > d3[j] -# qA[2,i] = j-1 -# break -# end -# end -# end -# qA[2,i] = max(qA[2,i],qA[1,i]) - -# @show qA - -# end diff --git a/src/predictor_corrector.jl b/src/predictor_corrector.jl index 79b2496..2d9c854 100644 --- a/src/predictor_corrector.jl +++ b/src/predictor_corrector.jl @@ -5,57 +5,56 @@ using GenericLinearAlgebra function predictor(solver::MySolver{T},halpha::Halpha) where {T} solver.predict = true - solver.Rp = solver.model.b + NLPModels.cons!(solver.model, solver.X, solver.Rp) + LinearAlgebra.rmul!(solver.Rp, -1) - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - solver.Rp -= solver.model.AA[i] * solver.X[i][:] - solver.Rd[i] .= solver.model.C[i] - solver.S[i] - mat(solver.model.AA[i]' * solver.y) - solver.Rc[i] .= solver.sigma .* solver.mu .* Matrix(I, length(solver.D[i]), 1) - solver.D[i] .^ 2 - end + for mat_idx = LRO.matrix_indices(solver.model) + i = mat_idx.value + solver.Rd[mat_idx] .= LRO.unsafe_dual_cons(solver.model, solver.y, mat_idx) + solver.Rc[i] .= solver.sigma .* solver.mu .* Matrix(I, length(solver.D[i]), 1) - solver.D[i] .^ 2 end - if solver.model.nlin > 0 - solver.Rp -= solver.model.C_lin * solver.X_lin[:] - solver.Rd_lin = solver.model.d_lin - solver.S_lin - solver.model.C_lin' * solver.y + if LRO.num_scalars(solver.model) > 0 + LRO.dual_cons!(solver.model, solver.y, solver.Rd[LRO.ScalarIndex], LRO.ScalarIndex) end + solver.Rd .-= solver.S if solver.kit == 0 # if direct solver; compute the Hessian matrix - # @timeit solver.to "BBBB" begin - if solver.model.nlmi > 0 - if solver.datarank == -1 - # if 1 == 0 - BBBB = makeBBBB_rank1(solver.model.n, solver.model.nlmi, solver.model.B, solver.G, solver.to) - else - BBBB = makeBBBBs(solver.model.n, solver.model.nlmi, solver.model.A, solver.model.AA, solver.W, solver.to, solver.model.qA, solver.model.sigmaA) - end - else - BBBB = zeros(T, solver.model.n, solver.model.n) - end - if solver.model.nlin > 0 - BBBB .+= solver.model.C_lin * spdiagm((solver.X_lin .* solver.S_lin_inv)[:]) * solver.model.C_lin' - end - BBBB = Hermitian(BBBB, :L) + LRO.schur_complement!(solver.model, solver.W, solver.BBBB) end # end - if solver.model.nlmi > 0 - h = makeRHS(solver.model.nlmi,solver.model.AA,solver.W,solver.S,solver.Rp,solver.Rd) - else - h = copy(solver.Rp) + # RHS for the Hessian equation + tmp = similar(solver.X) + if !isempty(tmp[LRO.ScalarIndex]) + tmp[LRO.ScalarIndex] .= spdiagm(solver.W[LRO.ScalarIndex]) * solver.Rd[LRO.ScalarIndex] + solver.X[LRO.ScalarIndex] end - if solver.model.nlin > 0 - h .+= solver.model.C_lin * (spdiagm((solver.X_lin .* solver.Si_lin)[:]) * solver.Rd_lin + solver.X_lin) + for i in LRO.matrix_indices(solver.model) + tmp[i] .= solver.W[i] * (solver.Rd[i] + solver.S[i]) * solver.W[i] end + h = solver.y_buffer + NLPModels.jprod!(solver.model, solver.X, tmp, h) + h .+= solver.Rp + + # solving the linear system() if solver.kit == 0 # direct solver + BBBB = LinearAlgebra.Hermitian(solver.BBBB) # @timeit solver.to "backslash" begin if ishermitian(BBBB) + if parent(BBBB) isa SparseMatrixCSC + # Convert to dense because + # 1. Cholesky is not implemented for `MultiFloat` for sparse + # 2. It causes issues like https://github.com/JuliaSparse/SparseArrays.jl/issues/630, although that issue could be fixed by densifying the vector `h`. + BBBB = LinearAlgebra.Hermitian(Matrix(parent(BBBB)), LinearAlgebra.sym_uplo(BBBB.uplo)) + end try - cholBBBB1, cholBBBB2 = cholesky(BBBB) - solver.cholBBBB = cholBBBB1 + solver.cholBBBB = cholesky(BBBB).L catch err + if !(err isa LinearAlgebra.PosDefException) + rethrow(err) + end if solver.verb > 0 println("Matrix H not positive definite, trying to regularize") end @@ -63,30 +62,29 @@ function predictor(solver::MySolver{T},halpha::Halpha) where {T} solver.regcount += 1 if solver.regcount > 5 if solver.verb > 0 - println("WARNING: too many regularizations of H, giving up") + @warn("too many regularizations of H, giving up") end solver.cholBBBB = I(size(BBBB, 1)) solver.status = 3 return end - while isposdef(BBBB) == false - BBBB = BBBB + 1e-4 .* I(size(BBBB, 1)) + while !isposdef(BBBB) + solver.BBBB .= solver.BBBB .+ 1e-4 .* I(size(solver.BBBB, 1)) + BBBB = LinearAlgebra.Hermitian(BBBB) icount = icount + 1 if icount > 1000 if solver.verb > 0 - println("WARNING: H cannot be made positive definite, giving up") + @warn("H cannot be made positive definite, giving up") end solver.cholBBBB = I(size(BBBB, 1)) solver.status = 3 return end end - solver.cholBBBB = cholesky(BBBB) - else - solver.cholBBBB = copy(solver.cholBBBB) + solver.cholBBBB = cholesky(BBBB).L end solver.dely = solver.cholBBBB \ h - solver.dely = solver.cholBBBB' \ (solver.cholBBBB \ h) + solver.dely = solver.cholBBBB' \ solver.dely # delyy = solver.dely else @warn("System matrix not Hermitian, stopping Loraine") @@ -115,15 +113,15 @@ function predictor(solver::MySolver{T},halpha::Halpha) where {T} # end else - A = MyA(solver.W,solver.model.AA,solver.model.nlin,solver.model.C_lin,solver.X_lin,solver.S_lin_inv,solver.to) + A = MyA(solver.W, solver.model, solver.to) if solver.preconditioner == 0 M = MyM_no(solver.to) elseif solver.preconditioner == 1 Prec_for_CG_tilS_prep(solver,halpha) - M = MyM(solver.model.AA, halpha.AAAATtau, halpha.Umat, halpha.Z, halpha.cholS) + M = MyM(solver.model, halpha.AAAATtau, halpha.Umat, halpha.Z, halpha.cholS) elseif solver.preconditioner == 2 || solver.preconditioner == 4 Prec_for_CG_beta(solver,halpha) - M = MyM_beta(solver.model.AA, halpha.AAAATtau) + M = MyM_beta(solver.model, halpha.AAAATtau) end # @timeit solver.to "CG predictor" begin @@ -155,40 +153,43 @@ function sigma_update(solver::MySolver{T}) where {T} else expon_used = max(1, min(solver.expon, T(3) * step_pred^2)) end - if btrace(solver.model.nlmi, solver.Xn, solver.Sn) .< 0 + dotXnSn = isempty(solver.Xn) ? zero(T) : dot(solver.Xn, solver.Sn) + if dotXnSn .< 0 solver.sigma = T(0.8) else - if solver.model.nlmi > 0 - tmp1 = btrace(solver.model.nlmi, solver.Xn, solver.Sn) + if LRO.num_scalars(solver.model) > 0 + tmp2 = dot(solver.Xn_lin', solver.Sn_lin) else - tmp1 = 0 + tmp2 = 0 end - if solver.model.nlin > 0 - tmp2 = dot(solver.Xn_lin', solver.Sn_lin) - else - tmp2 = 0 - end - tmp12 = (tmp1 + tmp2) / (sum(solver.model.msizes) + solver.model.nlin) + tmp12 = (dotXnSn + tmp2) / (LRO.num_scalars(solver.model) + sum(Base.Fix1(LRO.side_dimension, solver.model), LRO.matrix_indices(solver.model), init = 0)) tmp12 = convert(Float64, tmp12) mu = Float64(solver.mu) solver.sigma = min(1.0, ((tmp12) / mu) ^ Float64(expon_used)) end return solver.sigma -end +end -function corrector(solver,halpha) +function corrector(solver::MySolver{T},halpha) where {T} solver.predict = false - h = solver.Rp #RHS for the linear system() - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - h += solver.model.AA[i] * my_kron(solver.G[i], solver.G[i], (solver.G[i]' * solver.Rd[i] * solver.G[i] + spdiagm(solver.D[i]) - Diagonal((solver.sigma * solver.mu) ./ solver.D[i]) - solver.RNT[i])) # RHS using my_kron() - end - end - if solver.model.nlin > 0 + X = similar(solver.X) + if LRO.num_scalars(solver.model) > 0 tmp = (solver.delX_lin .* solver.delS_lin) .* (solver.Si_lin) - (solver.sigma * solver.mu) .* (solver.Si_lin) - h = h + solver.model.C_lin * (spdiagm((solver.X_lin .* solver.Si_lin)[:]) * solver.Rd_lin + solver.X_lin + tmp) + X[LRO.ScalarIndex] .= spdiagm((solver.X[LRO.ScalarIndex] .* solver.Si_lin)[:]) * solver.Rd[LRO.ScalarIndex] + solver.X[LRO.ScalarIndex] + tmp end + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value + W = solver.W[mat_idx] + X[mat_idx] .= my_kron( + W.factor, + W.factor, + W.factor' * solver.Rd[mat_idx] * W.factor + spdiagm(solver.D[i]) - Diagonal((solver.sigma * solver.mu) ./ solver.D[i]) - solver.RNT[i], + ) + end + h = solver.y_buffer + NLPModels.jprod!(solver.model, solver.X, X, h) + h .+= solver.Rp # solving the linear system() if solver.kit == 0 # direct solver @@ -221,16 +222,17 @@ function corrector(solver,halpha) # end # end else - A = MyA(solver.W,solver.model.AA,solver.model.nlin,solver.model.C_lin,solver.X_lin,solver.S_lin_inv,solver.to) + A = MyA(solver.W, solver.model, solver.to) if solver.preconditioner == 0 M = MyM_no(solver.to) elseif solver.preconditioner == 1 - M = MyM(solver.model.AA, halpha.AAAATtau, halpha.Umat, halpha.Z, halpha.cholS) + M = MyM(solver.model, halpha.AAAATtau, halpha.Umat, halpha.Z, halpha.cholS) else - M = MyM_beta(solver.model.AA, halpha.AAAATtau) + M = MyM_beta(solver.model, halpha.AAAATtau) end @timeit solver.to "CG corrector" begin + # `maxIter = 10000` fails on 32-bit, we need `maxIter = Int64(10000)` solver.dely, exit_code, num_iters = cg(A, h[:]; tol = Float64(solver.tol_cg), maxIter = Int64(10000), precon = M) end solver.cg_iter_cor += num_iters @@ -245,22 +247,23 @@ function corrector(solver,halpha) end function find_step(solver::MySolver{T}) where {T} - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi + if LRO.num_matrices(solver.model) > 0 + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value @timeit solver.to "find_step_A" begin - solver.delS[i] .= solver.Rd[i] .- mat(solver.model.AA[i]' * solver.dely) - Ξ = my_kron(solver.W[i], solver.W[i], solver.delS[i]) + solver.delS[i] .= solver.Rd[mat_idx] .- LRO.unsafe_jtprod(solver.model, solver.dely, mat_idx) + Ξ = vec(my_kron(solver.W[mat_idx].matrix, solver.W[mat_idx], solver.delS[i])) if solver.predict - solver.delX[i] .= mat(-solver.X[i][:] .- Ξ) + solver.delX[i] .= mat(-solver.X[mat_idx][:] .- Ξ) else - solver.delX[i] .= mat(((solver.sigma * solver.mu) .* solver.Si[i] .- solver.X[i])[:] .- Ξ .+ my_kron(solver.G[i], solver.G[i], solver.RNT[i])) + solver.delX[i] .= mat(((solver.sigma * solver.mu) .* solver.Si[i] .- solver.X[mat_idx])[:] .- Ξ .+ vec(my_kron(solver.W[mat_idx].factor, solver.W[mat_idx].factor, solver.RNT[i]))) end end # determining steplength to stay feasible @timeit solver.to "find_step_B" begin - delSb = solver.G[i]' * solver.delS[i] * solver.G[i] - delXb = solver.Gi[i] * solver.delX[i] * solver.Gi[i]' + delSb = solver.W[mat_idx].factor' * solver.delS[i] * solver.W[mat_idx].factor + delXb = solver.W[mat_idx].factor_inv * solver.delX[i] * solver.W[mat_idx].factor_inv' end @timeit solver.to "find_step_C" begin @@ -291,7 +294,7 @@ function find_step(solver::MySolver{T}) where {T} end end - if solver.model.nlin > 0 + if LRO.num_scalars(solver.model) > 0 find_step_lin(solver) else solver.alpha_lin = 1 @@ -300,25 +303,26 @@ function find_step(solver::MySolver{T}) where {T} if solver.predict # solution update - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - solver.Xn[i] = solver.X[i] + solver.alpha[i] .* solver.delX[i] - solver.Sn[i] = solver.S[i] + solver.beta[i] .* solver.delS[i] - deed = solver.D[i] * ones(1, Int(solver.model.msizes[i])) + ones(Int(solver.model.msizes[i]), 1) * solver.D[i]' - solver.RNT[i] = -(solver.Gi[i] * solver.delX[i] * solver.delS[i] * solver.G[i] + solver.G[i]' * solver.delS[i] * solver.delX[i] * solver.Gi[i]') ./ deed + if LRO.num_matrices(solver.model) > 0 + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value + solver.Xn[i] = solver.X[mat_idx] + solver.alpha[i] .* solver.delX[i] + solver.Sn[i] = solver.S[mat_idx] + solver.beta[i] .* solver.delS[i] + dim = LRO.side_dimension(solver.model, mat_idx) + deed = solver.D[i] * ones(dim)' + ones(LRO.side_dimension(solver.model, mat_idx)) * solver.D[i]' + solver.RNT[i] = -(solver.W[mat_idx].factor_inv * solver.delX[i] * solver.delS[i] * solver.W[mat_idx].factor + solver.W[mat_idx].factor' * solver.delS[i] * solver.delX[i] * solver.W[mat_idx].factor_inv') ./ deed end end else solver.yold = solver.y - solver.y = solver.y + minimum([solver.beta; solver.beta_lin]) * solver.dely - if solver.model.nlmi > 0 - for i = 1:solver.model.nlmi - solver.X[i] = solver.X[i] + minimum([solver.alpha; solver.alpha_lin]) .* solver.delX[i] - solver.X[i] = (solver.X[i] + solver.X[i]') ./ 2 - solver.S[i] = solver.S[i] + minimum([solver.beta; solver.beta_lin]) .* solver.delS[i] - solver.S[i] = (solver.S[i] + solver.S[i]') ./ 2 - end - end + solver.y .+= minimum([solver.beta; solver.beta_lin]) * solver.dely + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value + solver.X[mat_idx] .+= minimum([solver.alpha; solver.alpha_lin]) .* solver.delX[i] + solver.X[mat_idx] .= (solver.X[mat_idx] .+ solver.X[mat_idx]') ./ 2 + solver.S[mat_idx] .+= minimum([solver.beta; solver.beta_lin]) .* solver.delS[i] + solver.S[mat_idx] .= (solver.S[mat_idx] + solver.S[mat_idx]') ./ 2 + end end return @@ -326,19 +330,21 @@ end function find_step_lin(solver) - solver.delS_lin = solver.Rd_lin - solver.model.C_lin' * solver.dely + LRO.jtprod!(solver.model, solver.dely, solver.delS_lin, LRO.ScalarIndex) + solver.delS_lin .*= -1 + solver.delS_lin .+= solver.Rd[LRO.ScalarIndex] if solver.predict - solver.delX_lin = -solver.X_lin - (solver.X_lin) .* (solver.Si_lin) .* solver.delS_lin + solver.delX_lin = -solver.X[LRO.ScalarIndex] - (solver.X[LRO.ScalarIndex]) .* (solver.Si_lin) .* solver.delS_lin else - solver.delX_lin = -solver.X_lin - (solver.X_lin) .* (solver.Si_lin) .* solver.delS_lin + (solver.sigma * solver.mu) .* (solver.Si_lin) + solver.RNT_lin + solver.delX_lin = -solver.X[LRO.ScalarIndex] - (solver.X[LRO.ScalarIndex]) .* (solver.Si_lin) .* solver.delS_lin + (solver.sigma * solver.mu) .* (solver.Si_lin) + solver.RNT_lin end - mimiX_lin = minimum(solver.delX_lin ./ solver.X_lin) + mimiX_lin = minimum(solver.delX_lin ./ solver.X[LRO.ScalarIndex]) if mimiX_lin .> -1e-6 solver.alpha_lin = 0.99 else solver.alpha_lin = min(1, -solver.tau / mimiX_lin) end - mimiS_lin = minimum(solver.delS_lin ./ solver.S_lin) + mimiS_lin = minimum(solver.delS_lin ./ solver.S[LRO.ScalarIndex]) if mimiS_lin .> -1e-6 solver.beta_lin = 0.99 else @@ -347,16 +353,16 @@ function find_step_lin(solver) if solver.predict # solution update - solver.Xn_lin = solver.X_lin + solver.alpha_lin .* solver.delX_lin - solver.Sn_lin = solver.S_lin + solver.beta_lin .* solver.delS_lin + solver.Xn_lin = solver.X[LRO.ScalarIndex] + solver.alpha_lin .* solver.delX_lin + solver.Sn_lin = solver.S[LRO.ScalarIndex] + solver.beta_lin .* solver.delS_lin solver.RNT_lin = -(solver.delX_lin .* solver.delS_lin) .* solver.Si_lin else - # @show solver.X_lin + # @show solver.X[LRO.ScalarIndex] # @show mimiX_lin - solver.X_lin = solver.X_lin + minimum([solver.alpha; solver.alpha_lin]) .* solver.delX_lin - solver.S_lin = solver.S_lin + minimum([solver.beta; solver.beta_lin]) .* solver.delS_lin - solver.S_lin_inv = 1 ./ solver.S_lin + solver.X[LRO.ScalarIndex] .+= minimum([solver.alpha; solver.alpha_lin]) .* solver.delX_lin + solver.S[LRO.ScalarIndex] .+= minimum([solver.beta; solver.beta_lin]) .* solver.delS_lin + solver.S_lin_inv = inv.(solver.S[LRO.ScalarIndex]) end return diff --git a/src/prepare_W.jl b/src/prepare_W.jl index 5f3811b..84ae98b 100644 --- a/src/prepare_W.jl +++ b/src/prepare_W.jl @@ -2,38 +2,40 @@ using TimerOutputs using FameSVD using MultiFloats -function try_cholesky(solver, X, i::Integer, name::String) +function try_cholesky(solver, X, name::String) try - return cholesky(X[i]) + return cholesky(X) catch if solver.verb > 0 println("Matrix $name not positive definite, trying to regularize") end icount = 0 - while isposdef(X[i]) == false - X[i] += 1e-5 .* I(size(X[i], 1)) + while isposdef(X) == false + X .+= 1e-5 .* I(size(X, 1)) icount += 1 if icount > 1000 if solver.verb > 0 - println("WARNING: $name cannot be made positive definite, giving up") + @warn("$name cannot be made positive definite, giving up") end solver.status = 4 - return I(size(X[i], 1)) + return I(size(X, 1)) end end - return cholesky(X[i]) + return cholesky(X) end end function prepare_W(solver::MySolver{T}) where {T} # @timeit solver.to "prpr" begin - for i = 1:solver.model.nlmi + solver.W[LRO.ScalarIndex] .= solver.X[LRO.ScalarIndex] .* solver.S_lin_inv + for mat_idx = LRO.matrix_indices(solver.model) + i = mat_idx.value # @timeit to "prpr1" begin - Ctmp = try_cholesky(solver, solver.X, i, "X") - CtmpS = try_cholesky(solver, solver.S, i, "S") - # Ctmp = cholesky(solver.X[i]) - # CtmpS = cholesky(solver.S[i]) + Ctmp = try_cholesky(solver, solver.X[mat_idx], "X") + CtmpS = try_cholesky(solver, solver.S[mat_idx], "S") + # Ctmp = cholesky(solver.X[mat_idx]) + # CtmpS = cholesky(solver.S[mat_idx]) @timeit solver.to "prep W SVD" begin CCtmp = Matrix{T}(undef,size(CtmpS.L,1),size(CtmpS.L,1)) mul!(CCtmp, (CtmpS.L)' , Ctmp.L) @@ -51,29 +53,30 @@ function prepare_W(solver::MySolver{T}) where {T} Di2 = try Diagonal(1.0 ./ sqrt.(Dtmp)) catch err - println("WARNING: Numerical difficulties, giving up") + @warn("Numerical difficulties, giving up") solver.status = 4 Diagonal(I(size(solver.Dtmp, 1))) end # @timeit to "prpr3a" begin - solver.G[i] = Ctmp.L * V * Di2 + W = solver.W[mat_idx] + W.factor .= Ctmp.L * V * Di2 + W.factor_inv .= inv(W.factor) + LinearAlgebra.mul!(W.matrix, W.factor, W.factor') # end # @timeit to "prpr3" begin - solver.Gi[i] = inv(solver.G[i]) - solver.W[i] = solver.G[i] * solver.G[i]' # end # @timeit to "prpr4" begin # solver.Si[i] = inv(solver.S[i]) solver.Si[i] = (CtmpS.L)' \ ((CtmpS.L) \ (I(size(solver.Si[i],1)))) # S[i] inverse # DDtmp = (CtmpS.U * solver.G[i]) # DDtmp = DDtmp' * DDtmp - DDtmp = solver.G[i]' * solver.S[i] * solver.G[i] + DDtmp = solver.W[mat_idx].factor' * solver.S[mat_idx] * solver.W[mat_idx].factor DDtmp = (DDtmp + DDtmp') ./ 2.0 try solver.DDsi[i] = (1.0 ./ sqrt.(diag(DDtmp,0))) catch err - println("WARNING: Numerical difficulties, giving up") + @warn("Numerical difficulties, giving up") solver.DDsi[i] = diag(I(size(DDtmp, 1))) solver.status = 4 return @@ -82,13 +85,13 @@ function prepare_W(solver::MySolver{T}) where {T} end # end end - if solver.model.nlin > 0 - solver.Si_lin = 1.0 ./ solver.S_lin + if LRO.num_scalars(solver.model) > 0 + solver.Si_lin = inv.(solver.S[LRO.ScalarIndex]) else solver.Si_lin = [] end # end - return solver.D, solver.G, solver.Gi, solver.W, solver.Si, solver.DDsi, solver.Si_lin + return solver.D, solver.W, solver.Si, solver.DDsi, solver.Si_lin end diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index f55fca0..c61c96a 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -35,23 +35,12 @@ function tests() exclude = [ # No constraints r"test_solve_TerminationStatus_DUAL_INFEASIBLE$", - - # Unable to bridge RotatedSecondOrderCone to PSD because the dimension is too small: got 2, expected >= 3 - r"test_conic_SecondOrderCone_INFEASIBLE$", - r"test_constraint_PrimalStart_DualStart_SecondOrderCone$", - - # MathOptInterface.ITERATION_LIMIT == MathOptInterface.DUAL_INFEASIBLE - r"test_linear_DUAL_INFEASIBLE_2$", - r"test_conic_SecondOrderCone_no_initial_bound$", - r"test_conic_SecondOrderCone_negative_post_bound_2$", - r"test_conic_SecondOrderCone_negative_post_bound_3$", - - # PosDefException: matrix is not positive definite; Cholesky factorization failed. r"test_attribute_SolveTimeSec$", r"test_attribute_RawStatusString$", r"test_objective_ObjectiveFunction_blank$", + + # Warning: too many regularizations of H, giving up r"test_linear_transform$", - r"test_linear_DUAL_INFEASIBLE$", ], ) return From 57a45acf04e9827bd78aa177a81358203b5929b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 30 Jun 2026 12:40:05 +0200 Subject: [PATCH 2/8] Add support for SolveTimeSec --- src/Solvers.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Solvers.jl b/src/Solvers.jl index 4029e89..d29b8f4 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -240,6 +240,7 @@ function SolverCore.solve!( solve(solver.solver, solver.halpha) solver.stats.status = STATUS_MAP[solver.solver.status + 1] solver.stats.objective = NLPModels.obj(solver.solver.model, solver.solver.X) + solver.stats.elapsed_time = solver.solver.tottime return end From 10f20456f36f4782e3c2bdc2e5d3cae3e7d453af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 1 Jul 2026 17:24:11 +0200 Subject: [PATCH 3/8] Add support for BarrierIteration --- src/Solvers.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Solvers.jl b/src/Solvers.jl index d29b8f4..bef12ef 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -241,6 +241,7 @@ function SolverCore.solve!( solver.stats.status = STATUS_MAP[solver.solver.status + 1] solver.stats.objective = NLPModels.obj(solver.solver.model, solver.solver.X) solver.stats.elapsed_time = solver.solver.tottime + solver.stats.iter = solver.solver.iter return end From 9297a34c9122416888de02485a6fad982dd97173 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 1 Jul 2026 18:20:39 +0200 Subject: [PATCH 4/8] Fix initial point --- src/initial_point.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/initial_point.jl b/src/initial_point.jl index 2704e15..6e4d90b 100644 --- a/src/initial_point.jl +++ b/src/initial_point.jl @@ -36,7 +36,7 @@ function find_initial!(solver) Eta = n else mf = max(f, norm(LRO.grad(solver.model, mat_idx), 2)) - mf = (1 + mf) / dim + mf = (1 + mf) / sqrt(dim) Eta = sqrt(dim).* max(1, mf) end solver.S[mat_idx] .= Eta * Matrix(1.0I, dim, dim) From b01094dbf9b57e0b5b932fbee83b97014ba00b2a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 1 Jul 2026 21:18:14 +0200 Subject: [PATCH 5/8] Improve bench --- examples/benchmark_sdplib.jl | 45 +++++++++++++----------------------- examples/merge_results.jl | 20 ++++++++++++---- 2 files changed, 31 insertions(+), 34 deletions(-) diff --git a/examples/benchmark_sdplib.jl b/examples/benchmark_sdplib.jl index 32a2fdb..1a2e49a 100644 --- a/examples/benchmark_sdplib.jl +++ b/examples/benchmark_sdplib.jl @@ -21,10 +21,7 @@ # # Only options common to both `main` and `nlpmodel` are used, so the same # script produces comparable CSVs on either branch. The CSV is written with a -# plain `print` (no CSV.jl dependency). -# -# The active project must have `Chairmarks` (used on the worker for a robust -# minimum solve time) in addition to `JuMP` and `Loraine`. +# plain `print` (no CSV.jl dependency), and it needs only `JuMP` and `Loraine`. using Distributed import Printf @@ -63,14 +60,9 @@ end function start_worker() pid = addprocs(1; exeflags = `--project=$(Base.active_project())`)[1] ospid = remotecall_fetch(getpid, pid) - # Load in a separate eval from the function definition below: `@b` must be - # macro-expanded *after* `using Chairmarks` has taken effect. remotecall_wait(Core.eval, pid, Main, quote using JuMP import Loraine - using Chairmarks - end) - remotecall_wait(Core.eval, pid, Main, quote function solve_problem(path, kit) model = read_from_file(path) set_optimizer(model, Loraine.Optimizer{Float64}) @@ -80,17 +72,10 @@ function start_worker() set_attribute(model, "maxit", 100) set_attribute(model, "datasparsity", 8) optimize!(model) # warm up: absorb this problem's compilation - t = @elapsed optimize!(model) - # Millisecond solves are noise-dominated (GC/scheduler jitter), so - # take a robust minimum over repeated warm re-solves. Multi-second - # solves are measured once: they aren't noisy, and re-running them - # would risk blowing past the per-problem timeout. - if t < 0.5 - t = (@b optimize!(model) seconds = 1).time - end + wall = @elapsed optimize!(model) return ( - t, - solve_time(model), # solver's own elapsed time (no harness) + wall, + solve_time(model), # solver's own `SolveTimeSec` barrier_iterations(model), objective_value(model), string(termination_status(model)), @@ -121,9 +106,8 @@ function solve_capped(pid, path, kit, timeout) return (:timeout, nothing) end -function bench(; - out = get(ARGS, 1, "sdplib_results.csv"), - label = get(ARGS, 2, "loraine"), +function bench(label; + out = label * ".csv", maxn = parse(Int, get(ENV, "SDPLIB_MAXN", "1000")), maxm = parse(Int, get(ENV, "SDPLIB_MAXM", "3000")), kit = parse(Int, get(ENV, "SDPLIB_KIT", "0")), @@ -151,11 +135,16 @@ function bench(; warmup = joinpath(DATA, problems[1] * ".dat-s") fetch(launch(pid, warmup, kit)) - open(out, "w") do io - println( - io, - "label,problem,m,n,status,time_s,solve_time_s,iterations,objective,optimal,rel_gap", - ) + # Append so re-running accumulates more rows (one solve = one row); the + # minimum over rows is taken later, in `merge_results.jl`. + write_header = !isfile(out) || filesize(out) == 0 + open(out, "a") do io + if write_header + println( + io, + "label,problem,m,n,status,time_s,solve_time_s,iterations,objective,optimal,rel_gap", + ) + end ntot = length(problems) for (idx, name) in enumerate(problems) r = ref[name] @@ -215,5 +204,3 @@ function bench(; rmprocs(pid) println("\nWrote ", out) end - -bench() diff --git a/examples/merge_results.jl b/examples/merge_results.jl index 0843ff2..660b656 100644 --- a/examples/merge_results.jl +++ b/examples/merge_results.jl @@ -14,9 +14,18 @@ using CSV using DataFrames import Printf +# Load a CSV that may hold several rows per problem (from repeated/appended +# runs) and collapse each problem to its best run — the row with the smallest +# finite `solve_time_s` (falling back to the first row if none finished). function load(path) df = CSV.read(path, DataFrame) - return df, isempty(df.label) ? basename(path) : first(df.label) + label = isempty(df.label) ? basename(path) : first(df.label) + best = combine(groupby(df, :problem)) do sub + fin = findall(x -> !ismissing(x) && isfinite(x), sub.solve_time_s) + i = isempty(fin) ? 1 : fin[argmin(sub.solve_time_s[fin])] + sub[i, :] + end + return best, label end function fmt_time(t) @@ -43,14 +52,15 @@ function merge(baseline = ARGS[1], comparison = ARGS[2]) base, base_label = load(baseline) comp, comp_label = load(comparison) - # Compare on `time_s`: the robust (Chairmarks minimum) warm-solve time, - # free of compilation and one-shot noise. + # Compare on `solve_time_s`: the solver's own `SolveTimeSec` (minimum over + # the repeated warm re-solves), free of the worker harness. `time_s` (wall) + # is kept in the CSV as a cross-check. a = select( base, :problem, :n, :m, - :time_s => :t_base, + :solve_time_s => :t_base, :iterations => :it_base, :objective => :obj_base, :status => :st_base, @@ -58,7 +68,7 @@ function merge(baseline = ARGS[1], comparison = ARGS[2]) b = select( comp, :problem, - :time_s => :t_comp, + :solve_time_s => :t_comp, :iterations => :it_comp, :objective => :obj_comp, :status => :st_comp, From 580abcd722a1fade5b1c0d03da68244b3ae3afcf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 2 Jul 2026 08:41:45 +0200 Subject: [PATCH 6/8] All runs on files --- examples/benchmark_sdplib.jl | 100 +++++++++++++++++++++-------------- 1 file changed, 61 insertions(+), 39 deletions(-) diff --git a/examples/benchmark_sdplib.jl b/examples/benchmark_sdplib.jl index 1a2e49a..ef471dd 100644 --- a/examples/benchmark_sdplib.jl +++ b/examples/benchmark_sdplib.jl @@ -72,14 +72,24 @@ function start_worker() set_attribute(model, "maxit", 100) set_attribute(model, "datasparsity", 8) optimize!(model) # warm up: absorb this problem's compilation - wall = @elapsed optimize!(model) - return ( - wall, + sample() = ( solve_time(model), # solver's own `SolveTimeSec` barrier_iterations(model), objective_value(model), string(termination_status(model)), ) + # Fast problems are re-solved several times (one row per solve, so + # `merge` can take the minimum); each row records its own wall time + # and `SolveTimeSec`. `reps` targets ~2s of solving, capped at 30, + # so slow problems are solved once. + wall = @elapsed optimize!(model) + samples = [(wall, sample()...)] + reps = clamp(round(Int, 2.0 / wall), 1, 30) + for _ in 2:reps + w = @elapsed optimize!(model) + push!(samples, (w, sample()...)) + end + return samples end end) return pid, ospid @@ -106,15 +116,21 @@ function solve_capped(pid, path, kit, timeout) return (:timeout, nothing) end -function bench(label; +function bench( + label, + only = get(ENV, "SDPLIB_ONLY", nothing); out = label * ".csv", maxn = parse(Int, get(ENV, "SDPLIB_MAXN", "1000")), maxm = parse(Int, get(ENV, "SDPLIB_MAXM", "3000")), kit = parse(Int, get(ENV, "SDPLIB_KIT", "0")), timeout = parse(Float64, get(ENV, "SDPLIB_TIMEOUT", "60")), ) - only = haskey(ENV, "SDPLIB_ONLY") ? Set(split(ENV["SDPLIB_ONLY"], ",")) : - nothing + # `only`: `nothing` (all problems within the size caps), a single name or a + # comma-separated string (e.g. "truss1" or "truss1,qap8"), or any iterable + # of names. + only = + isnothing(only) ? nothing : + only isa AbstractString ? Set(split(only, ",")) : Set(string.(only)) ref = read_reference() problems = String[] @@ -148,57 +164,63 @@ function bench(label; ntot = length(problems) for (idx, name) in enumerate(problems) r = ref[name] - status, t, st, iters, obj = "ok", NaN, NaN, -1, NaN outcome, payload = solve_capped(pid, joinpath(DATA, name * ".dat-s"), kit, timeout) - if outcome === :ok - t, st, iters, obj, status = payload + # `samples`: one `(wall, solve_time, iters, obj, status)` per solve. + # Error/timeout collapse to a single synthetic sample. + samples = if outcome === :ok + payload elseif outcome === :error - status = "ERROR: " * sprint(showerror, payload)[1:min(end, 60)] + msg = "ERROR: " * sprint(showerror, payload)[1:min(end, 60)] + [(NaN, NaN, -1, NaN, msg)] else # Stuck: SIGKILL the worker's OS process, then respawn a fresh # one (and re-warm it) for the remaining problems. run(`kill -9 $ospid`) rmprocs(pid; waitfor = 0) - status, t = "TIMEOUT", timeout pid, ospid = start_worker() fetch(launch(pid, warmup, kit)) + [(timeout, NaN, -1, NaN, "TIMEOUT")] end - gap = isfinite(obj) && isfinite(r.opt) ? - abs(obj - r.opt) / max(1, abs(r.opt)) : NaN - println( - io, - join( - [ - label, - name, - r.m, - r.n, - "\"$status\"", - round(t, sigdigits = 5), - round(st, sigdigits = 5), - iters, - obj, - r.opt, - gap, - ], - ",", - ), - ) + for (t, st, iters, obj, status) in samples + gap = isfinite(obj) && isfinite(r.opt) ? + abs(obj - r.opt) / max(1, abs(r.opt)) : NaN + println( + io, + join( + [ + label, + name, + r.m, + r.n, + "\"$status\"", + round(t, sigdigits = 5), + round(st, sigdigits = 5), + iters, + obj, + r.opt, + gap, + ], + ",", + ), + ) + end + flush(io) + # console: minimum solve time over the samples, and the count + fin = filter(isfinite, [s[2] for s in samples]) Printf.@printf( - "[%2d/%2d] %-12s n=%-5d m=%-5d %8.3fs (solve %8.3fs, %3d it) obj=%-14.6g %s\n", + "[%2d/%2d] %-12s n=%-5d m=%-5d min solve %8.3fs (%2d runs, %3d it) obj=%-14.6g %s\n", idx, ntot, name, r.n, r.m, - t, - st, - iters, - obj, - status, + isempty(fin) ? NaN : minimum(fin), + length(samples), + samples[end][3], + samples[end][4], + samples[end][5], ) - flush(io) end end rmprocs(pid) From 04a65885e42c2c2ed3287949a72c115f3643be86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 2 Jul 2026 16:18:04 +0200 Subject: [PATCH 7/8] fix timeout --- examples/benchmark_sdplib.jl | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/examples/benchmark_sdplib.jl b/examples/benchmark_sdplib.jl index ef471dd..8c3531f 100644 --- a/examples/benchmark_sdplib.jl +++ b/examples/benchmark_sdplib.jl @@ -63,7 +63,8 @@ function start_worker() remotecall_wait(Core.eval, pid, Main, quote using JuMP import Loraine - function solve_problem(path, kit) + function solve_problem(path, kit, timeout) + t0 = time() model = read_from_file(path) set_optimizer(model, Loraine.Optimizer{Float64}) set_attribute(model, "kit", kit) @@ -78,14 +79,17 @@ function start_worker() objective_value(model), string(termination_status(model)), ) - # Fast problems are re-solved several times (one row per solve, so - # `merge` can take the minimum); each row records its own wall time - # and `SolveTimeSec`. `reps` targets ~2s of solving, capped at 30, - # so slow problems are solved once. + # Re-solve for more samples (one row per solve, so `merge` can take + # the minimum), but only while another solve fits the budget (85% of + # the timeout). This keeps the whole call under the per-problem + # timeout while letting a larger `timeout` buy more samples; a solve + # slower than ~half the timeout just yields a single sample. wall = @elapsed optimize!(model) samples = [(wall, sample()...)] - reps = clamp(round(Int, 2.0 / wall), 1, 30) - for _ in 2:reps + # The 32-sample cap only binds for fast problems (slow ones are + # limited by the budget, so a larger `timeout` buys them more rows). + budget = 0.85 * timeout + while length(samples) < 32 && (time() - t0) + wall < budget w = @elapsed optimize!(model) push!(samples, (w, sample()...)) end @@ -95,8 +99,8 @@ function start_worker() return pid, ospid end -launch(pid, path, kit) = - remotecall((p, k) -> Main.solve_problem(p, k), pid, path, kit) +launch(pid, path, kit, timeout) = + remotecall((p, k, t) -> Main.solve_problem(p, k, t), pid, path, kit, timeout) # Solve with a wall-clock cap. `isready` on a *remote* Future blocks while the # worker is busy (it has to query the worker), so we instead let an async task @@ -105,7 +109,7 @@ launch(pid, path, kit) = function solve_capped(pid, path, kit, timeout) ch = Channel{Any}(1) @async put!(ch, try - fetch(launch(pid, path, kit)) + fetch(launch(pid, path, kit, timeout)) catch err err end) @@ -149,7 +153,7 @@ function bench( pid, ospid = start_worker() # Warm up compilation on the smallest problem so timings exclude it. warmup = joinpath(DATA, problems[1] * ".dat-s") - fetch(launch(pid, warmup, kit)) + fetch(launch(pid, warmup, kit, timeout)) # Append so re-running accumulates more rows (one solve = one row); the # minimum over rows is taken later, in `merge_results.jl`. @@ -179,7 +183,7 @@ function bench( run(`kill -9 $ospid`) rmprocs(pid; waitfor = 0) pid, ospid = start_worker() - fetch(launch(pid, warmup, kit)) + fetch(launch(pid, warmup, kit, timeout)) [(timeout, NaN, -1, NaN, "TIMEOUT")] end for (t, st, iters, obj, status) in samples From 92e9f06d0ab2595fddfb38ceb15446d88e9d1c82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 27 Jul 2026 17:39:27 +0200 Subject: [PATCH 8/8] Decrease allocs --- src/Solvers.jl | 61 ++++++++----- src/kron_etc.jl | 8 ++ src/predictor_corrector.jl | 178 +++++++++++++++++++++++-------------- test/Project.toml | 1 + 4 files changed, 161 insertions(+), 87 deletions(-) diff --git a/src/Solvers.jl b/src/Solvers.jl index bef12ef..286563f 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -89,7 +89,8 @@ mutable struct MySolver{T,M} iter::Int64 DIMACS_error::T BBBB::Matrix{T} - cholBBBB + chol_work::Matrix{T} # dense scratch that `cholesky!` factorizes in place + cholBBBB::Cholesky{T,Matrix{T}} status::Int @@ -100,12 +101,12 @@ mutable struct MySolver{T,M} X::LRO.VectorizedSolution{T} S::LRO.VectorizedSolution{T} y::Vector{T} - yold - delX - delS - dely - Xn - Sn + yold::Vector{T} + delX::Vector{Matrix{T}} + delS::Vector{Matrix{T}} + dely::Vector{T} + Xn::Vector{Matrix{T}} + Sn::Vector{Matrix{T}} Si_lin::Vector{T} S_lin_inv::Vector{T} @@ -114,31 +115,37 @@ mutable struct MySolver{T,M} Xn_lin::Vector{T} Sn_lin::Vector{T} - D + D::Vector{Vector{T}} W::LRO.ShapedSolution{T,FactoredMatrix{T}} - Si - DDsi + Si::Vector{Matrix{T}} + DDsi::Vector{Vector{T}} Rp::Vector{T} Rd::LRO.VectorizedSolution{T} - Rc + Rc::Vector{Matrix{T}} - cg_iter_pre - cg_iter_cor - cg_iter_tot + cg_iter_pre::Int64 + cg_iter_cor::Int64 + cg_iter_tot::Int64 - alpha - beta - alpha_lin - beta_lin + alpha::Vector{T} + beta::Vector{T} + alpha_lin::T + beta_lin::T itertime - tottime + tottime::Float64 - RNT - RNT_lin + RNT::Vector{Matrix{T}} + RNT_lin::Vector{T} y_buffer::Vector{T} # Bufer of the same size as `y` + # Per-block dim×dim scratch reused by `my_kron!` and the matrix products in + # `predictor`/`corrector`/`find_step`, so those iterations don't allocate. + kron_tmp::Vector{Matrix{T}} + kron_tmp2::Vector{Matrix{T}} + kron_tmp3::Vector{Matrix{T}} + sol_buffer::LRO.VectorizedSolution{T} # reused RHS-of-Newton-system buffer function MySolver{T}( kit::Int64, @@ -430,6 +437,7 @@ function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} solver.Rd = similar(solver.X) solver.Rp = zeros(T, solver.model.meta.ncon) solver.y_buffer = similar(solver.Rp) + solver.dely = similar(solver.Rp) solver.delX = Matrix{T}[] solver.delS = Matrix{T}[] @@ -454,8 +462,13 @@ function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} solver.Sn = Matrix{T}[] solver.RNT = Matrix{T}[] + solver.kron_tmp = Matrix{T}[] + solver.kron_tmp2 = Matrix{T}[] + solver.kron_tmp3 = Matrix{T}[] + solver.sol_buffer = similar(solver.X) + solver.regcount = 0 - + solver.delS_lin = zeros(LRO.num_scalars(solver.model)) for mat_idx in LRO.matrix_indices(solver.model) dim = LRO.side_dimension(solver.model, mat_idx) @@ -468,6 +481,9 @@ function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} push!(solver.Xn,zeros(dim, dim)) push!(solver.Sn,zeros(dim, dim)) push!(solver.RNT,zeros(dim, dim)) + push!(solver.kron_tmp,zeros(dim, dim)) + push!(solver.kron_tmp2,zeros(dim, dim)) + push!(solver.kron_tmp3,zeros(dim, dim)) end halpha.Umat = Matrix{T}[] @@ -500,6 +516,7 @@ function setup_solver(solver::MySolver{T},halpha::Halpha) where {T} if solver.kit == 0 # if direct solver; compute the Hessian matrix solver.BBBB = zeros(T, ncon, ncon) + solver.chol_work = zeros(T, ncon, ncon) end end diff --git a/src/kron_etc.jl b/src/kron_etc.jl index 25ad9ad..32621fd 100644 --- a/src/kron_etc.jl +++ b/src/kron_etc.jl @@ -9,6 +9,14 @@ function my_kron(A::Matrix{T}, B, C) where {T} mul!(TMP,B,TMP1) return TMP end + +# In-place `my_kron`: writes `B * C * A'` into `dest`, using `tmp` (same size as +# `C * A'`) as scratch. No allocation. +function my_kron!(dest, A, B, C, tmp) + mul!(tmp, C, A') # tmp = C * A' + mul!(dest, B, tmp) # dest = B * (C * A') + return dest +end ########################################################################### function mat(vecA) n = isqrt(length(vecA)) diff --git a/src/predictor_corrector.jl b/src/predictor_corrector.jl index 2d9c854..de20428 100644 --- a/src/predictor_corrector.jl +++ b/src/predictor_corrector.jl @@ -2,6 +2,17 @@ using ConjugateGradients using GenericLinearAlgebra +# Fill `chol_work` with the identity and store its (trivial) Cholesky, so that +# `cholBBBB` stays a valid `Cholesky` even on the give-up path (status = 3). +function set_identity_chol!(solver::MySolver) + fill!(solver.chol_work, 0) + for i in axes(solver.chol_work, 1) + @inbounds solver.chol_work[i, i] = 1 + end + solver.cholBBBB = cholesky!(Hermitian(solver.chol_work)) + return +end + function predictor(solver::MySolver{T},halpha::Halpha) where {T} solver.predict = true @@ -25,12 +36,18 @@ function predictor(solver::MySolver{T},halpha::Halpha) where {T} # end # RHS for the Hessian equation - tmp = similar(solver.X) + tmp = solver.sol_buffer if !isempty(tmp[LRO.ScalarIndex]) tmp[LRO.ScalarIndex] .= spdiagm(solver.W[LRO.ScalarIndex]) * solver.Rd[LRO.ScalarIndex] + solver.X[LRO.ScalarIndex] end - for i in LRO.matrix_indices(solver.model) - tmp[i] .= solver.W[i] * (solver.Rd[i] + solver.S[i]) * solver.W[i] + for mat_idx in LRO.matrix_indices(solver.model) + i = mat_idx.value + # tmp[i] = W * (Rd[i] + S[i]) * W (W is a FactoredMatrix, in place) + RS = solver.kron_tmp[i] + RS .= solver.Rd[mat_idx] .+ solver.S[mat_idx] + mul!(solver.kron_tmp2[i], solver.W[mat_idx], RS) + mul!(solver.kron_tmp3[i], solver.kron_tmp2[i], solver.W[mat_idx]) + tmp[mat_idx] .= solver.kron_tmp3[i] end h = solver.y_buffer NLPModels.jprod!(solver.model, solver.X, tmp, h) @@ -42,56 +59,50 @@ function predictor(solver::MySolver{T},halpha::Halpha) where {T} if solver.kit == 0 # direct solver BBBB = LinearAlgebra.Hermitian(solver.BBBB) # @timeit solver.to "backslash" begin - if ishermitian(BBBB) - if parent(BBBB) isa SparseMatrixCSC - # Convert to dense because - # 1. Cholesky is not implemented for `MultiFloat` for sparse - # 2. It causes issues like https://github.com/JuliaSparse/SparseArrays.jl/issues/630, although that issue could be fixed by densifying the vector `h`. - BBBB = LinearAlgebra.Hermitian(Matrix(parent(BBBB)), LinearAlgebra.sym_uplo(BBBB.uplo)) + if parent(BBBB) isa SparseMatrixCSC + # Convert to dense because + # 1. Cholesky is not implemented for `MultiFloat` for sparse + # 2. It causes issues like https://github.com/JuliaSparse/SparseArrays.jl/issues/630, although that issue could be fixed by densifying the vector `h`. + BBBB = LinearAlgebra.Hermitian(Matrix(parent(BBBB)), LinearAlgebra.sym_uplo(BBBB.uplo)) + end + # Factorize into the reusable `chol_work` buffer (no per-iteration alloc). + # `check = false` returns an unsuccessful factorization instead of + # throwing, so we avoid building/catching a `PosDefException`. + uplo = LinearAlgebra.sym_uplo(BBBB.uplo) + copyto!(solver.chol_work, parent(BBBB)) + chol = cholesky!(Hermitian(solver.chol_work, uplo); check = false) + if !issuccess(chol) + if solver.verb > 0 + println("Matrix H not positive definite, trying to regularize") end - try - solver.cholBBBB = cholesky(BBBB).L - catch err - if !(err isa LinearAlgebra.PosDefException) - rethrow(err) - end + icount = 0 + solver.regcount += 1 + if solver.regcount > 5 if solver.verb > 0 - println("Matrix H not positive definite, trying to regularize") + @warn("too many regularizations of H, giving up") end - icount = 0 - solver.regcount += 1 - if solver.regcount > 5 + set_identity_chol!(solver) + solver.status = 3 + return + end + while true + solver.BBBB .= solver.BBBB .+ 1e-4 .* I(size(solver.BBBB, 1)) + copyto!(solver.chol_work, solver.BBBB) + chol = cholesky!(Hermitian(solver.chol_work, uplo); check = false) + issuccess(chol) && break + icount += 1 + if icount > 1000 if solver.verb > 0 - @warn("too many regularizations of H, giving up") + @warn("H cannot be made positive definite, giving up") end - solver.cholBBBB = I(size(BBBB, 1)) + set_identity_chol!(solver) solver.status = 3 return end - while !isposdef(BBBB) - solver.BBBB .= solver.BBBB .+ 1e-4 .* I(size(solver.BBBB, 1)) - BBBB = LinearAlgebra.Hermitian(BBBB) - icount = icount + 1 - if icount > 1000 - if solver.verb > 0 - @warn("H cannot be made positive definite, giving up") - end - solver.cholBBBB = I(size(BBBB, 1)) - solver.status = 3 - return - end - end - solver.cholBBBB = cholesky(BBBB).L end - solver.dely = solver.cholBBBB \ h - solver.dely = solver.cholBBBB' \ solver.dely - # delyy = solver.dely - else - @warn("System matrix not Hermitian, stopping Loraine") - solver.maxit = 1e10 - solver.status = 2 - solver.cholBBBB = 0 end + solver.cholBBBB = chol + ldiv!(solver.dely, solver.cholBBBB, h) # # Iterative refinement # resid = h - BBBB * solver.dely; # # @show norm(resid - (h[solver.cholBBBB.p] - solver.cholBBBB.L * solver.cholBBBB.U * solver.dely)) @@ -173,7 +184,7 @@ end function corrector(solver::MySolver{T},halpha) where {T} solver.predict = false - X = similar(solver.X) + X = solver.sol_buffer if LRO.num_scalars(solver.model) > 0 tmp = (solver.delX_lin .* solver.delS_lin) .* (solver.Si_lin) - (solver.sigma * solver.mu) .* (solver.Si_lin) X[LRO.ScalarIndex] .= spdiagm((solver.X[LRO.ScalarIndex] .* solver.Si_lin)[:]) * solver.Rd[LRO.ScalarIndex] + solver.X[LRO.ScalarIndex] + tmp @@ -181,11 +192,19 @@ function corrector(solver::MySolver{T},halpha) where {T} for mat_idx in LRO.matrix_indices(solver.model) i = mat_idx.value W = solver.W[mat_idx] - X[mat_idx] .= my_kron( - W.factor, - W.factor, - W.factor' * solver.Rd[mat_idx] * W.factor + spdiagm(solver.D[i]) - Diagonal((solver.sigma * solver.mu) ./ solver.D[i]) - solver.RNT[i], - ) + # C = W.factor' * Rd * W.factor + diag(D) - diag((σμ)/D) - RNT (in place) + C = solver.kron_tmp2[i] + mul!(solver.kron_tmp[i], W.factor', solver.Rd[mat_idx]) + mul!(C, solver.kron_tmp[i], W.factor) + C .-= solver.RNT[i] + d = solver.D[i] + sm = solver.sigma * solver.mu + @inbounds for k in eachindex(d) + C[k, k] += d[k] - sm / d[k] + end + # X[mat_idx] = W.factor * C * W.factor' + my_kron!(solver.kron_tmp3[i], W.factor, W.factor, C, solver.kron_tmp[i]) + X[mat_idx] .= solver.kron_tmp3[i] end h = solver.y_buffer NLPModels.jprod!(solver.model, solver.X, X, h) @@ -195,8 +214,7 @@ function corrector(solver::MySolver{T},halpha) where {T} if solver.kit == 0 # direct solver # @timeit to "corrector backsl" begin # solver.cholBBBB = cholesky(BBBB) - # solver.dely = solver.cholBBBB \ h - solver.dely = solver.cholBBBB' \ (solver.cholBBBB \ h) + ldiv!(solver.dely, solver.cholBBBB, h) # # Iterative refinement # # resid = h - BBBB * solver.dely; # resid = h - solver.cholBBBB * solver.cholBBBB' * solver.dely @@ -250,28 +268,43 @@ function find_step(solver::MySolver{T}) where {T} if LRO.num_matrices(solver.model) > 0 for mat_idx in LRO.matrix_indices(solver.model) i = mat_idx.value + W = solver.W[mat_idx] @timeit solver.to "find_step_A" begin solver.delS[i] .= solver.Rd[mat_idx] .- LRO.unsafe_jtprod(solver.model, solver.dely, mat_idx) - Ξ = vec(my_kron(solver.W[mat_idx].matrix, solver.W[mat_idx], solver.delS[i])) + # Ξ = W * delS * W.matrix' (into kron_tmp[i], scratch kron_tmp2[i]) + Ξ = solver.kron_tmp[i] + my_kron!(Ξ, W.matrix, W, solver.delS[i], solver.kron_tmp2[i]) if solver.predict - solver.delX[i] .= mat(-solver.X[mat_idx][:] .- Ξ) + # delX = symmetrize(-X - Ξ) + Ξ .= .-solver.X[mat_idx] .- Ξ + solver.delX[i] .= (Ξ .+ Ξ') ./ 2 else - solver.delX[i] .= mat(((solver.sigma * solver.mu) .* solver.Si[i] .- solver.X[mat_idx])[:] .- Ξ .+ vec(my_kron(solver.W[mat_idx].factor, solver.W[mat_idx].factor, solver.RNT[i]))) + # delX = symmetrize((σμ) Si - X - Ξ + W.factor RNT W.factor') + K2 = solver.kron_tmp3[i] + my_kron!(K2, W.factor, W.factor, solver.RNT[i], solver.kron_tmp2[i]) + sm = solver.sigma * solver.mu + Ξ .= sm .* solver.Si[i] .- solver.X[mat_idx] .- Ξ .+ K2 + solver.delX[i] .= (Ξ .+ Ξ') ./ 2 end end # determining steplength to stay feasible @timeit solver.to "find_step_B" begin - delSb = solver.W[mat_idx].factor' * solver.delS[i] * solver.W[mat_idx].factor - delXb = solver.W[mat_idx].factor_inv * solver.delX[i] * solver.W[mat_idx].factor_inv' + delSb = solver.kron_tmp2[i] + mul!(solver.kron_tmp[i], W.factor', solver.delS[i]) + mul!(delSb, solver.kron_tmp[i], W.factor) + delXb = solver.kron_tmp3[i] + mul!(solver.kron_tmp[i], W.factor_inv, solver.delX[i]) + mul!(delXb, solver.kron_tmp[i], W.factor_inv') end @timeit solver.to "find_step_C" begin - XXX = solver.DDsi[i]' .* delXb .* solver.DDsi[i] + XXX = solver.kron_tmp[i] + XXX .= solver.DDsi[i]' .* delXb .* solver.DDsi[i] XXX .= (XXX .+ XXX') ./ 2 end @timeit solver.to "find_step_D" begin - mimiX = eigmin(T.(XXX)) + mimiX = eigmin(XXX) end if mimiX .> -1e-6 solver.alpha[i] = 0.99 @@ -280,11 +313,12 @@ function find_step(solver::MySolver{T}) where {T} end @timeit solver.to "find_step_C" begin - XXX = solver.DDsi[i]' .* delSb .* solver.DDsi[i] + XXX = solver.kron_tmp[i] + XXX .= solver.DDsi[i]' .* delSb .* solver.DDsi[i] XXX .= (XXX .+ XXX') ./ 2 end @timeit solver.to "find_step_D" begin - mimiS = eigmin(T.(XXX)) + mimiS = eigmin(XXX) end if mimiS .> -1e-6 solver.beta[i] = 0.99 @@ -306,11 +340,25 @@ function find_step(solver::MySolver{T}) where {T} if LRO.num_matrices(solver.model) > 0 for mat_idx in LRO.matrix_indices(solver.model) i = mat_idx.value - solver.Xn[i] = solver.X[mat_idx] + solver.alpha[i] .* solver.delX[i] - solver.Sn[i] = solver.S[mat_idx] + solver.beta[i] .* solver.delS[i] - dim = LRO.side_dimension(solver.model, mat_idx) - deed = solver.D[i] * ones(dim)' + ones(LRO.side_dimension(solver.model, mat_idx)) * solver.D[i]' - solver.RNT[i] = -(solver.W[mat_idx].factor_inv * solver.delX[i] * solver.delS[i] * solver.W[mat_idx].factor + solver.W[mat_idx].factor' * solver.delS[i] * solver.delX[i] * solver.W[mat_idx].factor_inv') ./ deed + W = solver.W[mat_idx] + solver.Xn[i] .= solver.X[mat_idx] .+ solver.alpha[i] .* solver.delX[i] + solver.Sn[i] .= solver.S[mat_idx] .+ solver.beta[i] .* solver.delS[i] + # R = W.factor_inv delX delS W.factor + W.factor' delS delX W.factor_inv' + A = solver.kron_tmp[i] + B = solver.kron_tmp2[i] + R = solver.kron_tmp3[i] + mul!(A, W.factor_inv, solver.delX[i]) + mul!(B, A, solver.delS[i]) + mul!(R, B, W.factor) + mul!(A, W.factor', solver.delS[i]) + mul!(B, A, solver.delX[i]) + mul!(R, B, W.factor_inv', 1, 1) + # RNT[i] = -R ./ deed with deed[k,l] = D[i][k] + D[i][l] + d = solver.D[i] + RNTi = solver.RNT[i] + @inbounds for l in axes(R, 2), k in axes(R, 1) + RNTi[k, l] = -R[k, l] / (d[k] + d[l]) + end end end else diff --git a/test/Project.toml b/test/Project.toml index 9501b0b..ba6d3f6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Dualization = "191a621a-6537-11e9-281d-650236a99e60" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Loraine = "df0521af-66a0-42b3-822e-f30debbfe642" +LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"