Rename cond(LU) to the helper function _cond1Inf and allow for a norm

argument. The old version was broken and inefficient. While it would
be possible to fix just the bug, the performance problem is a consequence
of how condition number estimation works, i.e. a norm estimate is required
and that is not natural to produce from the factorized matrix. Hence,
the change to a helper function.

Fixes #21453

base/linalg/dense.jl

@@ -915,10 +915,12 @@ function cond(A::AbstractMatrix, p::Real=2)
         return maxv == 0.0 ? oftype(real(A[1,1]),Inf) : maxv / minimum(v)
     elseif p == 1 || p == Inf
         checksquare(A)
-        return cond(lufact(A), p)
+        return _cond1Inf(A, p)
     end
     throw(ArgumentError("p-norm must be 1, 2 or Inf, got $p"))
 end
+_cond1Inf(A::StridedMatrix{<:BlasFloat}, p::Real) = _cond1Inf(lufact(A), p, norm(A, p))
+_cond1Inf(A::AbstractMatrix, p::Real)             = norm(A, p)*norm(inv(A), p)
 
 ## Lyapunov and Sylvester equation
 

base/linalg/lu.jl

@@ -112,7 +112,6 @@ The relationship between `F` and `A` is
 |:---------------------------------|:-----|:-----------------------|
 | [`/`](@ref)                      | ✓    |                        |
 | [`\\`](@ref)                     | ✓    | ✓                      |
-| [`cond`](@ref)                   | ✓    |                        |
 | [`det`](@ref)                    | ✓    | ✓                      |
 | [`logdet`](@ref)                 | ✓    | ✓                      |
 | [`logabsdet`](@ref)              | ✓    | ✓                      |
@@ -309,9 +308,12 @@ inv!(A::LU{<:BlasFloat,<:StridedMatrix}) =
 inv(A::LU{<:BlasFloat,<:StridedMatrix}) =
     inv!(LU(copy(A.factors), copy(A.ipiv), copy(A.info)))
 
-cond(A::LU{<:BlasFloat,<:StridedMatrix}, p::Number) =
-    inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm((A[:L]*A[:U])[A[:p],:], p)))
-cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
+function _cond(A::LU{<:BlasFloat,<:StridedMatrix}, p::Number, normA::Real)
+    if p != 1 && p != Inf
+        throw(ArgumentError("p must be either 1 or Inf"))
+    end
+    return inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, normA))
+end
 
 # Tridiagonal
 

test/linalg/lu.jl

@@ -208,3 +208,6 @@ end
 
 @test @inferred(logdet(Complex64[1.0f0 0.5f0; 0.5f0 -1.0f0])) === 0.22314355f0 + 3.1415927f0im
 @test_throws DomainError logdet([1 1; 1 -1])
+
+# Issue 21453.
+@test_throws ArgumentError LinAlg._cond1Inf(lufact(randn(5,5)), 2, 2.0)