From ef013bbd709470f00549fba5e6854f99c4a54c9b Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <virals@gmail.com>
Date: Mon, 26 Oct 2026 23:08:19 +0000
Subject: [PATCH] Catch only the factorization's `SingularException` in sparse
 `cond`

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/linalg.jl          | 30 +++++++++++++++++++++---------
 test/solvers/issues.jl |  5 ++++-
 2 files changed, 25 insertions(+), 10 deletions(-)

diff --git a/src/linalg.jl b/src/linalg.jl
index caf696c3..9caf66d0 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -993,15 +993,16 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2)
     # as for a dense matrix; the estimator of the inverse's norm needs at least one column
     (p == 1 || p == Inf) && checksquare(A) == 0 && return zero(T)
     if p == 1 || p == Inf
-        # a singular matrix has an infinite condition number, as for a dense matrix; a
-        # triangular or diagonal `factorize(A)` reports it only in the solve
-        normAinv = try
-            p == 1 ? opnormestinv(A) : opnormestinv(copy(A'))
+        B = p == 1 ? A : copy(A')
+        # a singular matrix has an infinite condition number, as for a dense matrix
+        F = try
+            factorize(B)
         catch e
             e isa SingularException || rethrow()
             return convert(T, Inf)
         end
-        return convert(T, opnorm(A, p) * normAinv)
+        _issingularfactor(F) && return convert(T, Inf)
+        return convert(T, opnorm(A, p) * _opnormestinv(B, F, min(2, size(B, 1))))
     elseif p == 2
         throw(ArgumentError("only 1- and Inf-norm condition numbers are implemented for sparse matrices, for 2-norm try cond(Array(A), 2) instead"))
     else
@@ -1009,11 +1010,22 @@ function cond(A::AbstractSparseMatrixCSC, p::Real=2)
     end
 end
 
-function opnormestinv(A::AbstractSparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) where T
+# `factorize` throws for a singular LU or Hermitian factorization, but leaves a triangular
+# or diagonal one to the solve
+_issingularfactor(F) = false
+_issingularfactor(D::Diagonal) = any(iszero, D.diag)
+_issingularfactor(L::LowerTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(false, parent(L)))
+_issingularfactor(U::UpperTriangular{<:Any,<:AbstractSparseMatrixCSC}) = !iszero(_trisingular(true, parent(U)))
+
+function opnormestinv(A::AbstractSparseMatrixCSC, t::Integer = min(2,maximum(size(A))))
+    checksquare(A)
+    return _opnormestinv(A, factorize(A), t)
+end
+
+# the estimate for the factorization `F` of `A`
+function _opnormestinv(A::AbstractSparseMatrixCSC{T}, F, t::Integer) where T
     maxiter = 5
-    # Check the input
-    n = checksquare(A)
-    F = factorize(A)
+    n = size(A, 1)
     if t <= 0
         throw(ArgumentError("number of blocks must be a positive integer"))
     end
diff --git a/test/solvers/issues.jl b/test/solvers/issues.jl
index 0d50f795..2eff6dfa 100644
--- a/test/solvers/issues.jl
+++ b/test/solvers/issues.jl
@@ -256,8 +256,11 @@ end
 
 @testset "cond of a singular matrix is Inf" begin
     @test cond(sparse([1.0 2; 3 6]), 1) === Inf
-    # the triangular solve, not the factorization, finds this one singular
+    # `factorize` leaves the zero on the diagonal of a triangular or diagonal factor
     @test cond(sparse([1.0 0; 2 0]), Inf) === Inf
+    @test cond(sparse([1.0 0; 0 0]), 1) === Inf
+    # Hermitian, so `factorize` tries Cholesky and LDLt before the LU finds it singular
+    @test cond(sparse([1.0 1; 1 1]), 1) === Inf
     # the LU works in double precision, but the result is of the matrix's precision
     @test cond(sparse(Float32[1 2; 3 4]), 1) isa Float32
 end