From 0ed75b520583315c2f306b1aabfe08ae2bcc52f9 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Tue, 3 Feb 2026 09:38:07 -0500 Subject: [PATCH 1/3] Add efficient O(nnz) isdiag for sparse matrices MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add a specialized `isdiag` method for `AbstractSparseMatrixCSC` that traverses the CSC structure directly in O(nnz) time, compared to the generic fallback which is O(n²) for an n×n matrix. This addresses performance issues in packages like OrdinaryDiffEq.jl where `isdiag` checks on large sparse mass matrices (e.g., 259k variables in DAE systems) caused initialization times of ~60 seconds. With this optimization, initialization drops to <1ms. Reference: https://github.com/SciML/OrdinaryDiffEq.jl/pull/3011 Co-Authored-By: Chris Rackauckas --- src/SparseArrays.jl | 2 +- src/sparsematrix.jl | 36 ++++++++++++++++++++++++++++++++++++ test/sparsematrix_ops.jl | 28 ++++++++++++++++++++++++++++ 3 files changed, 65 insertions(+), 1 deletion(-) diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 9c0bfd19..2c8038a8 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -17,7 +17,7 @@ using LinearAlgebra: AdjOrTrans, AdjointFactorization, TransposeFactorization, m import Base: +, -, *, \, /, ==, zero import Base: Matrix, Vector import LinearAlgebra: mul!, ldiv!, rdiv!, cholesky, adjoint!, diag, eigen, dot, - issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded, + issymmetric, istril, istriu, lu, tr, transpose!, tril!, triu!, isbanded, isdiag, cond, diagm, factorize, ishermitian, norm, opnorm, lmul!, rmul!, tril, triu, matprod_dest, generic_matvecmul!, generic_matmatmul!, copytrito! diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index ae7292c7..95ab39be 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -4200,6 +4200,42 @@ function istril(A::AbstractSparseMatrixCSC, k::Integer=0) return true end +""" + isdiag(A::AbstractSparseMatrixCSC) + +Test whether a sparse matrix is diagonal. Returns `true` if all non-zero elements +are on the main diagonal. + +This implementation is O(nnz) by traversing the CSC structure directly, compared +to the generic fallback which is O(n²) for an n×n matrix. + +# Examples +```jldoctest +julia> using SparseArrays + +julia> isdiag(sparse([1, 2, 3], [1, 2, 3], [1, 2, 3])) +true + +julia> isdiag(sparse([1, 2], [1, 3], [1, 2])) +false +``` +""" +function isdiag(A::AbstractSparseMatrixCSC) + m, n = size(A) + m != n && return false + colptr = getcolptr(A) + rowval = rowvals(A) + nzval = nonzeros(A) + @inbounds for col in 1:n + for k in colptr[col]:(colptr[col + 1] - 1) + if rowval[k] != col && _isnotzero(nzval[k]) + return false + end + end + end + return true +end + _nnz(v::AbstractSparseVector) = nnz(v) _nnz(v::AbstractVector) = length(v) diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl index d99418a8..361c05fa 100644 --- a/test/sparsematrix_ops.jl +++ b/test/sparsematrix_ops.jl @@ -639,4 +639,32 @@ end end end +@testset "isdiag" begin + # Diagonal matrices should return true + @test isdiag(sparse(Diagonal(1:4))) + @test isdiag(sparse(Diagonal([1.0, 2.0, 3.0]))) + @test isdiag(spzeros(5, 5)) # Empty matrix is diagonal + + # Non-diagonal matrices should return false + @test !isdiag(sparse(Tridiagonal(1:3, 1:4, 1:3))) + @test !isdiag(sparse(Bidiagonal(1:4, 1:3, :U))) + @test !isdiag(sparse(Bidiagonal(1:4, 1:3, :L))) + @test !isdiag(sparse([1 2; 3 4])) + + # Non-square matrices should return false + @test !isdiag(sparse([1 0 0; 0 2 0])) + @test !isdiag(sparse([1 0; 0 2; 0 0])) + + # Consistency with dense isdiag + for T in Any[Diagonal(1:4), Tridiagonal(1:3, 1:4, 1:3), + Bidiagonal(1:4, 1:3, :U), diagm(-1=>1:3, 1=>1:3)] + S = sparse(T) + @test isdiag(S) == isdiag(T) + end + + # Explicit zeros on off-diagonal should still be diagonal + S = sparse([1, 2, 1], [1, 2, 2], [1.0, 2.0, 0.0]) + @test isdiag(S) +end + end # module From 0d2a9c68b70915e4169e9ba870a46fa66fdd197b Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Tue, 3 Feb 2026 14:44:40 -0500 Subject: [PATCH 2/3] Fix isdiag to support non-square matrices The generic isdiag implementation in LinearAlgebra returns true for rectangular diagonal matrices (e.g., `diagm(4,5,0 => ones(4))`). This commit removes the square matrix check to be consistent with the generic implementation. Updated tests to verify: - Non-square diagonal matrices return true - Non-square non-diagonal matrices return false Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.5 --- src/sparsematrix.jl | 1 - test/sparsematrix_ops.jl | 12 +++++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 95ab39be..0f6187eb 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -4222,7 +4222,6 @@ false """ function isdiag(A::AbstractSparseMatrixCSC) m, n = size(A) - m != n && return false colptr = getcolptr(A) rowval = rowvals(A) nzval = nonzeros(A) diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl index 361c05fa..ab2065cb 100644 --- a/test/sparsematrix_ops.jl +++ b/test/sparsematrix_ops.jl @@ -651,9 +651,15 @@ end @test !isdiag(sparse(Bidiagonal(1:4, 1:3, :L))) @test !isdiag(sparse([1 2; 3 4])) - # Non-square matrices should return false - @test !isdiag(sparse([1 0 0; 0 2 0])) - @test !isdiag(sparse([1 0; 0 2; 0 0])) + # Non-square diagonal matrices should return true (consistent with generic isdiag) + @test isdiag(sparse([1 0 0; 0 2 0])) # 2x3 diagonal + @test isdiag(sparse([1 0; 0 2; 0 0])) # 3x2 diagonal + @test isdiag(spzeros(3, 5)) # Empty non-square matrix is diagonal + @test isdiag(spzeros(5, 3)) + + # Non-square non-diagonal matrices should return false + @test !isdiag(sparse([1 1 0; 0 2 0])) # Off-diagonal element + @test !isdiag(sparse([1 0; 0 2; 1 0])) # Off-diagonal element # Consistency with dense isdiag for T in Any[Diagonal(1:4), Tridiagonal(1:3, 1:4, 1:3), From ea96eddbd05ecc032a15b5495ab960120cd65d1d Mon Sep 17 00:00:00 2001 From: Chris Rackauckas - Beep Boop Edition Date: Tue, 3 Feb 2026 15:28:03 -0500 Subject: [PATCH 3/3] Apply suggestion from @ChrisRackauckas Co-authored-by: Christopher Rackauckas --- src/sparsematrix.jl | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl index 0f6187eb..c799fe75 100644 --- a/src/sparsematrix.jl +++ b/src/sparsematrix.jl @@ -4200,26 +4200,6 @@ function istril(A::AbstractSparseMatrixCSC, k::Integer=0) return true end -""" - isdiag(A::AbstractSparseMatrixCSC) - -Test whether a sparse matrix is diagonal. Returns `true` if all non-zero elements -are on the main diagonal. - -This implementation is O(nnz) by traversing the CSC structure directly, compared -to the generic fallback which is O(n²) for an n×n matrix. - -# Examples -```jldoctest -julia> using SparseArrays - -julia> isdiag(sparse([1, 2, 3], [1, 2, 3], [1, 2, 3])) -true - -julia> isdiag(sparse([1, 2], [1, 3], [1, 2])) -false -``` -""" function isdiag(A::AbstractSparseMatrixCSC) m, n = size(A) colptr = getcolptr(A)