From ad0c56ec23d887e4b55db41131882be73a37a6c7 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreasnoackjensen@gmail.com>
Date: Tue, 25 Apr 2017 12:59:16 -0400
Subject: [PATCH] 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 |  4 +++-
 base/linalg/lu.jl    | 10 ++++++----
 test/linalg/lu.jl    |  3 +++
 3 files changed, 12 insertions(+), 5 deletions(-)

diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index 7beccb480789ac..476c67fe28bff1 100644
--- a/base/linalg/dense.jl
+++ b/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
 
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index 6f715ea76474d7..ba3422e18ce74d 100644
--- a/base/linalg/lu.jl
+++ b/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
 
diff --git a/test/linalg/lu.jl b/test/linalg/lu.jl
index 674a5e71cf6367..fd6a59621a9c0b 100644
--- a/test/linalg/lu.jl
+++ b/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)