diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 479b1470042af..c710c1fe81842 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -250,7 +250,8 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular
inv{T<:BlasFloat}(A::LU{T})=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
-cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond{T<:BlasFloat}(A::LU{T}, p::Number) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p)))
+cond(A::LU, p::Number) = norm(A[:L]*A[:U],p)*norm(inv(A),p)
####################
# QR Factorization #
diff --git a/test/linalg1.jl b/test/linalg1.jl
index ba59b8842555c..a2956c47b3ff1 100644
--- a/test/linalg1.jl
+++ b/test/linalg1.jl
@@ -20,14 +20,16 @@ areal = randn(n,n)/2
aimg = randn(n,n)/2
breal = randn(n,2)/2
bimg = randn(n,2)/2
-for eltya in (Float16, Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
- for eltyb in (Float16, Float32, Float64, Complex32, Complex64, Complex128, Int)
- a = eltya == Int ? rand(1:5, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
+for eltya in (Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
+ for eltyb in (Float32, Float64, Complex32, Complex64, Complex128, Int)
+ a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
asym = a'+a # symmetric indefinite
apd = a'*a # symmetric positive-definite
- ε = max(eps(abs(float(one(eltya)))),eps(abs(float(one(eltyb)))))
+ εa = eps(abs(float(one(eltya))))
+ εb = eps(abs(float(one(eltyb))))
+ ε = max(εa,εb)
debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
@@ -50,8 +52,8 @@ debug && println("(Automatic) upper Cholesky factor")
@test norm(apd*x-b)/norm(b) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ
@test_approx_eq apd * inv(capd) eye(n)
- @test_approx_eq_eps a*(capd\(a'*b)) b 8000ε
- @test_approx_eq det(capd) det(apd)
+ @test norm(a*(capd\(a'*b)) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
+ @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit
@test_approx_eq logdet(capd) log(det(capd)) # logdet is less likely to overflow
debug && println("lower Cholesky factor")
@@ -64,7 +66,7 @@ debug && println("pivoted Choleksy decomposition")
cpapd = cholfact(apd, pivot=true)
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
- @test_approx_eq_eps b apd * (cpapd\b) 15000ε
+ @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
if isreal(apd)
@test_approx_eq apd * inv(cpapd) eye(n)
end
@@ -83,12 +85,13 @@ debug && println("Bunch-Kaufman factors of a pos-def matrix")
end
debug && println("(Automatic) Square LU decomposition")
+ κ = cond(a,1)
lua = factorize(a)
l,u,p = lua[:L], lua[:U], lua[:p]
@test_approx_eq l*u a[p,:]
@test_approx_eq l[invperm(p),:]*u a
@test_approx_eq a * inv(lua) eye(n)
- @test_approx_eq_eps a*(lua\b) b 80ε
+ @test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
debug && println("Thin LU")
lua = lufact(a[:,1:5])
@@ -188,10 +191,11 @@ debug && println("Generalized svd")
end
debug && println("Solve square general system of equations")
+ κ = cond(a,1)
x = a \ b
- @test_approx_eq_eps a*x b 80ε
@test_throws b'\b
@test_throws b\b'
+ @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit!
debug && println("Solve upper triangular system")
x = triu(a) \ b
diff --git a/test/linalg2.jl b/test/linalg2.jl
index 9ddb296788171..3e1410196c7c6 100644
--- a/test/linalg2.jl
+++ b/test/linalg2.jl
@@ -43,6 +43,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
V = convert(Matrix{elty}, V)
C = convert(Matrix{elty}, C)
end
+ ε = eps(abs2(float(one(elty))))
T = Tridiagonal(dl, d, du)
@test size(T, 1) == n
@test size(T) == (n, n)
@@ -95,8 +96,8 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
F = full(W)
@test_approx_eq W*v F*v
iFv = F\v
- @test_approx_eq W\v iFv
- @test_approx_eq det(W) det(F)
+ @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong
+ @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
iWv = similar(iFv)
if elty != Int
Base.LinAlg.solve!(iWv, W, v)
@@ -208,7 +209,7 @@ end
#Triangular matrices
n=12
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
A = convert(Matrix{elty}, randn(n, n))
b = convert(Matrix{elty}, randn(n, 2))
if elty <: Complex
@@ -245,7 +246,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt
end
#Tridiagonal matrices
-for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty})
+for relty in (Float32, Float64), elty in (relty, Complex{relty})
a = convert(Vector{elty}, randn(n-1))
b = convert(Vector{elty}, randn(n))
c = convert(Vector{elty}, randn(n-1))
@@ -263,11 +264,10 @@ for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty})
end
#SymTridiagonal (symmetric tridiagonal) matrices
-for relty in (Float16, Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work
+for relty in (Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work
a = convert(Vector{elty}, randn(n))
b = convert(Vector{elty}, randn(n-1))
if elty <: Complex
- relty==Float16 && continue
a += im*convert(Vector{elty}, randn(n))
b += im*convert(Vector{elty}, randn(n-1))
end
@@ -306,7 +306,7 @@ for elty in (Float32, Float64)
end
#Bidiagonal matrices
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
dv = convert(Vector{elty}, randn(n))
ev = convert(Vector{elty}, randn(n-1))
b = convert(Matrix{elty}, randn(n, 2))
@@ -361,7 +361,7 @@ end
#Diagonal matrices
n=12
-for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty})
+for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
d=convert(Vector{elty}, randn(n))
v=convert(Vector{elty}, randn(n))
U=convert(Matrix{elty}, randn(n,n))
@@ -529,7 +529,7 @@ end
nnorm = 1000
mmat = 100
nmat = 80
-for elty in (Float16, Float32, Float64, BigFloat, Complex{Float16}, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
+for elty in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}, Complex{BigFloat}, Int32, Int64, BigInt)
debug && println(elty)
## Vector