From c347c8ffde1c4a8e893456bb766c6b9b0869fe04 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Wed, 29 Apr 2026 17:48:44 +0200 Subject: [PATCH 01/10] Review return types from solves --- src/linalg.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/linalg.jl b/src/linalg.jl index 3115ff79..dc306683 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1309,7 +1309,11 @@ function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, :: C end -function (\)(A::Union{UpperTriangular,LowerTriangular}, B::AbstractSparseMatrixCSC) +# the following matrix types have typically dense inverses, even when they wrap sparse matrices +# so they should return dense arrays +const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} + +function (\)(A::StructuredWithDenseInverse, B::AbstractSparseMatrixCSC) require_one_based_indexing(B) TAB = promote_op(\, eltype(A), eltype(B)) ldiv!(Matrix{TAB}(undef, size(B)), A, B) @@ -1319,6 +1323,17 @@ function (\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSpars TAB = LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B)) ldiv!(Matrix{TAB}(undef, size(B)), A, B) end +function (/)(A::AbstractSparseMatrixCSC, B::StructuredWithDenseInverse) + require_one_based_indexing(A) + TAB = promote_op(/, eltype(A), eltype(B)) + LinearAlgebra._rdiv!(Matrix{TAB}(undef, size(A)), A, B) +end +function (/)(A::AbstractSparseMatrixCSC, B::Union{UnitUpperTriangular,UnitLowerTriangular}) + require_one_based_indexing(A) + TAB = LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B)) + LinearAlgebra._rdiv!(Matrix{TAB}(undef, size(A)), A, B) +end + # (*)(L::DenseTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) ## end of triangular From d279cab5109612341d3e4e3cb1159e7072e7eff7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Thu, 30 Apr 2026 11:43:56 +0200 Subject: [PATCH 02/10] only define destination array type --- src/linalg.jl | 31 ++++++++----------------------- 1 file changed, 8 insertions(+), 23 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index dc306683..035d47db 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1313,29 +1313,14 @@ end # so they should return dense arrays const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} -function (\)(A::StructuredWithDenseInverse, B::AbstractSparseMatrixCSC) - require_one_based_indexing(B) - TAB = promote_op(\, eltype(A), eltype(B)) - ldiv!(Matrix{TAB}(undef, size(B)), A, B) -end -function (\)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::AbstractSparseMatrixCSC) - require_one_based_indexing(B) - TAB = LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B)) - ldiv!(Matrix{TAB}(undef, size(B)), A, B) -end -function (/)(A::AbstractSparseMatrixCSC, B::StructuredWithDenseInverse) - require_one_based_indexing(A) - TAB = promote_op(/, eltype(A), eltype(B)) - LinearAlgebra._rdiv!(Matrix{TAB}(undef, size(A)), A, B) -end -function (/)(A::AbstractSparseMatrixCSC, B::Union{UnitUpperTriangular,UnitLowerTriangular}) - require_one_based_indexing(A) - TAB = LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B)) - LinearAlgebra._rdiv!(Matrix{TAB}(undef, size(A)), A, B) -end - -# (*)(L::DenseTriangular, B::AbstractSparseMatrixCSC) = lmul!(L, Array(B)) - +matldiv_dest(A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = + Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) +matldiv_dest(A::LinearAlgebra.UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = + Matrix{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, size(B)) +matrdiv_dest(A::QuasiSparseMatrix, B::StructuredWithDenseInverse) = + Matrix{promote_op(/, eltype(A), eltype(B))}(undef, size(A)) +matrdiv_dest(A::QuasiSparseMatrix, B::LinearAlgebra.UnitUpperOrUnitLowerTriangular) = + Matrix{LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B))}(undef, size(A)) ## end of triangular # symmetric/Hermitian From a6b1e7312e8cd13bd30a7de859400579b4c69a79 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 8 May 2026 12:03:15 +0200 Subject: [PATCH 03/10] update --- src/linalg.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 035d47db..8d10cd1b 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1313,13 +1313,13 @@ end # so they should return dense arrays const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} -matldiv_dest(A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = +matop_dest(::typeof(\), A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) -matldiv_dest(A::LinearAlgebra.UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = +matop_dest(::typeof(\), A::LinearAlgebra.UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = Matrix{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, size(B)) -matrdiv_dest(A::QuasiSparseMatrix, B::StructuredWithDenseInverse) = +matop_dest(::typeof(/), A::QuasiSparseMatrix, B::StructuredWithDenseInverse) = Matrix{promote_op(/, eltype(A), eltype(B))}(undef, size(A)) -matrdiv_dest(A::QuasiSparseMatrix, B::LinearAlgebra.UnitUpperOrUnitLowerTriangular) = +matop_dest(::typeof(/), A::QuasiSparseMatrix, B::LinearAlgebra.UnitUpperOrUnitLowerTriangular) = Matrix{LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B))}(undef, size(A)) ## end of triangular From a2a92874f6dc55a576a1d80e683f2e7ea3e02ac2 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 8 May 2026 16:52:58 +0200 Subject: [PATCH 04/10] include vector solves --- src/SparseArrays.jl | 2 +- src/linalg.jl | 8 ++++++-- src/sparsevector.jl | 11 ----------- 3 files changed, 7 insertions(+), 14 deletions(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 6c0e83b4..135fbbb9 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -11,7 +11,7 @@ using Base.Order: Forward using LinearAlgebra using LinearAlgebra: AdjOrTrans, AdjointFactorization, TransposeFactorization, matprod, AbstractQ, AdjointQ, HessenbergQ, QRCompactWYQ, QRPackedQ, LQPackedQ, MulAddMul, - UpperOrLowerTriangular, @stable_muladdmul + UpperOrLowerTriangular, UnitUpperOrUnitLowerTriangular, @stable_muladdmul import Base: +, -, *, \, /, ==, zero diff --git a/src/linalg.jl b/src/linalg.jl index 8d10cd1b..4b5aa708 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1313,13 +1313,17 @@ end # so they should return dense arrays const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} +matop_dest(::typeof(\), A::StructuredWithDenseInverse, b::AbstractSparseVector) = + Vector{promote_op(\, eltype(A), eltype(B))}(undef, length(b)) +matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, b::AbstractSparseVector) = + Vector{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, length(b)) matop_dest(::typeof(\), A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) -matop_dest(::typeof(\), A::LinearAlgebra.UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = +matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = Matrix{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, size(B)) matop_dest(::typeof(/), A::QuasiSparseMatrix, B::StructuredWithDenseInverse) = Matrix{promote_op(/, eltype(A), eltype(B))}(undef, size(A)) -matop_dest(::typeof(/), A::QuasiSparseMatrix, B::LinearAlgebra.UnitUpperOrUnitLowerTriangular) = +matop_dest(::typeof(/), A::QuasiSparseMatrix, B::UnitUpperOrUnitLowerTriangular) = Matrix{LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B))}(undef, size(A)) ## end of triangular diff --git a/src/sparsevector.jl b/src/sparsevector.jl index 1d7daf0b..8151ec48 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -2160,17 +2160,6 @@ for isunittri in (true, false), islowertri in (true, false) halfstr = islowertri ? "Lower" : "Upper" tritype = :(LinearAlgebra.$(Symbol(unitstr, halfstr, "Triangular"))) - # build out-of-place left-division operations - # broad method where elements are Numbers - @eval function \(A::$tritype{<:TA,<:AbstractMatrix}, b::AbstractCompressedVector{Tb}) where {TA<:Number,Tb<:Number} - TAb = $(isunittri ? - :(typeof(zero(TA)*zero(Tb) + zero(TA)*zero(Tb))) : - :(typeof((zero(TA)*zero(Tb) + zero(TA)*zero(Tb))/one(TA))) ) - return LinearAlgebra.ldiv!(convert(AbstractArray{TAb}, A), convert(Array{TAb}, b)) - end - # fallback where elements are not Numbers - @eval \(A::$tritype, b::AbstractCompressedVector) = LinearAlgebra.ldiv!(A, copy(b)) - # faster method requiring good view support of the # triangular matrix type. hence the StridedMatrix restriction. for (istrans, applyxform, xformtype, xformop) in ( From fb46a43caf00d338eaea12acbfba5868991ee2cc Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 16 May 2026 16:08:58 +0200 Subject: [PATCH 05/10] fix typos --- src/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index 0d30a476..decb1a02 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1314,9 +1314,9 @@ end const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} matop_dest(::typeof(\), A::StructuredWithDenseInverse, b::AbstractSparseVector) = - Vector{promote_op(\, eltype(A), eltype(B))}(undef, length(b)) + Vector{promote_op(\, eltype(A), eltype(b))}(undef, length(b)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, b::AbstractSparseVector) = - Vector{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, length(b)) + Vector{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(b))}(undef, length(b)) matop_dest(::typeof(\), A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = From 6f6052d1348faab5b0ef6bdea000f4b963cc463d Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 16 May 2026 17:08:48 +0200 Subject: [PATCH 06/10] add tests --- test/linalg.jl | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/test/linalg.jl b/test/linalg.jl index c50de457..3700a406 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -147,6 +147,42 @@ end end end +@testset "destination array density in solves" begin + O = diagm(-1 => fill(-1, 9), 0 => fill(2, 10), 1 => fill(-1, 9)) + wrappers = (a -> Bidiagonal(a, :U), + a -> Bidiagonal(a, :L), + a -> SymTridiagonal(2diag(a), -diag(a, 1)), + a -> Tridiagonal(-diag(a, -1), 2diag(a), -diag(a, 1)), + LowerTriangular, + UnitLowerTriangular, + UpperTriangular, + UnitUpperTriangular, + UpperHessenberg) + for T in wrappers + A = T(O) + bs = sprandn(10, 0.3) + bd = Array(bs) + x = A \ bs + @test x ≈ A \ bd + @test !issparse(x) + Bs = sprandn(10, 3, 0.2) + Bd = Matrix(Bs) + X = A \ Bs + @test X ≈ A \ Bd + @test !issparse(X) + Cs = copy(Bs') + Cd = Matrix(Cs) + Y = Cs / A + @test Y ≈ Cd / A + @test !issparse(Y) + end + b, B = ones(Int, 10), ones(Int, 10, 10) + for T in (UnitLowerTriangular, UnitUpperTriangular) + A = T(O) + @test eltype(A \ b) == eltype(A \ B) == eltype(B / A) == Int + end +end + @testset "multiplication of special sparse with dense matrix" begin # this results in a call of the most generic multiplication code in LinearAlgebra.jl A = randn(2, 2) @@ -228,6 +264,7 @@ begin AW = tr(wr(A)) MAW = tr(wr(MA)) @test AW \ B ≈ MAW \ B + @test !issparse(AW \ B) # and for SparseMatrixCSCView - a view of all rows and unit range of cols vAW = tr(wr(view([zero(A)+I A], :, (n+1):2n))) @test vAW \ B ≈ AW \ B From d91ef7a1570b74225addbd0736bcdb84476b37d1 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 17 May 2026 17:43:13 +0200 Subject: [PATCH 07/10] fix tests, choose dense arrays for any solves with sparse arrays --- src/linalg.jl | 10 +++------- test/linalg.jl | 13 +++++++------ 2 files changed, 10 insertions(+), 13 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index decb1a02..f5592f24 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1309,19 +1309,15 @@ function LinearAlgebra.generic_trimatdiv!(C::StridedVecOrMat, uploc, isunitc, :: C end -# the following matrix types have typically dense inverses, even when they wrap sparse matrices -# so they should return dense arrays -const StructuredWithDenseInverse = Union{Bidiagonal,SymTridiagonal,Tridiagonal,LowerTriangular,UpperTriangular,UpperHessenberg} - -matop_dest(::typeof(\), A::StructuredWithDenseInverse, b::AbstractSparseVector) = +matop_dest(::typeof(\), A, b::AbstractSparseVector) = Vector{promote_op(\, eltype(A), eltype(b))}(undef, length(b)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, b::AbstractSparseVector) = Vector{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(b))}(undef, length(b)) -matop_dest(::typeof(\), A::StructuredWithDenseInverse, B::QuasiSparseMatrix) = +matop_dest(::typeof(\), A, B::QuasiSparseMatrix) = Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = Matrix{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, size(B)) -matop_dest(::typeof(/), A::QuasiSparseMatrix, B::StructuredWithDenseInverse) = +matop_dest(::typeof(/), A::QuasiSparseMatrix, B) = Matrix{promote_op(/, eltype(A), eltype(B))}(undef, size(A)) matop_dest(::typeof(/), A::QuasiSparseMatrix, B::UnitUpperOrUnitLowerTriangular) = Matrix{LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B))}(undef, size(A)) diff --git a/test/linalg.jl b/test/linalg.jl index 3700a406..715c166d 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -151,30 +151,31 @@ end O = diagm(-1 => fill(-1, 9), 0 => fill(2, 10), 1 => fill(-1, 9)) wrappers = (a -> Bidiagonal(a, :U), a -> Bidiagonal(a, :L), - a -> SymTridiagonal(2diag(a), -diag(a, 1)), - a -> Tridiagonal(-diag(a, -1), 2diag(a), -diag(a, 1)), + # SymTridiagonal, + Tridiagonal, LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, - UpperHessenberg) + # UpperHessenberg + ) for T in wrappers A = T(O) bs = sprandn(10, 0.3) bd = Array(bs) x = A \ bs @test x ≈ A \ bd - @test !issparse(x) + @test !issparse(x) broken=(T in (SymTridiagonal, UpperHessenberg)) Bs = sprandn(10, 3, 0.2) Bd = Matrix(Bs) X = A \ Bs @test X ≈ A \ Bd - @test !issparse(X) + @test !issparse(X) broken=(T in (SymTridiagonal, UpperHessenberg)) Cs = copy(Bs') Cd = Matrix(Cs) Y = Cs / A @test Y ≈ Cd / A - @test !issparse(Y) + @test !issparse(Y) broken=(T in (SymTridiagonal, UpperHessenberg)) end b, B = ones(Int, 10), ones(Int, 10, 10) for T in (UnitLowerTriangular, UnitUpperTriangular) From 3f35cc529e6c14d681125ad80b9dbf5722a2a41f Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 17 May 2026 20:08:28 +0200 Subject: [PATCH 08/10] fix ambiguities --- src/linalg.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/linalg.jl b/src/linalg.jl index f5592f24..69e2a0d5 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -1313,14 +1313,20 @@ matop_dest(::typeof(\), A, b::AbstractSparseVector) = Vector{promote_op(\, eltype(A), eltype(b))}(undef, length(b)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, b::AbstractSparseVector) = Vector{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(b))}(undef, length(b)) +matop_dest(::typeof(\), A::Diagonal, b::AbstractSparseVector) = + similar(b , promote_op(\, eltype(A), eltype(b))) matop_dest(::typeof(\), A, B::QuasiSparseMatrix) = Matrix{promote_op(\, eltype(A), eltype(B))}(undef, size(B)) +matop_dest(::typeof(\), A::Diagonal, B::QuasiSparseMatrix) = + similar(B , promote_op(\, eltype(A), eltype(B)), size(B)) matop_dest(::typeof(\), A::UnitUpperOrUnitLowerTriangular, B::QuasiSparseMatrix) = Matrix{LinearAlgebra._inner_type_promotion(\, eltype(A), eltype(B))}(undef, size(B)) matop_dest(::typeof(/), A::QuasiSparseMatrix, B) = Matrix{promote_op(/, eltype(A), eltype(B))}(undef, size(A)) matop_dest(::typeof(/), A::QuasiSparseMatrix, B::UnitUpperOrUnitLowerTriangular) = Matrix{LinearAlgebra._inner_type_promotion(/, eltype(A), eltype(B))}(undef, size(A)) +matop_dest(::typeof(/), A::QuasiSparseMatrix, B::Diagonal) = + similar(A , promote_op(/, eltype(A), eltype(B)), size(A)) ## end of triangular # symmetric/Hermitian From 0be6de56a8e8ef334750e316aef8cccef2ec8db7 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 18 May 2026 13:01:01 +0200 Subject: [PATCH 09/10] extend tests --- test/linalg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index 715c166d..d587c610 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -151,13 +151,13 @@ end O = diagm(-1 => fill(-1, 9), 0 => fill(2, 10), 1 => fill(-1, 9)) wrappers = (a -> Bidiagonal(a, :U), a -> Bidiagonal(a, :L), - # SymTridiagonal, + SymTridiagonal, Tridiagonal, LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, - # UpperHessenberg + UpperHessenberg ) for T in wrappers A = T(O) From 1cec95e910082ae2f0044d493e3dde956f4fe5ba Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 18 May 2026 18:09:01 +0200 Subject: [PATCH 10/10] fix tests --- test/linalg.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/linalg.jl b/test/linalg.jl index d587c610..758ed3cb 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -157,7 +157,8 @@ end UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, - UpperHessenberg + # UpperHessenberg, + a -> UpperHessenberg(float(a)) ) for T in wrappers A = T(O) @@ -165,17 +166,17 @@ end bd = Array(bs) x = A \ bs @test x ≈ A \ bd - @test !issparse(x) broken=(T in (SymTridiagonal, UpperHessenberg)) + @test !issparse(x) Bs = sprandn(10, 3, 0.2) Bd = Matrix(Bs) X = A \ Bs @test X ≈ A \ Bd - @test !issparse(X) broken=(T in (SymTridiagonal, UpperHessenberg)) + @test !issparse(X) Cs = copy(Bs') Cd = Matrix(Cs) Y = Cs / A @test Y ≈ Cd / A - @test !issparse(Y) broken=(T in (SymTridiagonal, UpperHessenberg)) + @test !issparse(Y) end b, B = ones(Int, 10), ones(Int, 10, 10) for T in (UnitLowerTriangular, UnitUpperTriangular)