From 779b40965abb9b9c24dea316794f70f3cf57ee8b Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Fri, 28 Aug 2026 14:02:44 -0400 Subject: [PATCH 1/7] use LAPACK in compactBFGS system solve --- ext/HSL/ma57_workspace.jl | 18 +++++++++++------- ext/LDLFactorizations/ldlt.jl | 22 +++++++++++++--------- src/linear_algebra/K2.jl | 2 -- src/linear_algebra/mumps.jl | 18 ++++++++++-------- test/test-cutest.jl | 4 +++- 5 files changed, 37 insertions(+), 27 deletions(-) diff --git a/ext/HSL/ma57_workspace.jl b/ext/HSL/ma57_workspace.jl index b4421df9..400a47a7 100644 --- a/ext/HSL/ma57_workspace.jl +++ b/ext/HSL/ma57_workspace.jl @@ -2,6 +2,7 @@ mutable struct PenaltyMA57Workspace{ WP<:Ma57, K2<:AbstractMatrix, V<:AbstractVector, + VI<:Union{Nothing,AbstractVector}, T<:Real, } <: AbstractHSLWorkspace M::WP @@ -10,6 +11,7 @@ mutable struct PenaltyMA57Workspace{ work::V _qn_work::V dx::V + _ipiv::VI # For CompactBFGS LU factorization σ::T n::Int m::Int @@ -42,6 +44,7 @@ function construct_ma57_workspace( similar(u1, 4*length(u1)), similar(u1, 0), similar(u1), + nothing zero(T), n, m, @@ -72,6 +75,7 @@ function construct_ma57_workspace( similar(u1, 4*length(u1)), similar(u1, 2*length(u1)*H.B._mem), similar(u1), + Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), zero(T), n, m, @@ -214,12 +218,12 @@ function solve_system!( B = workspace.H.B n, m = workspace.n, workspace.m p = min(B._insert - 1, B._mem) - x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2 + x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1 Z1, Z2 = H.Z1, H.Z2 Uk = @view B.Uk[:, 1:p] Vk = @view B.Vk[:, 1:p] - # Step 0: Write (#TODO: we can use easily use QRMumps instead of LDLFactorization here...) + # Step 0: Write # [B Aᵀ] = [σI+ξI Aᵀ] + [-U V]([U V])ᵀ # [A -αI] = [A -αI] + [ 0 0]([0 0]) # Hence, @@ -316,9 +320,9 @@ function solve_system!( # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] # ( [A -αI] ) [y₁] # using Julia LinearALgebra's lu! - F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ? - @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)]) - if any(isnan, y2) + @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) + @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) + if any(isnan, @view y1[1:(2*p)]) workspace.status = :failed return end @@ -326,8 +330,8 @@ function solve_system!( # Step 6: Compute # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2))) + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) # Step 7: Solve # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] diff --git a/ext/LDLFactorizations/ldlt.jl b/ext/LDLFactorizations/ldlt.jl index 3cd46e1d..ad4495dd 100644 --- a/ext/LDLFactorizations/ldlt.jl +++ b/ext/LDLFactorizations/ldlt.jl @@ -2,6 +2,7 @@ mutable struct PenaltyLDLTWorkspace{ WP<:LDLFactorization, K2<:AbstractMatrix, V<:AbstractVector, + VI<:Union{Nothing,AbstractVector}, T<:Real, } <: AbstractLDLTWorkspace M::WP @@ -9,6 +10,7 @@ mutable struct PenaltyLDLTWorkspace{ x::V dx::V r::V + _ipiv::VI # For CompactBFGS LU factorization σ::T n::Int m::Int @@ -39,6 +41,7 @@ function construct_ldlt_workspace( similar(u1), similar(u1), similar(u1), + nothing, zero(T), n, m, @@ -60,6 +63,7 @@ function construct_ldlt_workspace( similar(u1), similar(u1), similar(u1), + Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), zero(T), n, m, @@ -268,12 +272,12 @@ function solve_system!( B = workspace.H.B n, m = workspace.n, workspace.m p = min(B._insert - 1, B._mem) - x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2 + x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1 Z1, Z2 = H.Z1, H.Z2 Uk = @view B.Uk[:, 1:p] Vk = @view B.Vk[:, 1:p] - # Step 0: Write (#TODO: we can use easily use QRMumps instead of LDLFactorization here...) + # Step 0: Write # [B Aᵀ] = [σI+ξI Aᵀ] + [-U V]([U V])ᵀ # [A -αI] = [A -αI] + [ 0 0]([0 0]) # Hence, @@ -292,7 +296,7 @@ function solve_system!( return end - # Step 2: Compute # TODO: allow for iterative refinement + # Step 2: Compute # [x₁] = [σI+ξI Aᵀ]⁻¹[u] # [x₁] = [A -αI] [u] ldiv!(x1, workspace.M, u) @@ -339,10 +343,10 @@ function solve_system!( # Step 5: Solve # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] # ( [A -αI] ) [y₁] - # using Julia LinearALgebra's lu! - F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ? - @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)]) - if any(isnan, y2) + # using LAPACK + @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) + @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) + if any(isnan, @view y1[1:(2*p)]) workspace.status = :failed return end @@ -350,8 +354,8 @@ function solve_system!( # Step 6: Compute # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2))) + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) # Step 7: Solve # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] diff --git a/src/linear_algebra/K2.jl b/src/linear_algebra/K2.jl index c7e84255..2381d285 100644 --- a/src/linear_algebra/K2.jl +++ b/src/linear_algebra/K2.jl @@ -15,7 +15,6 @@ mutable struct CompactBFGSK2{ x2::V x3::V y1::V - y2::V end function K2( @@ -128,7 +127,6 @@ function K2( zeros(T, n+m), zeros(T, n+m), zeros(T, 2*B._mem), - zeros(T, 2*B._mem), ) end diff --git a/src/linear_algebra/mumps.jl b/src/linear_algebra/mumps.jl index c98e6f4a..ca539738 100644 --- a/src/linear_algebra/mumps.jl +++ b/src/linear_algebra/mumps.jl @@ -2,11 +2,13 @@ mutable struct PenaltyMUMPSWorkspace{ WP<:Mumps, K2<:AbstractMatrix, V<:AbstractVector, + VI<:Union{Nothing,AbstractVector}, T<:Real, } <: AbstractMUMPSWorkspace M::WP H::K2 x::V + _ipiv::VI # For CompactBFGS LU factorization σ::T n::Int m::Int @@ -84,7 +86,7 @@ function construct_mumps_workspace( S.rhs = pointer(x) S._y_gc_haven = x - return PenaltyMUMPSWorkspace(S, H, x, zero(T), n, m, :uninitialized, false, 0) + return PenaltyMUMPSWorkspace(S, H, x, nothing, zero(T), n, m, :uninitialized, false, 0) end function construct_mumps_workspace( @@ -145,7 +147,7 @@ function construct_mumps_workspace( S.rhs = pointer(x) S._y_gc_haven = x - return PenaltyMUMPSWorkspace(S, H, x, zero(T), n, m, :uninitialized, false, 0) + return PenaltyMUMPSWorkspace(S, H, x, Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), zero(T), n, m, :uninitialized, false, 0) end function update_workspace!( @@ -280,7 +282,7 @@ function solve_system!( mumps = workspace.M n, m = workspace.n, workspace.m p = min(B._insert - 1, B._mem) - x1, x2, x3, y1, y2 = H.x1, H.x2, H.x3, H.y1, H.y2 + x1, x2, x3, y1 = H.x1, H.x2, H.x3, H.y1 Z1, Z2 = H.Z1, H.Z2 Uk = @view B.Uk[:, 1:p] @@ -388,9 +390,9 @@ function solve_system!( # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] # ( [A -αI] ) [y₁] # using Julia LinearALgebra's lu! - F = lu!(Z2[1:(2*p), 1:(2*p)], check = false) # FIXME ? - @views ldiv!(y2[1:(2*p)], F, y1[1:(2*p)]) - if any(isnan, y2) + @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) + @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) + if any(isnan, @view y1[1:(2*p)]) workspace.status = :failed return end @@ -398,8 +400,8 @@ function solve_system!( # Step 6: Compute # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y2[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y2[1:p], -one(eltype(y2)), one(eltype(y2))) + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) # Step 7: Solve # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] diff --git a/test/test-cutest.jl b/test/test-cutest.jl index d4bfcf22..05c732d7 100644 --- a/test/test-cutest.jl +++ b/test/test-cutest.jl @@ -102,7 +102,9 @@ function test_problem( solver = L2PenaltySolver(LBFGS_model, linear_solver = linear_solver) stats_optimized = PeneloptExecutionStats(LBFGS_model) - solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0) + @test @wrappedallocs( + solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0) + ) == 0 stats_optimized.solution = recover_full_solution(LBFGS_model, stats_optimized.solution) From 142a530cbd2aae3278f2d45ed4611886ecc139fe Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Fri, 28 Aug 2026 14:22:15 -0400 Subject: [PATCH 2/7] Use unsafe associate rhs --- src/linear_algebra/mumps.jl | 6 +++--- test/test-cutest.jl | 8 ++++++-- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/linear_algebra/mumps.jl b/src/linear_algebra/mumps.jl index ca539738..f8a65205 100644 --- a/src/linear_algebra/mumps.jl +++ b/src/linear_algebra/mumps.jl @@ -334,7 +334,7 @@ function solve_system!( # [x₁] = [A -αI] [u] x1 .= u - MUMPS.associate_rhs!(mumps, x1) + MUMPS.associate_rhs!(mumps, x1; unsafe = true) MUMPS.mumps_solve!(x1, mumps; rhs_changed = true) # MUMPS infog(1): a negative value is an error in the factorization. @@ -362,7 +362,7 @@ function solve_system!( @views Z1[1:n, 1:p] .= Uk .* (-1) @views Z1[1:n, (p+1):(2*p)] .= Vk - MUMPS.associate_rhs!(mumps, Z1) + MUMPS.associate_rhs!(mumps, Z1; unsafe = true) MUMPS.mumps_solve!(Z1, mumps; rhs_changed = true) # MUMPS infog(1): a negative value is an error in the factorization. @@ -407,7 +407,7 @@ function solve_system!( # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] # [x₃] = [A -αI] [x₂] x3 .= x2 - MUMPS.associate_rhs!(mumps, x3) + MUMPS.associate_rhs!(mumps, x3; unsafe = true) MUMPS.mumps_solve!(x3, mumps; rhs_changed = true) # MUMPS infog(1): a negative value is an error in the factorization. diff --git a/test/test-cutest.jl b/test/test-cutest.jl index 05c732d7..63ba7cc1 100644 --- a/test/test-cutest.jl +++ b/test/test-cutest.jl @@ -102,9 +102,13 @@ function test_problem( solver = L2PenaltySolver(LBFGS_model, linear_solver = linear_solver) stats_optimized = PeneloptExecutionStats(LBFGS_model) - @test @wrappedallocs( + if linear_solver != "mumps" + @test @wrappedallocs( + solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0) + ) == 0 + else solve!(solver, LBFGS_model, stats_optimized, atol = 1e-3, rtol = 1e-3, τ0 = 1.0) - ) == 0 + end stats_optimized.solution = recover_full_solution(LBFGS_model, stats_optimized.solution) From 8d1fc8950c8e62fbadc5571f7e38fce46f9721d0 Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Fri, 28 Aug 2026 14:39:39 -0400 Subject: [PATCH 3/7] Prevent storing Jacobian twice --- src/types/PenalizedProblem.jl | 2 -- src/types/ShiftedPenalizedProblem.jl | 21 +++++++++------------ src/types/norm/CompositeNormL2.jl | 7 +------ src/types/norm/ShiftedCompositeNormL2.jl | 14 ++------------ 4 files changed, 12 insertions(+), 32 deletions(-) diff --git a/src/types/PenalizedProblem.jl b/src/types/PenalizedProblem.jl index bbcef241..da75f830 100644 --- a/src/types/PenalizedProblem.jl +++ b/src/types/PenalizedProblem.jl @@ -37,14 +37,12 @@ function L2PenalizedProblem(nlp::AbstractNLPModel{T,S}) where {T,S} A = SparseMatrixCOO(nlp.meta.ncon, nlp.meta.nvar, rows, cols, vals) b = similar(x0, eltype(x0), nlp.meta.ncon) - store_previous_jacobian = isa(nlp, QuasiNewtonModel) ? true : false penalty = CompositeNormL2( one(T), (c, x) -> cons!(nlp, x, c), (j, x) -> jac_coord!(nlp, x, j.vals), A, b, - store_previous_jacobian = store_previous_jacobian, ) return L2PenalizedProblem(nlp, penalty, nlp.meta) end diff --git a/src/types/ShiftedPenalizedProblem.jl b/src/types/ShiftedPenalizedProblem.jl index d3dc0779..acbe71ca 100644 --- a/src/types/ShiftedPenalizedProblem.jl +++ b/src/types/ShiftedPenalizedProblem.jl @@ -31,7 +31,6 @@ mutable struct ShiftedL2PenalizedProblem{ h::H parent::P meta::meta - _qn_∇f_prev::SN _qn_y::SN _qn_x_prev::SN _is_first_shift::Bool @@ -56,7 +55,6 @@ function ShiftedL2PenalizedProblem( ψ, penalty_nlp, penalty_nlp.meta, - zero(∇f), similar(∇f), zero(∇f), true, @@ -163,23 +161,28 @@ function shift!( } nlp, h = shifted_penalty_nlp.parent.model, shifted_penalty_nlp.parent.h φ, ψ = shifted_penalty_nlp.model, shifted_penalty_nlp.h - qn_y, qn_g_prev, qn_x_prev = shifted_penalty_nlp._qn_y, - shifted_penalty_nlp._qn_∇f_prev, + qn_y, qn_x_prev = shifted_penalty_nlp._qn_y, shifted_penalty_nlp._qn_x_prev is_first_shift = shifted_penalty_nlp._is_first_shift qn_s = qn_x_prev g, B = φ.data.c, φ.data.H + if !is_first_shift + qn_y .= g + if !isnothing(y) + mul!(qn_y, ψ.A', y, -one(T), -one(T)) # y = - g_prev - J(x)_prev^T λ + end + end + isnothing(∇f) ? grad!(nlp, x, g) : (g .= ∇f) shift!(ψ, x, J = J, c = c) # Update the approximation. if !is_first_shift - @. qn_y = g - qn_g_prev + @. qn_y .+= g if !isnothing(y) mul!(qn_y, ψ.A', y, one(T), one(T)) # y = y + J(x)^T λ - mul!(qn_y, ψ.A_prev', y, -one(T), one(T)) # y = y - J(x)_prev^T λ end qn_s .= x .- qn_x_prev @@ -189,11 +192,7 @@ function shift!( shifted_penalty_nlp._is_first_shift = false end - # Copy the gradient and Jacobian. - qn_g_prev .= g - ψ.A_prev.vals .= ψ.A.vals qn_x_prev .= x - end function shift!( @@ -254,9 +253,7 @@ function reset!( nlp, h = shifted_penalty_nlp.parent.model, shifted_penalty_nlp.parent.h φ, ψ = shifted_penalty_nlp.model, shifted_penalty_nlp.h x_prev = shifted_penalty_nlp._qn_x_prev .= 0 - g_prev = shifted_penalty_nlp._qn_∇f_prev .= 0 - ψ.A_prev.vals .= ψ.A.vals LinearOperators.reset!(φ.data.H) shifted_penalty_nlp._is_first_shift = true end diff --git a/src/types/norm/CompositeNormL2.jl b/src/types/norm/CompositeNormL2.jl index 12a45949..734bf99c 100644 --- a/src/types/norm/CompositeNormL2.jl +++ b/src/types/norm/CompositeNormL2.jl @@ -1,7 +1,7 @@ abstract type AbstractCompositeNorm end @doc raw""" - CompositeNormL2(λ, c!, J!, A, b; store_previous_jacobian::Bool = false) + CompositeNormL2(λ, c!, J!, A, b) Returns function `c` composed with the `ℓ₂` norm: ```math @@ -21,8 +21,6 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`. c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real}) J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real}) ``` -Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`. -This is particularly useful for quasi-Newton updates in the context of constrained optimization. """ mutable struct CompositeNormL2{ T<:Real, @@ -36,7 +34,6 @@ mutable struct CompositeNormL2{ J!::F1 A::M b::V - store_previous_jacobian::Bool function CompositeNormL2( λ::T, @@ -44,7 +41,6 @@ mutable struct CompositeNormL2{ J!::Function, A::AbstractMatrix{T}, b::AbstractVector{T}; - store_previous_jacobian::Bool = false, ) where {T<:Real} λ > 0 || error("CompositeNormL2: λ should be positive") length(b) == size(A, 1) || error( @@ -56,7 +52,6 @@ mutable struct CompositeNormL2{ J!, A, b, - store_previous_jacobian, ) end end diff --git a/src/types/norm/ShiftedCompositeNormL2.jl b/src/types/norm/ShiftedCompositeNormL2.jl index 4b3126ad..d4bb13d8 100644 --- a/src/types/norm/ShiftedCompositeNormL2.jl +++ b/src/types/norm/ShiftedCompositeNormL2.jl @@ -1,7 +1,7 @@ abstract type AbstractShiftedCompositeNorm end @doc raw""" - ShiftedCompositeNormL2(h, c!, J!, A, b; store_previous_jacobian::Bool = false) + ShiftedCompositeNormL2(h, c!, J!, A, b) Returns the shift of a function `c` composed with the `ℓ₂` norm (see CompositeNormL2.jl). Here, `c` is linearized i.e, `c(x + s) ≈ c(x) + J(x)s`. @@ -22,23 +22,18 @@ such that `J` is the Jacobian of `c`. It is expected that `m ≤ n`. c!(b <: AbstractVector{Real}, xk <: AbstractVector{Real}) J!(A <: AbstractSparseMatrixCOO{Real, Integer}, xk <: AbstractVector{Real}) ``` -Moreover, if you want shifted instances of the operator to store the previous Jacobian on each shift, you can specify `store_previous_jacobian = true`. -In this case, each time a shift is performed, the previous Jacobian is stored in the `A_prev` field. -This is particularly useful for quasi-Newton updates in the context of constrained optimization. """ mutable struct ShiftedCompositeNormL2{ T<:Real, F0<:Function, F1<:Function, M<:AbstractMatrix{T}, - N<:Union{Nothing,M}, V<:AbstractVector{T}, } <: AbstractShiftedCompositeNorm h::NormL2{T} c!::F0 J!::F1 A::M - A_prev::N # (Optional) can be used to store the previous Jacobian, useful for quasi-Newton approximations b::V g::V function ShiftedCompositeNormL2( @@ -47,7 +42,6 @@ mutable struct ShiftedCompositeNormL2{ J!::Function, A::AbstractMatrix{T}, b::AbstractVector{T}; - store_previous_jacobian::Bool = false, ) where {T<:Real} if length(b) != size(A, 1) error( @@ -55,15 +49,13 @@ mutable struct ShiftedCompositeNormL2{ ) end - A_prev = store_previous_jacobian ? copy(A) : nothing g = similar(b) - new{T,typeof(c!),typeof(J!),typeof(A),typeof(A_prev),typeof(b)}( + new{T,typeof(c!),typeof(J!),typeof(A),typeof(b)}( NormL2(λ), c!, J!, A, - A_prev, b, g, ) @@ -83,7 +75,6 @@ shifted( ψ.J!, A, b, - store_previous_jacobian = ψ.store_previous_jacobian, ) end @@ -103,7 +94,6 @@ function shift!( J = nothing, c = nothing, ) where {R<:Real} - !isnothing(ψ.A_prev) && (ψ.A_prev.vals .= ψ.A.vals) # Update previous Jacobian if necessary isnothing(c) ? ψ.c!(ψ.b, shift) : (ψ.b .= c) isnothing(J) ? ψ.J!(ψ.A, shift) : (ψ.A.vals .= J.vals) return ψ From 55c8b414ac5b6788450a2732c81a8e65ac42d70a Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Fri, 28 Aug 2026 14:52:36 -0400 Subject: [PATCH 4/7] Fix constructors --- src/types/ShiftedPenalizedProblem.jl | 1 - test/instances/instance-generator.jl | 1 - test/instances/instance-reader.jl | 1 - 3 files changed, 3 deletions(-) diff --git a/src/types/ShiftedPenalizedProblem.jl b/src/types/ShiftedPenalizedProblem.jl index acbe71ca..83842e21 100644 --- a/src/types/ShiftedPenalizedProblem.jl +++ b/src/types/ShiftedPenalizedProblem.jl @@ -89,7 +89,6 @@ function ShiftedL2PenalizedProblem( penalty_nlp.meta, nothing, nothing, - nothing, true, ) end diff --git a/test/instances/instance-generator.jl b/test/instances/instance-generator.jl index 3c1d48a0..8620d640 100644 --- a/test/instances/instance-generator.jl +++ b/test/instances/instance-generator.jl @@ -84,7 +84,6 @@ function generate_instance( model.meta, nothing, nothing, - nothing, true, ), Dict(:u => u, :y => y, :tau => tau) diff --git a/test/instances/instance-reader.jl b/test/instances/instance-reader.jl index b52c5d83..41ca3e61 100644 --- a/test/instances/instance-reader.jl +++ b/test/instances/instance-reader.jl @@ -52,7 +52,6 @@ function read_instance(file::String; type = Float64, Hessian_modifier = H -> H) model.meta, nothing, nothing, - nothing, true, ) end From 5e03008922e92b60efa89d4b5c92963bca52ef92 Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Sun, 30 Aug 2026 14:25:38 -0400 Subject: [PATCH 5/7] Start Lapack --- src/Penelopt.jl | 1 + src/linear_algebra/lapack.jl | 47 ++++++++++++++++++++++++++++++++++++ 2 files changed, 48 insertions(+) create mode 100644 src/linear_algebra/lapack.jl diff --git a/src/Penelopt.jl b/src/Penelopt.jl index c5aac526..917efef6 100644 --- a/src/Penelopt.jl +++ b/src/Penelopt.jl @@ -47,6 +47,7 @@ include("types/pre-processing/FixedVariable.jl") include("linear_algebra/K2.jl") include("linear_algebra/construct_workspace.jl") include("linear_algebra/mumps.jl") +include("linear_algebra/lapack.jl") include("types/PenalizedProblem.jl") include("types/ShiftedPenalizedProblem.jl") diff --git a/src/linear_algebra/lapack.jl b/src/linear_algebra/lapack.jl new file mode 100644 index 00000000..22ef09df --- /dev/null +++ b/src/linear_algebra/lapack.jl @@ -0,0 +1,47 @@ +mutable struct PenaltyLacpackWorkspace{T <: Real} <: PenaltyDirectWorkspace + m::BlasInt # rows of A + n::BlasInt # cols of A (= rows/cols since we assume square) + a::Matrix{T} # the factorized matrix, overwritten in-place by getrf! + lda::BlasInt # leading dimension of a + ipiv::Vector{BlasInt} # pivot indices, length min(m,n) + info::Base.RefValue{BlasInt} # LAPACK status output + nrhs::BlasInt # number of right-hand sides for getrs! + ldb::BlasInt # leading dimension of b +end + +function PenaltyLacpackWorkspace(A::Matrix{T}; nrhs::Integer = 1) where {T <: Real} + m, n = size(A) + m == n || throw(ArgumentError("getrs! requires a square matrix, got $m x $n")) + lda = max(1, m) + ipiv = Vector{BlasInt}(undef, m) + info = Ref{BlasInt}(0) + ldb = max(1, n) + return PenaltyLacpackWorkspace{T}(m, n, A, lda, ipiv, info, BlasInt(nrhs), BlasInt(ldb)) +end + +for (getrf, getrs, T) in + ((:sgetrf_, :sgetrs_, :Float32), (:dgetrf_, :dgetrs_, :Float64)) + @eval begin + + # getrf + function getrf!( + workspace::PenaltyLacpackWorkspace{ST}, + ) where{ST <: $T} + return ccall((@blasfunc($getrf), libblastrampoline), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}), + workspace.m, workspace.n, workspace.a, workspace.lda, workspace.ipiv, workspace.info) + end + + # getrs + function getrs!( + trans::Char, + workspace::PenaltyLacpackWorkspace{ST}, + b::AbstractVector{ST}, + ) where{ST <: $T} + return ccall((@blasfunc($getrs), libblastrampoline), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, + Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong), + trans, workspace.n, nrhs, workspace.a, workspace.lda, workspace.ipiv, b, workspace.ldb, workspace.info, 1) + end + end +end \ No newline at end of file From 36bfae1d6cb8031747090cfba9828b619993b2f1 Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Sun, 30 Aug 2026 15:23:46 -0400 Subject: [PATCH 6/7] Use homemade LAPACK functions --- ext/HSL/ma57_workspace.jl | 150 +++++++++++++++++++--------------- ext/LDLFactorizations/ldlt.jl | 134 +++++++++++++++++------------- ext/PeneloptHSLExt.jl | 5 ++ src/Penelopt.jl | 4 + src/linear_algebra/lapack.jl | 60 ++++++-------- src/linear_algebra/mumps.jl | 144 +++++++++++++++++--------------- 6 files changed, 274 insertions(+), 223 deletions(-) diff --git a/ext/HSL/ma57_workspace.jl b/ext/HSL/ma57_workspace.jl index 400a47a7..d4ae1dfe 100644 --- a/ext/HSL/ma57_workspace.jl +++ b/ext/HSL/ma57_workspace.jl @@ -10,6 +10,7 @@ mutable struct PenaltyMA57Workspace{ x::V work::V _qn_work::V + _info::Base.RefValue{BlasInt} # For CompactBFGS LU factorization dx::V _ipiv::VI # For CompactBFGS LU factorization σ::T @@ -44,7 +45,8 @@ function construct_ma57_workspace( similar(u1, 4*length(u1)), similar(u1, 0), similar(u1), - nothing + nothing, + Ref{BlasInt}(0), zero(T), n, m, @@ -76,6 +78,7 @@ function construct_ma57_workspace( similar(u1, 2*length(u1)*H.B._mem), similar(u1), Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), + Ref{BlasInt}(0), zero(T), n, m, @@ -277,79 +280,94 @@ function solve_system!( return end - # Step 3: Compute - # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] - # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] - @views mul!(y1[1:p], Uk', x1[1:n]) - @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) + if p > 0 + # Step 3: Compute + # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] + # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] + @views mul!(y1[1:p], Uk', x1[1:n]) + @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) - # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) - # ( [A -αI] ) - # Step 4.1: Compute - # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] - # Z₁ = [A -αI] E = [A -αI] [ 0 0] - Z1 .= 0 - @views Z1[1:n, 1:p] .= Uk .* (-1) - @views Z1[1:n, (p+1):(2*p)] .= Vk - try - ma57_solve!(workspace.M, Z1, workspace._qn_work) - catch e - !(e isa HSL.Ma57Exception) && rethrow(e) - end - if any(isnan, Z1) || workspace.M.info.info[1] < 0 - workspace.status = :failed - return - end + # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) + # ( [A -αI] ) + # Step 4.1: Compute + # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] + # Z₁ = [A -αI] E = [A -αI] [ 0 0] + Z1 .= 0 - # Step 4.2: Compute - # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] - # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] - Z2 .= 0 - @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) - @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) - - # Step 4.3: Compute - # Z₂ = I + Z₂ - for i = 1:(2*p) - Z2[i, i] += 1 - end + @views Z1[1:n, 1:p] .= Uk .* (-1) + @views Z1[1:n, (p+1):(2*p)] .= Vk + try + ma57_solve!(workspace.M, Z1, workspace._qn_work) + catch e + !(e isa HSL.Ma57Exception) && rethrow(e) + end + if any(isnan, Z1) || workspace.M.info.info[1] < 0 + workspace.status = :failed + return + end - # Step 5: Solve - # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] - # ( [A -αI] ) [y₁] - # using Julia LinearALgebra's lu! - @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) - @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) - if any(isnan, @view y1[1:(2*p)]) - workspace.status = :failed - return - end + # Step 4.2: Compute + # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] + # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] + Z2 .= 0 + @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) + @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) + + # Step 4.3: Compute + # Z₂ = I + Z₂ + for i = 1:(2*p) + Z2[i, i] += 1 + end + + # Step 5: Solve + # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] + # ( [A -αI] ) [y₁] + # using LAPACK + @views info_f = getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info) + if info_f != 0 + workspace.status = :failed + return + end - # Step 6: Compute - # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] - # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) + @views info_s = getrs!('N', BlasInt(2p), BlasInt(1), Z2[1:(2p), 1:(2p)], stride(Z2, 2), + workspace._ipiv, y1[1:(2p)], BlasInt(2p), workspace._info) + if info_s != 0 + workspace.status = :failed + return + end + if any(isnan, @view y1[1:(2*p)]) + workspace.status = :failed + return + end - # Step 7: Solve - # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] - # [x₃] = [A -αI] [x₂] - try - ma57_solve!(workspace.M, x2, x3, workspace.dx, workspace.work, 10) - catch e - !(e isa HSL.Ma57Exception) && rethrow(e) - end - if any(isnan, x3) || workspace.M.info.info[1] < 0 - workspace.status = :failed - return - end + # Step 6: Compute + # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] + # x₂ = E[y₂] = [ 0 0][y₂] = [0] + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) - # Step 8: - # [B Aᵀ]⁻¹[u] = x₁ - x₃ - # [A -αI] [u] = x₁ - x₃ - workspace.x .= x1 .- x3 + # Step 7: Solve + # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] + # [x₃] = [A -αI] [x₂] + try + ma57_solve!(workspace.M, x2, x3, workspace.dx, workspace.work, 10) + catch e + !(e isa HSL.Ma57Exception) && rethrow(e) + end + if any(isnan, x3) || workspace.M.info.info[1] < 0 + workspace.status = :failed + return + end + + # Step 8: + # [B Aᵀ]⁻¹[u] = x₁ - x₃ + # [A -αI] [u] = x₁ - x₃ + workspace.x .= x1 .- x3 + else + workspace.x .= x1 + end end function get_solution!( diff --git a/ext/LDLFactorizations/ldlt.jl b/ext/LDLFactorizations/ldlt.jl index ad4495dd..3119eb3e 100644 --- a/ext/LDLFactorizations/ldlt.jl +++ b/ext/LDLFactorizations/ldlt.jl @@ -11,6 +11,7 @@ mutable struct PenaltyLDLTWorkspace{ dx::V r::V _ipiv::VI # For CompactBFGS LU factorization + _info::Base.RefValue{BlasInt} # For CompactBFGS LU factorization σ::T n::Int m::Int @@ -42,6 +43,7 @@ function construct_ldlt_workspace( similar(u1), similar(u1), nothing, + Ref{BlasInt}(0), zero(T), n, m, @@ -64,6 +66,7 @@ function construct_ldlt_workspace( similar(u1), similar(u1), Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), + Ref{BlasInt}(0), zero(T), n, m, @@ -305,71 +308,86 @@ function solve_system!( return end - # Step 3: Compute - # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] - # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] - @views mul!(y1[1:p], Uk', x1[1:n]) - @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) + if p > 0 + # Step 3: Compute + # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] + # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] + @views mul!(y1[1:p], Uk', x1[1:n]) + @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) - # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) - # ( [A -αI] ) - # Step 4.1: Compute - # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] - # Z₁ = [A -αI] E = [A -αI] [ 0 0] - Z1 .= 0 - @views Z1[1:n, 1:p] .= Uk .* (-1) - @views Z1[1:n, (p+1):(2*p)] .= Vk - ldiv!(workspace.M, Z1) - if any(isnan, Z1) - workspace.status = :failed - return - end + # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) + # ( [A -αI] ) + # Step 4.1: Compute + # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] + # Z₁ = [A -αI] E = [A -αI] [ 0 0] + Z1 .= 0 - # Step 4.2: Compute - # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] - # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] - Z2 .= 0 - @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) - @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) - - # Step 4.3: Compute - # Z₂ = I + Z₂ - for i = 1:(2*p) - Z2[i, i] += 1 - end + @views Z1[1:n, 1:p] .= Uk .* (-1) + @views Z1[1:n, (p+1):(2*p)] .= Vk + ldiv!(workspace.M, Z1) + if any(isnan, Z1) + workspace.status = :failed + return + end - # Step 5: Solve - # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] - # ( [A -αI] ) [y₁] - # using LAPACK - @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) - @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) - if any(isnan, @view y1[1:(2*p)]) - workspace.status = :failed - return - end + # Step 4.2: Compute + # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] + # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] + Z2 .= 0 + @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) + @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) + + # Step 4.3: Compute + # Z₂ = I + Z₂ + for i = 1:(2*p) + Z2[i, i] += 1 + end - # Step 6: Compute - # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] - # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) - - # Step 7: Solve - # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] - # [x₃] = [A -αI] [x₂] - ldiv!(x3, workspace.M, x2) - if any(isnan, x3) - workspace.status = :failed - return - end + # Step 5: Solve + # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] + # ( [A -αI] ) [y₁] + # using LAPACK + @views info_f = getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info) + if info_f != 0 + workspace.status = :failed + return + end - # Step 8: - # [B Aᵀ]⁻¹[u] = x₁ - x₃ - # [A -αI] [u] = x₁ - x₃ - workspace.x .= x1 .- x3 + @views info_s = getrs!('N', BlasInt(2p), BlasInt(1), Z2[1:(2p), 1:(2p)], stride(Z2, 2), + workspace._ipiv, y1[1:(2p)], BlasInt(2p), workspace._info) + if info_s != 0 + workspace.status = :failed + return + end + if any(isnan, @view y1[1:(2*p)]) + workspace.status = :failed + return + end + + # Step 6: Compute + # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] + # x₂ = E[y₂] = [ 0 0][y₂] = [0] + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) + + # Step 7: Solve + # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] + # [x₃] = [A -αI] [x₂] + ldiv!(x3, workspace.M, x2) + if any(isnan, x3) + workspace.status = :failed + return + end + + # Step 8: + # [B Aᵀ]⁻¹[u] = x₁ - x₃ + # [A -αI] [u] = x₁ - x₃ + workspace.x .= x1 .- x3 + else + workspace.x .= x1 + end end function get_solution!(x::V, workspace::PenaltyLDLTWorkspace) where {V<:AbstractVector} diff --git a/ext/PeneloptHSLExt.jl b/ext/PeneloptHSLExt.jl index b6c36327..5c6f5dfb 100644 --- a/ext/PeneloptHSLExt.jl +++ b/ext/PeneloptHSLExt.jl @@ -5,9 +5,14 @@ using Penelopt using LinearAlgebra, SparseMatricesCOO +# Import BLAS functions +import LinearAlgebra.BLAS: @blasfunc +import LinearAlgebra: BlasInt, libblastrampoline + import Penelopt: AbstractHSLWorkspace import Penelopt: construct_ma57_workspace, solve_system!, update_workspace! import Penelopt: get_inertia, get_solution!, get_status +import Penelopt: getrf!, getrs! import Penelopt: set_dual_inertia!, set_primal_inertia! function __init__() diff --git a/src/Penelopt.jl b/src/Penelopt.jl index 917efef6..8c1c1d66 100644 --- a/src/Penelopt.jl +++ b/src/Penelopt.jl @@ -22,6 +22,10 @@ using NLPModels, NLPModelsModifiers using LinearOperators, QuadraticModels, SolverCore, SparseMatricesCOO using MPI, MUMPS +# Import BLAS functions +import LinearAlgebra.BLAS: @blasfunc +import LinearAlgebra: BlasInt, libblastrampoline + import SolverCore: get_status, reset! function __init__() diff --git a/src/linear_algebra/lapack.jl b/src/linear_algebra/lapack.jl index 22ef09df..c1598942 100644 --- a/src/linear_algebra/lapack.jl +++ b/src/linear_algebra/lapack.jl @@ -1,47 +1,37 @@ -mutable struct PenaltyLacpackWorkspace{T <: Real} <: PenaltyDirectWorkspace - m::BlasInt # rows of A - n::BlasInt # cols of A (= rows/cols since we assume square) - a::Matrix{T} # the factorized matrix, overwritten in-place by getrf! - lda::BlasInt # leading dimension of a - ipiv::Vector{BlasInt} # pivot indices, length min(m,n) - info::Base.RefValue{BlasInt} # LAPACK status output - nrhs::BlasInt # number of right-hand sides for getrs! - ldb::BlasInt # leading dimension of b -end - -function PenaltyLacpackWorkspace(A::Matrix{T}; nrhs::Integer = 1) where {T <: Real} - m, n = size(A) - m == n || throw(ArgumentError("getrs! requires a square matrix, got $m x $n")) - lda = max(1, m) - ipiv = Vector{BlasInt}(undef, m) - info = Ref{BlasInt}(0) - ldb = max(1, n) - return PenaltyLacpackWorkspace{T}(m, n, A, lda, ipiv, info, BlasInt(nrhs), BlasInt(ldb)) -end - for (getrf, getrs, T) in ((:sgetrf_, :sgetrs_, :Float32), (:dgetrf_, :dgetrs_, :Float64)) @eval begin - # getrf function getrf!( - workspace::PenaltyLacpackWorkspace{ST}, - ) where{ST <: $T} - return ccall((@blasfunc($getrf), libblastrampoline), Cvoid, - (Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}), - workspace.m, workspace.n, workspace.a, workspace.lda, workspace.ipiv, workspace.info) + m::BlasInt, + n::BlasInt, + a::AbstractMatrix{$T}, + lda::BlasInt, + ipiv::AbstractVector{BlasInt}, + info::Base.RefValue{BlasInt}, + ) + ccall((@blasfunc($getrf), libblastrampoline), Cvoid, + (Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt}), + m, n, a, lda, ipiv, info) + return info[] end - # getrs function getrs!( trans::Char, - workspace::PenaltyLacpackWorkspace{ST}, - b::AbstractVector{ST}, - ) where{ST <: $T} - return ccall((@blasfunc($getrs), libblastrampoline), Cvoid, - (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, - Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong), - trans, workspace.n, nrhs, workspace.a, workspace.lda, workspace.ipiv, b, workspace.ldb, workspace.info, 1) + n::BlasInt, + nrhs::BlasInt, + a::AbstractMatrix{$T}, + lda::BlasInt, + ipiv::AbstractVector{BlasInt}, + b::AbstractVecOrMat{$T}, + ldb::BlasInt, + info::Base.RefValue{BlasInt}, + ) + ccall((@blasfunc($getrs), libblastrampoline), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{BlasInt}, + Ptr{$T}, Ref{BlasInt}, Ref{BlasInt}, Clong), + trans, n, nrhs, a, lda, ipiv, b, ldb, info, 1) + return info[] end end end \ No newline at end of file diff --git a/src/linear_algebra/mumps.jl b/src/linear_algebra/mumps.jl index f8a65205..81d065cf 100644 --- a/src/linear_algebra/mumps.jl +++ b/src/linear_algebra/mumps.jl @@ -9,6 +9,7 @@ mutable struct PenaltyMUMPSWorkspace{ H::K2 x::V _ipiv::VI # For CompactBFGS LU factorization + _info::Base.RefValue{BlasInt} # For CompactBFGS LU factorization σ::T n::Int m::Int @@ -86,7 +87,7 @@ function construct_mumps_workspace( S.rhs = pointer(x) S._y_gc_haven = x - return PenaltyMUMPSWorkspace(S, H, x, nothing, zero(T), n, m, :uninitialized, false, 0) + return PenaltyMUMPSWorkspace(S, H, x, nothing, Ref{BlasInt}(0), zero(T), n, m, :uninitialized, false, 0) end function construct_mumps_workspace( @@ -147,7 +148,7 @@ function construct_mumps_workspace( S.rhs = pointer(x) S._y_gc_haven = x - return PenaltyMUMPSWorkspace(S, H, x, Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), zero(T), n, m, :uninitialized, false, 0) + return PenaltyMUMPSWorkspace(S, H, x, Vector{LinearAlgebra.BlasInt}(undef, 2 * H.B._mem), Ref{BlasInt}(0), zero(T), n, m, :uninitialized, false, 0) end function update_workspace!( @@ -345,83 +346,98 @@ function solve_system!( update_pivtol!(workspace) - # Step 3: Compute - # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] - # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] - @views mul!(y1[1:p], Uk', x1[1:n]) - @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) + if p > 0 + # Step 3: Compute + # y₁ = Fᵀx₁ = [Uᵀx₁(1:n)] + # y₁ = Fᵀx₁ = [Vᵀx₁(1:n)] + @views mul!(y1[1:p], Uk', x1[1:n]) + @views mul!(y1[(p+1):(2*p)], Vk', x1[1:n]) - # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) - # ( [A -αI] ) - # Step 4.1: Compute - # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] - # Z₁ = [A -αI] E = [A -αI] [ 0 0] - Z1 .= 0 - @views Z1[1:n, 1:p] .= Uk .* (-1) - @views Z1[1:n, (p+1):(2*p)] .= Vk + # Step 4: Assemble Schur complement (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E ) + # ( [A -αI] ) + # Step 4.1: Compute + # Z₁ = [σI+ξI Aᵀ]⁻¹ E = [σI+ξI Aᵀ]⁻¹[-U V] + # Z₁ = [A -αI] E = [A -αI] [ 0 0] + Z1 .= 0 - MUMPS.associate_rhs!(mumps, Z1; unsafe = true) - MUMPS.mumps_solve!(Z1, mumps; rhs_changed = true) + @views Z1[1:n, 1:p] .= Uk .* (-1) + @views Z1[1:n, (p+1):(2*p)] .= Vk - # MUMPS infog(1): a negative value is an error in the factorization. - if any(isnan, Z1) || mumps.infog[1] < 0 - workspace.status = :failed - return - end + MUMPS.associate_rhs!(mumps, Z1; unsafe = true) + MUMPS.mumps_solve!(Z1, mumps; rhs_changed = true) + + # MUMPS infog(1): a negative value is an error in the factorization. + if any(isnan, Z1) || mumps.infog[1] < 0 + workspace.status = :failed + return + end update_pivtol!(workspace) - # Step 4.2: Compute - # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] - # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] - Z2 .= 0 - @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) - @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) - - # Step 4.3: Compute - # Z₂ = I + Z₂ - for i = 1:(2*p) - Z2[i, i] += 1 - end + # Step 4.2: Compute + # Z₂ = FᵀZ₁ = UᵀZ₁[1:n] + # Z₂ = FᵀZ₁ = VᵀZ₁[1:n] + Z2 .= 0 + @views mul!(Z2[1:p, 1:(2*p)], Uk', Z1[1:n, (1:(2*p))]) + @views mul!(Z2[(p+1):(2*p), 1:(2*p)], Vk', Z1[1:n, (1:(2*p))]) + + # Step 4.3: Compute + # Z₂ = I + Z₂ + for i = 1:(2*p) + Z2[i, i] += 1 + end - # Step 5: Solve - # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] - # ( [A -αI] ) [y₁] - # using Julia LinearALgebra's lu! - @views LinearAlgebra.LAPACK.getrf!(Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)]) - @views LinearAlgebra.LAPACK.getrs!('N', Z2[1:(2*p), 1:(2*p)], workspace._ipiv[1:(2*p)], y1[1:(2*p)]) - if any(isnan, @view y1[1:(2*p)]) - workspace.status = :failed - return - end + # Step 5: Solve + # (I + Fᵀ [σI+ξI Aᵀ]⁻¹ E )⁻¹[y₁] + # ( [A -αI] ) [y₁] + # using LAPACK + @views info_f = getrf!(BlasInt(2p), BlasInt(2p), Z2, stride(Z2, 2), workspace._ipiv, workspace._info) + if info_f != 0 + workspace.status = :failed + return + end - # Step 6: Compute - # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] - # x₂ = E[y₂] = [ 0 0][y₂] = [0] - @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) - @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) + @views info_s = getrs!('N', BlasInt(2p), BlasInt(1), Z2[1:(2p), 1:(2p)], stride(Z2, 2), + workspace._ipiv, y1[1:(2p)], BlasInt(2p), workspace._info) + if info_s != 0 + workspace.status = :failed + return + end + if any(isnan, @view y1[1:(2*p)]) + workspace.status = :failed + return + end - # Step 7: Solve - # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] - # [x₃] = [A -αI] [x₂] - x3 .= x2 - MUMPS.associate_rhs!(mumps, x3; unsafe = true) - MUMPS.mumps_solve!(x3, mumps; rhs_changed = true) + # Step 6: Compute + # x₂ = E[y₂] = [-U V][y₂] = [-Uy₂ + Vy₂] + # x₂ = E[y₂] = [ 0 0][y₂] = [0] + @views mul!(x2[1:n], Vk, y1[(p+1):(2*p)]) + @views mul!(x2[1:n], Uk, y1[1:p], -one(eltype(y1)), one(eltype(y1))) - # MUMPS infog(1): a negative value is an error in the factorization. - if any(isnan, x3) || mumps.infog[1] < 0 - workspace.status = :failed - return - end + # Step 7: Solve + # [x₃] = [σI+ξI Aᵀ]⁻¹[x₂] + # [x₃] = [A -αI] [x₂] + x3 .= x2 + MUMPS.associate_rhs!(mumps, x3; unsafe = true) + MUMPS.mumps_solve!(x3, mumps; rhs_changed = true) + + # MUMPS infog(1): a negative value is an error in the factorization. + if any(isnan, x3) || mumps.infog[1] < 0 + workspace.status = :failed + return + end update_pivtol!(workspace) - # Step 8: - # [B Aᵀ]⁻¹[u] = x₁ - x₃ - # [A -αI] [u] = x₁ - x₃ - workspace.x .= x1 .- x3 + # Step 8: + # [B Aᵀ]⁻¹[u] = x₁ - x₃ + # [A -αI] [u] = x₁ - x₃ + workspace.x .= x1 .- x3 + else + workspace.x .= x1 + end end function get_solution!(x::V, workspace::PenaltyMUMPSWorkspace) where {V<:AbstractVector} From 50cc052b34349d642043a85059564940ee5d76da Mon Sep 17 00:00:00 2001 From: MaxenceGollier Date: Tue, 1 Sep 2026 12:10:27 -0400 Subject: [PATCH 7/7] use blas in LDLFactorizations --- ext/PeneloptLDLFactorizationsExt.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/ext/PeneloptLDLFactorizationsExt.jl b/ext/PeneloptLDLFactorizationsExt.jl index c30ab3ad..d8990bb8 100644 --- a/ext/PeneloptLDLFactorizationsExt.jl +++ b/ext/PeneloptLDLFactorizationsExt.jl @@ -5,9 +5,14 @@ using Penelopt using LinearAlgebra, SparseArrays, SparseMatricesCOO +# Import BLAS functions +import LinearAlgebra.BLAS: @blasfunc +import LinearAlgebra: BlasInt, libblastrampoline + import Penelopt: AbstractLDLTWorkspace import Penelopt: construct_ldlt_workspace, solve_system!, update_workspace! import Penelopt: get_inertia, get_solution!, get_status +import Penelopt: getrf!, getrs! import Penelopt: set_dual_inertia!, set_primal_inertia! include("LDLFactorizations/ldlt.jl")