From 2a16e449df4a229b74db14187113a38d53e64bfd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 1 Jul 2026 08:50:15 +0200 Subject: [PATCH 1/4] Faster add_sub_jprod --- src/buffer.jl | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/buffer.jl b/src/buffer.jl index 0338724..d3c5a4b 100644 --- a/src/buffer.jl +++ b/src/buffer.jl @@ -131,9 +131,22 @@ function add_sub_jprod!( Jv::AbstractVector, I, ) - # `view(cache, I)` would be terribly slow, only the number of elements of `I` matter here + # `Jv[k] += ⟨A[:, I[k]], vec(V)⟩` for each `k`. Going through + # `mul!(Jv, view(A, :, I)', vec(V))` hits the generic (non-BLAS, sparse + # `getindex`-based) matvec since the adjoint of a `SparseMatrixCSC` column + # subset has no specialized method; iterate the CSC columns directly. A = model.jprod_buffer[i.value] - return _add_jprod!(V, Jv, view(A, :, I)) + v = _vec(V) + rows = SparseArrays.rowvals(A) + vals = SparseArrays.nonzeros(A) + @inbounds for (k, j) in enumerate(I) + acc = zero(eltype(Jv)) + for p in SparseArrays.nzrange(A, j) + acc += vals[p] * v[rows[p]] + end + Jv[k] += acc + end + return Jv end function add_jprod!( From b247577fcf3b11b96deab1f0a4218e9739484425 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 27 Jul 2026 19:04:09 +0200 Subject: [PATCH 2/4] better --- src/buffer.jl | 21 ++++-------------- src/factorization.jl | 53 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 17 deletions(-) diff --git a/src/buffer.jl b/src/buffer.jl index d3c5a4b..4d04632 100644 --- a/src/buffer.jl +++ b/src/buffer.jl @@ -120,8 +120,8 @@ _vec(x::FillArrays.Zeros{T}) where {T} = FillArrays.Zeros{T}(length(x)) _vec(x::AbstractArray) = UnsafeArrays.uview(x, :) _vec(x::Base.ReshapedArray) = _vec(parent(x)) -function _add_jprod!(V, Jv::AbstractArray{T}, A) where {T} - return LinearAlgebra.mul!(Jv, A', _vec(V), true, true) +function _add_jprod!(V, Jv::AbstractArray, A) + return _add_mul!(Jv, A', _vec(V), true) end function add_sub_jprod!( @@ -131,22 +131,9 @@ function add_sub_jprod!( Jv::AbstractVector, I, ) - # `Jv[k] += ⟨A[:, I[k]], vec(V)⟩` for each `k`. Going through - # `mul!(Jv, view(A, :, I)', vec(V))` hits the generic (non-BLAS, sparse - # `getindex`-based) matvec since the adjoint of a `SparseMatrixCSC` column - # subset has no specialized method; iterate the CSC columns directly. + # `view(cache, I)` would be terribly slow, only the number of elements of `I` matter here A = model.jprod_buffer[i.value] - v = _vec(V) - rows = SparseArrays.rowvals(A) - vals = SparseArrays.nonzeros(A) - @inbounds for (k, j) in enumerate(I) - acc = zero(eltype(Jv)) - for p in SparseArrays.nzrange(A, j) - acc += vals[p] * v[rows[p]] - end - Jv[k] += acc - end - return Jv + return _add_jprod!(V, Jv, view(A, :, I)) end function add_jprod!( diff --git a/src/factorization.jl b/src/factorization.jl index a0884b0..848619b 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -457,6 +457,59 @@ function _add_mul!( end end +function _add_mul!( + res::AbstractVector, + Ft::LinearAlgebra.Adjoint{<:Any,<:SparseArrays.SparseMatrixCSC}, + C::AbstractVector, + α, +) + F = parent(Ft) + @assert axes(C, 1) == axes(F, 1) + @assert axes(res, 1) == axes(F, 2) + for col in axes(F, 2) + acc = zero(eltype(res)) + for i in SparseArrays.nzrange(F, col) + acc += SparseArrays.nonzeros(F)[i] * C[SparseArrays.rowvals(F)[i]] + end + res[col] += acc * α + end + return res +end + +# The adjoint of a column subset of a `SparseMatrixCSC` has no specialized +# `mul!` method so it would hit the generic `getindex`-based fallback +# that does a binary search for each entry of the full dense index space; +# iterate the stored entries of the selected columns of the parent instead. +function _add_mul!( + res::AbstractVector, + Ft::LinearAlgebra.Adjoint{ + <:Any, + <:SubArray{ + <:Any, + 2, + <:SparseArrays.SparseMatrixCSC, + <:Tuple{Base.Slice,Any}, + }, + }, + C::AbstractVector, + α, +) + F = parent(parent(Ft)) + cols = parentindices(parent(Ft))[2] + @assert axes(C, 1) == axes(F, 1) + @assert axes(res, 1) == axes(cols, 1) + rows = SparseArrays.rowvals(F) + vals = SparseArrays.nonzeros(F) + @inbounds for (k, col) in enumerate(cols) + acc = zero(eltype(res)) + for i in SparseArrays.nzrange(F, col) + acc += vals[i] * C[rows[i]] + end + res[k] += acc * α + end + return res +end + function _add_mul!( res::AbstractMatrix, F::SparseArrays.SparseMatrixCSC, From c99e5ff42b5e5c62b5cdf12e6d07ee0573680ad5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 27 Jul 2026 19:13:13 +0200 Subject: [PATCH 3/4] Fix --- src/factorization.jl | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 848619b..d587ff0 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -457,23 +457,18 @@ function _add_mul!( end end +# `SparseArrays` has a specialized (CSC-traversing) `mul!` method for the +# adjoint of a `SparseMatrixCSC` so we can safely delegate. It requires a +# strided destination, so we restrict `res` to `StridedVector` to get a +# `MethodError` for a non-strided destination rather than silently hitting +# the generic `getindex`-based fallback. function _add_mul!( - res::AbstractVector, + res::StridedVector, Ft::LinearAlgebra.Adjoint{<:Any,<:SparseArrays.SparseMatrixCSC}, C::AbstractVector, α, ) - F = parent(Ft) - @assert axes(C, 1) == axes(F, 1) - @assert axes(res, 1) == axes(F, 2) - for col in axes(F, 2) - acc = zero(eltype(res)) - for i in SparseArrays.nzrange(F, col) - acc += SparseArrays.nonzeros(F)[i] * C[SparseArrays.rowvals(F)[i]] - end - res[col] += acc * α - end - return res + return LinearAlgebra.mul!(res, Ft, C, α, true) end # The adjoint of a column subset of a `SparseMatrixCSC` has no specialized From 4340e963d9319060b56cc42d28149a9122dc93aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 27 Jul 2026 19:19:49 +0200 Subject: [PATCH 4/4] Fix format --- test/BurerMonteiro.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 9fb325a..67f5c01 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -134,7 +134,7 @@ end; set_attribute(model, "max_iter", 0) optimize!(model) @test termination_status(model) == MOI.ITERATION_LIMIT - nlp = unsafe_backend(model).model; + nlp = unsafe_backend(model).model @test nlp.C isa Vector{SparseMatrixCSC{T,Int}} @test nlp.C[1] == [3 1; 1 4] @test nlp.A isa Matrix{ @@ -210,7 +210,7 @@ end optimize!(model) @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) - nlp = unsafe_backend(model).model; + nlp = unsafe_backend(model).model T = Float64 MT = LRO.Factorization{T,Matrix{T},LRO.Ones{T}} @test nlp.C isa Vector{MT}