From aa09d8091c1ee710ccdb6f1062aaa3351ece8d17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Wed, 14 Jan 2026 20:05:45 +0100 Subject: [PATCH 1/9] Introduce AVI workspace to allow for hot starts --- interfaces/daqp-julia/src/avi.jl | 155 ++++++++++++++++++++----------- 1 file changed, 101 insertions(+), 54 deletions(-) diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index 9b562ed..d21fab8 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -1,3 +1,101 @@ +Base.@kwdef mutable struct AVISettings + iter_limit::Int = 1000 + inner_iter_limit::Int = 10000 + primal_tol::Float64 = 1e-6 + dual_tol::Float64 = 1e-12 + alpha::Float64 = 1.0 + beta::Float64 = 0.25 +end + +struct AVIWorkspace + H::Matrix{Float64} + f::Vector{Float64} + A::Matrix{Float64} + bu::Vector{Float64} + bl::Vector{Float64} + + H1pI::Matrix{Float64} + H2mI::Matrix{Float64} + H2pI::LU{Float64} # LU + H2::Matrix{Float64} + H2xpf::Vector{Float64} + x::Vector{Float64} + y::Vector{Float64} + daqp_workspace::DAQPBase.Model + kkt_buffer::Vector{Float64} + + settings::AVISettings +end + + +function setup_avi(H,f,A,bu,bl;x0 = zeros(0), daqp_workspace=nothing, settings=AVISettings()) + n,m = length(f),length(bu) + bl = isempty(bl) ? fill(-1e30,m) : bl + + H1 = (H+H') / 4 + H2 = H1+(H-H') / 2 + H1pI,H2mI,H2pI = H1+I, H2-I, lu!(H2+I) + + x = isempty(x0) ? zeros(n) : copy(x0) + y,H2xpf = zeros(n),zeros(n) + + kkt_buffer = Float64[] + sizehint!(kkt_buffer,4*(n+1)*(n+1)) + + daqp_workspace = isnothing(daqp_workspace) ? Model() : daqp_workspace; + DAQPBase.settings(daqp_workspace, Dict(:iter_limit => settings.inner_iter_limit, + :primal_tol => settings.primal_tol, + :dual_tol => settings.dual_tol)) + exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl) + + return exitflag, AVIWorkspace(H,f,A,bu,bl, + H1pI,H2mI,H2pI,H2,H2xpf, + x,y,daqp_workspace,kkt_buffer, + settings) +end + +function solve(ws::AVIWorkspace) + exitflag = -4 + λstar, ASstar = zeros(0), Int[] + n = length(ws.f) + ϵp, ϵd = ws.settings.primal_tol, ws.settings.dual_tol + α,β = ws.settings.alpha, ws.settings.beta + tot_iter,outer_iter = 0,0 + @inbounds for k in 1:ws.settings.iter_limit + ws.H2xpf .= ws.f + mul!(ws.H2xpf, ws.H2mI, ws.x, 1.0, 1.0) + DAQPBase.update(ws.daqp_workspace, nothing, ws.H2xpf, nothing, nothing,nothing) + ws.y[:], _, exitflag, info = DAQPBase.solve(ws.daqp_workspace) + exitflag < 0 && break + tot_iter += info.iterations + + # Same AS -> Check if KKT conditions are satisfied + if info.iterations == 1 + ASu,ASl = _get_AS(info.λ,ϵd) + z,AS = _solve_kkt(ws.H,ws.f,ws.A,ws.bu,ws.bl,ASu,ASl,ws.kkt_buffer) + xt = @view z[1:n] + λ = @view z[n+1:end] + nu = length(ASu) + if all(λ[i] > -ϵd for i in 1:nu) && all(λ[i] < ϵd for i in nu+1:length(λ)) + Ax = ws.A*xt + if all(Ax-ws.bu .<= ϵp) && all(Ax-ws.bl .>=-ϵp) + ws.x .= xt + λstar,ASstar = copy(λ),AS + exitflag, outer_iter = 1,k + break + end + end + axpby!(β,xt,1.0-β,ws.x) + else + axpby!(1.0-α,ws.x,α,ws.y) + mul!(ws.y, ws.H2, ws.x, 1.0, 1.0) + ldiv!(ws.x, ws.H2pI, ws.y) + end + end + info = (status=flag2status[exitflag], AS=ASstar, outer_iterations=outer_iter, iterations=tot_iter) + return ws.x,λstar,exitflag,info +end + function _get_AS(λ,dual_tol) ASu,ASl = Int[],Int[] for i in 1:length(λ) @@ -47,58 +145,7 @@ x ∈ {x: blower ≤ Ax ≤ bupper} such that: * `info` - tuple containing extra information from the solver such as status, active set and number of iterations. """ -function solve_avi(H::AbstractMatrix,f::AbstractVector,A::AbstractMatrix,bupper::AbstractVector,blower::AbstractVector=Float64[]; - alpha=1.0,beta=0.25,primal_tol=1e-6, dual_tol = 1e-12, x0 = zeros(0), iter_limit=1000, inner_iter_limit=10000) - - # Solves VI(Hx +h, C) where C={x| blower <= Ax <= bupper} using Douglas-Rachford - n,m = length(f),length(bupper) - blower = isempty(blower) ? fill(-1e30,m) : blower - - H1 = (H+H') / 4 - H2 = H1+(H-H') / 2 - H1pI,H2mI,H2pI = H1+I, H2-I, lu!(H2+I) - - x = isempty(x0) ? zeros(n) : copy(x0) - y,H2xpf = zeros(n),zeros(n),zeros(n) - - backstep = DAQPBase.Model() - DAQPBase.settings(backstep, Dict(:iter_limit =>inner_iter_limit, :primal_tol=>primal_tol, :dual_tol=>dual_tol)) - exitflag,tsetup = DAQPBase.setup(backstep, H1pI, H2xpf, A, bupper, blower) - exitflag < 0 && return nothing,nothing, (status=:Infeasible, exitflag=exitflag, outer_iterations=0, iterations=0) - - tot_iter = 0 - - kkt_buffer = Float64[] - sizehint!(kkt_buffer,4*(n+1)*(n+1)) - - @inbounds for k in 1:iter_limit - H2xpf .= f - mul!(H2xpf, H2mI, x, 1.0, 1.0) - DAQPBase.update(backstep, nothing, H2xpf, nothing, nothing,nothing) - - y[:], _, exitflag, info = DAQPBase.solve(backstep) - exitflag < 0 && return nothing,nothing, (status=:Infeasible, exitflag=exitflag, outer_iterations=k, iterations=tot_iter) - tot_iter += info.iterations - - # Same AS -> Check if KKT conditions are satisfied - if info.iterations == 1 - ASu,ASl = _get_AS(info.λ,dual_tol) - z,AS = _solve_kkt(H,f,A,bupper,blower,ASu,ASl,kkt_buffer) - xt = @view z[1:n] - λ = @view z[n+1:end] - if all(λ[i] > -dual_tol for i in 1:length(ASu)) && all(λ[i] < dual_tol for i in length(ASu)+1:length(λ)) - Ax = A*xt - if all(Ax-bupper .<= primal_tol) && all(Ax-blower .>=-primal_tol) - return copy(xt),copy(λ),(AS=AS, status=:Solved, outer_iterations=k,iterations=tot_iter) - end - end - axpby!(beta,xt,1.0-beta,x) - else - axpby!(1.0-alpha,x,alpha,y) - mul!(y, H2, x, 1.0, 1.0) - ldiv!(x, H2pI, y) - end - end - return x,nothing, (status=:MaximumIterations, outer_iterations=iter_limit, iterations=tot_iter) +function solve_avi(H::AbstractMatrix,f::AbstractVector,A::AbstractMatrix,bupper::AbstractVector,blower::AbstractVector=Float64[]; x0 = zeros(0), settings=AVISettings()) + exitflag, ws = setup_avi(H,f,A,bupper,blower;x0=zeros(0),settings) + return solve(ws) end - From ef675749d9f8134f0d7e991379ea19e6ca96ed09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= <55484604+darnstrom@users.noreply.github.com> Date: Fri, 26 Dec 2025 14:55:33 +0100 Subject: [PATCH 2/9] DAQPBase.jl v0.3.2 (#103) * Avoid allocations in DAQPBase.solve * Add option to reset C Workspace from Julia * Bump to 0.3.2 --- interfaces/daqp-julia/Project.toml | 2 +- interfaces/daqp-julia/src/api.jl | 42 ++++++++++++++---------------- 2 files changed, 20 insertions(+), 24 deletions(-) diff --git a/interfaces/daqp-julia/Project.toml b/interfaces/daqp-julia/Project.toml index 02fe803..dff76b1 100644 --- a/interfaces/daqp-julia/Project.toml +++ b/interfaces/daqp-julia/Project.toml @@ -1,7 +1,7 @@ name = "DAQPBase" uuid = "833476c3-a8c0-4073-9b64-6473509843fe" authors = ["Daniel Arnström "] -version = "0.3.1" +version = "0.3.2" [deps] DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d" diff --git a/interfaces/daqp-julia/src/api.jl b/interfaces/daqp-julia/src/api.jl index 9f6ca9d..ff83096 100644 --- a/interfaces/daqp-julia/src/api.jl +++ b/interfaces/daqp-julia/src/api.jl @@ -130,6 +130,8 @@ mutable struct Model qpc::QPc qpc_ptr::Ptr{DAQPBase.QPc} has_model::Bool + x::Vector{Float64} + λ::Vector{Float64} function Model() # Setup initial model work = Libc.calloc(1,sizeof(DAQPBase.Workspace)) @@ -176,14 +178,8 @@ function setup(daqp::DAQPBase.Model, qp::DAQPBase.QPj) settings(daqp,old_settings) else daqp.has_model = true - # Quick fix to initialize x for proximal - # will be fixed in DAQP_jll 0.3.2 - if((isempty(qp.H)&&!isempty(qp.f)) || old_settings.eps_prox != 0) - workspace = unsafe_load(daqp.work); - xinit = zeros(workspace.n) - unsafe_copyto!(workspace.x,pointer(xinit),workspace.n); - end - # should be unsafe_copy_to + daqp.x = Vector{Float64}(undef, qp.n) + daqp.λ = Vector{Float64}(undef, qp.m) end return exitflag, setup_time @@ -197,21 +193,19 @@ end function solve(daqp::DAQPBase.Model) if(!daqp.has_model) return zeros(0), NaN, -10, [] end - xstar = zeros(Float64,daqp.qpc.n); - lam = zeros(Float64,daqp.qpc.m); - result= Ref(DAQPResult(xstar,lam)); + result= Ref(DAQPResult(daqp.x,daqp.λ)); exitflag=ccall((:daqp_solve, DAQPBase.libdaqp), Cint, (Ref{DAQPBase.DAQPResult},Ref{DAQPBase.Workspace}), result,daqp.work) - info = (x = xstar, λ=lam, fval=result[].fval, + info = (x = daqp.x, λ=daqp.λ, fval=result[].fval, exitflag=result[].exitflag, status = DAQPBase.flag2status[result[].exitflag], solve_time = result[].solve_time, setup_time = result[].setup_time, iterations= result[].iter, nodes = result[].nodes) - return xstar,result[].fval,result[].exitflag,info + return copy(daqp.x),result[].fval,result[].exitflag,info end @@ -287,6 +281,15 @@ function update(daqp::DAQPBase.Model, H,f,A,bupper,blower,sense=nothing,break_po update_mask, daqp.work,Ref(daqp.qpc)); end +function reset(p::Ptr{DAQPBase.Workspace}) + ccall((:deactivate_constraints,libdaqp),Cvoid,(Ptr{DAQPBase.Workspace},),p); + ccall((:reset_daqp_workspace,libdaqp),Cvoid,(Ptr{DAQPBase.Workspace},),p); +end + +function reset(d::DAQPBase.Model) + reset(d.work) +end + using Downloads function codegen(d::DAQPBase.Model; fname="daqp_workspace", dir="codegen", src=false) @assert(d.has_model, "setup the model before code generation") @@ -294,9 +297,7 @@ function codegen(d::DAQPBase.Model; fname="daqp_workspace", dir="codegen", src=f dir[end] != '/' && (dir*="/") ## Make sure it is correct directory path isdir(dir) || mkdir(dir) - # Make sure workspace is cleared - ccall((:deactivate_constraints,libdaqp),Cvoid,(Ptr{DAQPBase.Workspace},),d.work); - ccall((:reset_daqp_workspace,libdaqp),Cvoid,(Ptr{DAQPBase.Workspace},),d.work); + reset(d) # Make sure workspace is cleared exitflag = ccall((:render_daqp_workspace, libdaqp),Cvoid, (Ptr{DAQPBase.Workspace},Cstring,Cstring,), d.work,fname,dir); @@ -370,13 +371,8 @@ function isfeasible(p::Ptr{DAQPBase.Workspace}, m=nothing, ms=nothing ;validate= @warn "Couldn't validate infeas. with Frakas (err:$(err), fval=$(daqp_ws.fval))" end end - - # Reset the workspace - ccall((:deactivate_constraints,DAQPBase.libdaqp),Cvoid,(Ptr{Cvoid},),p); - ccall((:reset_daqp_workspace,DAQPBase.libdaqp),Cvoid,(Ptr{Cvoid},),p); - #if(exitflag != 1 && exitflag != -1) - # @warn "exitflag: $exitflag" - #end + # Make sure workspace is clean for next solve + reset(p) return exitflag == 1 end From 6e82101fcef79c7d2b0576e3c4aedc5f62048257 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= <55484604+darnstrom@users.noreply.github.com> Date: Fri, 26 Dec 2025 15:29:46 +0100 Subject: [PATCH 3/9] Correct formatting of constraints in README --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 2aaa008..989318d 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ DAQP is a dual active-set solver that solves convex quadratic programs of the fo minimize 0.5 x' H x + f' x subject to l <= x <= u - bl <= Ax <= bu. + bl <= Ax <= bu. ``` Binary constraints of the form $A x \in \lbrace b_l, b_u \rbrace$ are also supported, allowing for mixed-integer quadratic programs to be solved. From da85bb4b52b17475d30215dd779a423a2aa33630 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= <55484604+darnstrom@users.noreply.github.com> Date: Sun, 4 Jan 2026 00:18:36 +0100 Subject: [PATCH 4/9] Make sure new daqp_solve works correctly for Julia <= 1.10 (#104) * Make sure new daqp_solve works correctly for Julia <= 1.10 --- interfaces/daqp-julia/Project.toml | 2 +- interfaces/daqp-julia/src/api.jl | 24 +++++++++++++----------- interfaces/daqp-julia/src/avi.jl | 2 +- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/interfaces/daqp-julia/Project.toml b/interfaces/daqp-julia/Project.toml index dff76b1..1bffee6 100644 --- a/interfaces/daqp-julia/Project.toml +++ b/interfaces/daqp-julia/Project.toml @@ -1,7 +1,7 @@ name = "DAQPBase" uuid = "833476c3-a8c0-4073-9b64-6473509843fe" authors = ["Daniel Arnström "] -version = "0.3.2" +version = "0.3.3" [deps] DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d" diff --git a/interfaces/daqp-julia/src/api.jl b/interfaces/daqp-julia/src/api.jl index ff83096..504f2f9 100644 --- a/interfaces/daqp-julia/src/api.jl +++ b/interfaces/daqp-julia/src/api.jl @@ -195,17 +195,19 @@ function solve(daqp::DAQPBase.Model) if(!daqp.has_model) return zeros(0), NaN, -10, [] end result= Ref(DAQPResult(daqp.x,daqp.λ)); - exitflag=ccall((:daqp_solve, DAQPBase.libdaqp), Cint, - (Ref{DAQPBase.DAQPResult},Ref{DAQPBase.Workspace}), - result,daqp.work) - - info = (x = daqp.x, λ=daqp.λ, fval=result[].fval, - exitflag=result[].exitflag, - status = DAQPBase.flag2status[result[].exitflag], - solve_time = result[].solve_time, - setup_time = result[].setup_time, - iterations= result[].iter, nodes = result[].nodes) - return copy(daqp.x),result[].fval,result[].exitflag,info + GC.@preserve result begin + ccall((:daqp_solve, DAQPBase.libdaqp), Nothing, + (Ref{DAQPBase.DAQPResult},Ref{DAQPBase.Workspace}), + result,daqp.work) + + info = (x = daqp.x, λ=daqp.λ, fval=result[].fval, + exitflag=result[].exitflag, + status = DAQPBase.flag2status[result[].exitflag], + solve_time = result[].solve_time, + setup_time = result[].setup_time, + iterations= result[].iter, nodes = result[].nodes) + return copy(daqp.x),result[].fval,result[].exitflag,info + end end diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index d21fab8..df74638 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -113,7 +113,7 @@ function _solve_kkt(H,f,A,bu,bl,ASu,ASl,kkt_buffer) n = length(f) nkkt = n+length(AS) resize!(kkt_buffer,nkkt^2) - K = reshape(kkt_buffer,nkkt,nkkt) + K = reshape(view(kkt_buffer,1:nkkt^2),(nkkt,nkkt)) @views K[1:n,1:n] = H @views K[n+1:end,1:n] = A[AS,:] @views K[1:n,n+1:end] = A[AS,:]' From c869db6de26609630cd3f5165db6b153aa8d6047 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= <55484604+darnstrom@users.noreply.github.com> Date: Mon, 5 Jan 2026 02:01:44 +0100 Subject: [PATCH 5/9] Avoid compiler to optimize away DAQPResult struct (#106) --- interfaces/daqp-julia/Project.toml | 2 +- interfaces/daqp-julia/src/api.jl | 29 ++++++++++++++--------------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/interfaces/daqp-julia/Project.toml b/interfaces/daqp-julia/Project.toml index 1bffee6..fc79260 100644 --- a/interfaces/daqp-julia/Project.toml +++ b/interfaces/daqp-julia/Project.toml @@ -1,7 +1,7 @@ name = "DAQPBase" uuid = "833476c3-a8c0-4073-9b64-6473509843fe" authors = ["Daniel Arnström "] -version = "0.3.3" +version = "0.3.4" [deps] DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d" diff --git a/interfaces/daqp-julia/src/api.jl b/interfaces/daqp-julia/src/api.jl index 504f2f9..d783114 100644 --- a/interfaces/daqp-julia/src/api.jl +++ b/interfaces/daqp-julia/src/api.jl @@ -193,21 +193,20 @@ end function solve(daqp::DAQPBase.Model) if(!daqp.has_model) return zeros(0), NaN, -10, [] end - result= Ref(DAQPResult(daqp.x,daqp.λ)); - - GC.@preserve result begin - ccall((:daqp_solve, DAQPBase.libdaqp), Nothing, - (Ref{DAQPBase.DAQPResult},Ref{DAQPBase.Workspace}), - result,daqp.work) - - info = (x = daqp.x, λ=daqp.λ, fval=result[].fval, - exitflag=result[].exitflag, - status = DAQPBase.flag2status[result[].exitflag], - solve_time = result[].solve_time, - setup_time = result[].setup_time, - iterations= result[].iter, nodes = result[].nodes) - return copy(daqp.x),result[].fval,result[].exitflag,info - end + result_ptr= Ref(DAQPResult(daqp.x,daqp.λ)); + + ccall((:daqp_solve, DAQPBase.libdaqp), Nothing, + (Ref{DAQPBase.DAQPResult},Ref{DAQPBase.Workspace}), + result_ptr,daqp.work) + + result = unsafe_load(Base.unsafe_convert(Ptr{DAQPResult}, result_ptr)) + info = (x = daqp.x, λ=daqp.λ, fval=result.fval, + exitflag=result.exitflag, + status = DAQPBase.flag2status[result.exitflag], + solve_time = result.solve_time, + setup_time = result.setup_time, + iterations= result.iter, nodes = result.nodes) + return copy(daqp.x),result.fval,result.exitflag,info end From 87d09230c29a652f609ea06801d05a2232edfaa4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Wed, 14 Jan 2026 20:19:17 +0100 Subject: [PATCH 6/9] Handle simple bounds --- interfaces/daqp-julia/src/avi.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index df74638..73383d0 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -48,6 +48,10 @@ function setup_avi(H,f,A,bu,bl;x0 = zeros(0), daqp_workspace=nothing, settings=A :dual_tol => settings.dual_tol)) exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl) + if m > size(A,1) + A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds + end + return exitflag, AVIWorkspace(H,f,A,bu,bl, H1pI,H2mI,H2pI,H2,H2xpf, x,y,daqp_workspace,kkt_buffer, From 2424f44c6e64e5e37b93affa92a1ed1964974aaa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Thu, 15 Jan 2026 19:26:26 +0100 Subject: [PATCH 7/9] Include sense, solve kkt using Schur, start from unconstrained minimum --- interfaces/daqp-julia/src/avi.jl | 78 +++++++++++++++++++++++--------- 1 file changed, 57 insertions(+), 21 deletions(-) diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index 73383d0..37656a8 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -14,6 +14,9 @@ struct AVIWorkspace bu::Vector{Float64} bl::Vector{Float64} + H_factor::LU{Float64} + x_unc::Vector{Float64} + H1pI::Matrix{Float64} H2mI::Matrix{Float64} H2pI::LU{Float64} # LU @@ -27,16 +30,25 @@ struct AVIWorkspace settings::AVISettings end +AVIWorkspace() = AVIWorkspace(zeros(0,0),zeros(0),zeros(0,0),zeros(0),zeros(0), + lu(zeros(0,0)),zeros(0), + zeros(0,0),zeros(0,0),lu(zeros(0,0)),zeros(0,0), zeros(0), + zeros(0),zeros(0), Model(), zeros(0), AVISettings()) + -function setup_avi(H,f,A,bu,bl;x0 = zeros(0), daqp_workspace=nothing, settings=AVISettings()) +function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, settings=AVISettings()) n,m = length(f),length(bu) bl = isempty(bl) ? fill(-1e30,m) : bl + sense = isempty(sense) ? zeros(Cint,m) : sense H1 = (H+H') / 4 H2 = H1+(H-H') / 2 H1pI,H2mI,H2pI = H1+I, H2-I, lu!(H2+I) - x = isempty(x0) ? zeros(n) : copy(x0) + H_factor = lu(H) + x_unc = -(H_factor\f) + + x = isempty(x0) ? copy(x_unc) : copy(x0) y,H2xpf = zeros(n),zeros(n) kkt_buffer = Float64[] @@ -46,25 +58,26 @@ function setup_avi(H,f,A,bu,bl;x0 = zeros(0), daqp_workspace=nothing, settings=A DAQPBase.settings(daqp_workspace, Dict(:iter_limit => settings.inner_iter_limit, :primal_tol => settings.primal_tol, :dual_tol => settings.dual_tol)) - exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl) + exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl, sense) if m > size(A,1) A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds end return exitflag, AVIWorkspace(H,f,A,bu,bl, + H_factor,x_unc, H1pI,H2mI,H2pI,H2,H2xpf, x,y,daqp_workspace,kkt_buffer, settings) end function solve(ws::AVIWorkspace) - exitflag = -4 + exitflag = 0 λstar, ASstar = zeros(0), Int[] n = length(ws.f) ϵp, ϵd = ws.settings.primal_tol, ws.settings.dual_tol α,β = ws.settings.alpha, ws.settings.beta - tot_iter,outer_iter = 0,0 + tot_iter,outer_iter,nkkt = 0,0,0 @inbounds for k in 1:ws.settings.iter_limit ws.H2xpf .= ws.f mul!(ws.H2xpf, ws.H2mI, ws.x, 1.0, 1.0) @@ -75,10 +88,9 @@ function solve(ws::AVIWorkspace) # Same AS -> Check if KKT conditions are satisfied if info.iterations == 1 + nkkt += 1 ASu,ASl = _get_AS(info.λ,ϵd) - z,AS = _solve_kkt(ws.H,ws.f,ws.A,ws.bu,ws.bl,ASu,ASl,ws.kkt_buffer) - xt = @view z[1:n] - λ = @view z[n+1:end] + xt,λ,AS = _solve_kkt_schur(ws,ASu,ASl) nu = length(ASu) if all(λ[i] > -ϵd for i in 1:nu) && all(λ[i] < ϵd for i in nu+1:length(λ)) Ax = ws.A*xt @@ -95,8 +107,9 @@ function solve(ws::AVIWorkspace) mul!(ws.y, ws.H2, ws.x, 1.0, 1.0) ldiv!(ws.x, ws.H2pI, ws.y) end + k == ws.settings.iter_limit && (exitflag = -4) end - info = (status=flag2status[exitflag], AS=ASstar, outer_iterations=outer_iter, iterations=tot_iter) + info = (status=flag2status[exitflag], AS=ASstar, outer_iterations=outer_iter, iterations=tot_iter, nkkt=nkkt) return ws.x,λstar,exitflag,info end @@ -112,17 +125,40 @@ function _get_AS(λ,dual_tol) return ASu,ASl end -function _solve_kkt(H,f,A,bu,bl,ASu,ASl,kkt_buffer) +function _solve_kkt_schur(ws::AVIWorkspace,ASu,ASl) + # Prepare some helpers AS = [ASu;ASl] - n = length(f) - nkkt = n+length(AS) - resize!(kkt_buffer,nkkt^2) - K = reshape(view(kkt_buffer,1:nkkt^2),(nkkt,nkkt)) - @views K[1:n,1:n] = H - @views K[n+1:end,1:n] = A[AS,:] - @views K[1:n,n+1:end] = A[AS,:]' - K[n+1:end,n+1:end] .= 0 - return K\[-f;bu[ASu];bl[ASl]], AS + n,nAS = length(ws.f),length(AS) + + # Juggle some buffers + resize!(ws.kkt_buffer,nAS*(nAS+n+1)) # Is this necessary? + invH_At = reshape(view(ws.kkt_buffer,1:nAS*n),(n,nAS)) + offset = nAS*n + S = reshape(view(ws.kkt_buffer,offset+1:offset+nAS^2),(nAS,nAS)) + offset +=nAS^2 + λ = view(ws.kkt_buffer,offset+1:offset+nAS) + + AAS = ws.A[AS,:] + + # Schur complement + @inbounds for (i,id) in enumerate(AS) + @views invH_At[:,i] .= ws.A[id,:] + end + #transpose!(invH_At, AAS) + ldiv!(ws.H_factor,invH_At) + mul!(S, AAS, invH_At) + S_fact = lu!(S) + + # Solve for λ + @views λ .= AAS * ws.x_unc + @views λ[1:length(ASu)] .-= ws.bu[ASu] + @views λ[length(ASu)+1:end] .-= ws.bl[ASl] + ldiv!(S_fact,λ) + + # Solve for x + x = ws.x_unc - invH_At * λ + + return x,λ,AS end """ @@ -149,7 +185,7 @@ x ∈ {x: blower ≤ Ax ≤ bupper} such that: * `info` - tuple containing extra information from the solver such as status, active set and number of iterations. """ -function solve_avi(H::AbstractMatrix,f::AbstractVector,A::AbstractMatrix,bupper::AbstractVector,blower::AbstractVector=Float64[]; x0 = zeros(0), settings=AVISettings()) - exitflag, ws = setup_avi(H,f,A,bupper,blower;x0=zeros(0),settings) +function solve_avi(H::AbstractMatrix,f::AbstractVector,A::AbstractMatrix,bupper::AbstractVector,blower::AbstractVector=Float64[], sense=Cint[]; x0 = zeros(0), settings=AVISettings()) + exitflag, ws = setup_avi(H,f,A,bupper,blower,sense;x0,settings) return solve(ws) end From 2901d73e4260b5c70a29156c23934e46b8daf700 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Thu, 15 Jan 2026 21:20:13 +0100 Subject: [PATCH 8/9] Use simple kkt solve for now --- interfaces/daqp-julia/src/avi.jl | 95 +++++++++++++++++++------------- 1 file changed, 57 insertions(+), 38 deletions(-) diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index 37656a8..2a19ca9 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -45,10 +45,10 @@ function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, sett H2 = H1+(H-H') / 2 H1pI,H2mI,H2pI = H1+I, H2-I, lu!(H2+I) - H_factor = lu(H) - x_unc = -(H_factor\f) + #H_factor = lu(H) + H_factor = lu(zeros(0,0)) - x = isempty(x0) ? copy(x_unc) : copy(x0) + x = isempty(x0) ? fill(NaN,n) : copy(x0) y,H2xpf = zeros(n),zeros(n) kkt_buffer = Float64[] @@ -65,7 +65,7 @@ function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, sett end return exitflag, AVIWorkspace(H,f,A,bu,bl, - H_factor,x_unc, + H_factor,zeros(n), H1pI,H2mI,H2pI,H2,H2xpf, x,y,daqp_workspace,kkt_buffer, settings) @@ -78,6 +78,7 @@ function solve(ws::AVIWorkspace) ϵp, ϵd = ws.settings.primal_tol, ws.settings.dual_tol α,β = ws.settings.alpha, ws.settings.beta tot_iter,outer_iter,nkkt = 0,0,0 + isnan(ws.x[1]) && (ws.x .= -ws.H\ws.f) # Use unconstrained optimum as x0 @inbounds for k in 1:ws.settings.iter_limit ws.H2xpf .= ws.f mul!(ws.H2xpf, ws.H2mI, ws.x, 1.0, 1.0) @@ -90,7 +91,8 @@ function solve(ws::AVIWorkspace) if info.iterations == 1 nkkt += 1 ASu,ASl = _get_AS(info.λ,ϵd) - xt,λ,AS = _solve_kkt_schur(ws,ASu,ASl) + #xt,λ,AS = _solve_kkt_schur(ws,ASu,ASl) + xt,λ,AS = _solve_kkt(ws,ASu,ASl) nu = length(ASu) if all(λ[i] > -ϵd for i in 1:nu) && all(λ[i] < ϵd for i in nu+1:length(λ)) Ax = ws.A*xt @@ -125,40 +127,57 @@ function _get_AS(λ,dual_tol) return ASu,ASl end -function _solve_kkt_schur(ws::AVIWorkspace,ASu,ASl) - # Prepare some helpers +# TODO kkt_schur is faster in isolation, but leads to slower total execution +#function _solve_kkt_schur(ws::AVIWorkspace,ASu,ASl) +# # Prepare some helpers +# AS = [ASu;ASl] +# n,nAS = length(ws.f),length(AS) +# +# # Juggle some buffers +# resize!(ws.kkt_buffer,nAS*(nAS+n+1)+n) # Is this necessary? +# invH_At = reshape(view(ws.kkt_buffer,1:nAS*n),(n,nAS)) +# offset = nAS*n +# S = reshape(view(ws.kkt_buffer,offset+1:offset+nAS^2),(nAS,nAS)) +# offset +=nAS^2 +# λ = view(ws.kkt_buffer,offset+1:offset+nAS) +# offset +=nAS +# x = view(ws.kkt_buffer,offset+1:offset+n) +# +# AAS = ws.A[AS,:] +# +# # Schur complement +# @inbounds for (i,id) in enumerate(AS) +# @views invH_At[:,i] .= ws.A[id,:] +# end +# ldiv!(ws.H_factor,invH_At) +# mul!(S, AAS, invH_At) +# S_fact = lu!(S) +# +# # Solve for λ +# @views λ .= AAS * ws.x_unc +# @views λ[1:length(ASu)] .-= ws.bu[ASu] +# @views λ[length(ASu)+1:end] .-= ws.bl[ASl] +# ldiv!(S_fact,λ) +# +# # Solve for x +# x .= ws.x_unc .- invH_At * λ +# +# return x,λ,AS +#end +function _solve_kkt(ws,ASu,ASl) AS = [ASu;ASl] - n,nAS = length(ws.f),length(AS) - - # Juggle some buffers - resize!(ws.kkt_buffer,nAS*(nAS+n+1)) # Is this necessary? - invH_At = reshape(view(ws.kkt_buffer,1:nAS*n),(n,nAS)) - offset = nAS*n - S = reshape(view(ws.kkt_buffer,offset+1:offset+nAS^2),(nAS,nAS)) - offset +=nAS^2 - λ = view(ws.kkt_buffer,offset+1:offset+nAS) - - AAS = ws.A[AS,:] - - # Schur complement - @inbounds for (i,id) in enumerate(AS) - @views invH_At[:,i] .= ws.A[id,:] - end - #transpose!(invH_At, AAS) - ldiv!(ws.H_factor,invH_At) - mul!(S, AAS, invH_At) - S_fact = lu!(S) - - # Solve for λ - @views λ .= AAS * ws.x_unc - @views λ[1:length(ASu)] .-= ws.bu[ASu] - @views λ[length(ASu)+1:end] .-= ws.bl[ASl] - ldiv!(S_fact,λ) - - # Solve for x - x = ws.x_unc - invH_At * λ - - return x,λ,AS + n = length(ws.f) + nkkt = n+length(AS) + resize!(ws.kkt_buffer,nkkt^2) + K = reshape(view(ws.kkt_buffer,1:nkkt^2),(nkkt,nkkt)) + @views K[1:n,1:n] = ws.H + @views K[n+1:end,1:n] = ws.A[AS,:] + @views K[1:n,n+1:end] = ws.A[AS,:]' + K[n+1:end,n+1:end] .= 0 + z = K\[-ws.f;ws.bu[ASu];ws.bl[ASl]] + x = @view z[1:n] + λ = @view z[n+1:end] + return x,λ,AS end """ From 5e07f08d1fe6bb79ee7e438d1429126cd1f67c9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Thu, 15 Jan 2026 21:21:01 +0100 Subject: [PATCH 9/9] Bunmp version --- interfaces/daqp-julia/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/daqp-julia/Project.toml b/interfaces/daqp-julia/Project.toml index fc79260..80a1148 100644 --- a/interfaces/daqp-julia/Project.toml +++ b/interfaces/daqp-julia/Project.toml @@ -1,7 +1,7 @@ name = "DAQPBase" uuid = "833476c3-a8c0-4073-9b64-6473509843fe" authors = ["Daniel Arnström "] -version = "0.3.4" +version = "0.3.5" [deps] DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d"