Skip to content

Commit

Permalink
Fixed broken exception in TriDiag LU and added tests
Browse files Browse the repository at this point in the history
  • Loading branch information
kshyatt committed May 6, 2015
1 parent c4ab817 commit c9f0efd
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 2 deletions.
4 changes: 2 additions & 2 deletions base/linalg/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -289,7 +289,7 @@ end

function At_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
n = size(A,1)
n == size(B,1) || throw(DimentionsMismatch(""))
n == size(B,1) || throw(DimensionMismatch(""))
nrhs = size(B,2)
dl = A.factors.dl
d = A.factors.d
Expand Down Expand Up @@ -322,7 +322,7 @@ end
# Ac_ldiv_B!{T<:Real}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) = At_ldiv_B!(A,B)
function Ac_ldiv_B!{T}(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat)
n = size(A,1)
n == size(B,1) || throw(DimentionsMismatch(""))
n == size(B,1) || throw(DimensionMismatch(""))
nrhs = size(B,2)
dl = A.factors.dl
d = A.factors.d
Expand Down
29 changes: 29 additions & 0 deletions test/linalg/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,23 @@ areal = randn(n,n)/2
aimg = randn(n,n)/2
breal = randn(n,2)/2
bimg = randn(n,2)/2
creal = randn(n)/2
cimg = randn(n)/2
dureal = randn(n-1)/2
duimg = randn(n-1)/2
dlreal = randn(n-1)/2
dlimg = randn(n-1)/2
dreal = randn(n)/2
dimg = randn(n)/2

for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
d = eltya == Int ? Tridiagonal(rand(1:7, n-1), rand(1:7, n), rand(1:7, n-1)) : convert(Tridiagonal{eltya}, eltya <: Complex ? Tridiagonal(complex(dlreal, dlimg), complex(dreal, dimg), complex(dureal, duimg)) : Tridiagonal(dlreal, dreal, dureal))
ε = εa = eps(abs(float(one(eltya))))

for eltyb in (Float32, Float64, Complex64, Complex128, Int)
b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
c = eltyb == Int ? rand(1:5, n) : convert(Vector{eltyb}, eltyb <: Complex ? complex(creal, cimg) : creal)
εb = eps(abs(float(one(eltyb))))
ε = max(εa,εb)

Expand All @@ -36,8 +46,27 @@ debug && println("(Automatic) Square LU decomposition")
@test norm(a*(lua\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
@test norm(a'*(lua'\b) - b, 1) < ε*κ*n*2 # Two because the right hand side has two columns
@test norm(a'*(lua'\a') - a', 1) < ε*κ*n^2
@test norm(a*(lua\c) - c, 1) < ε*κ*n # c is a vector
@test norm(a'*(lua'\c) - c, 1) < ε*κ*n # c is a vector
if eltya <: Real && eltyb <: Real
@test norm(a.'*(lua.'\b) - b,1) < ε*κ*n*2 # Two because the right hand side has two columns
@test norm(a.'*(lua.'\c) - c,1) < ε*κ*n
end
@test_approx_eq full(lua) a

debug && println("Tridiagonal LU")
κd = cond(full(d),1)
lud = lufact(d)
@test norm(d*(lud\b) - b, 1) < ε*κd*n*2 # Two because the right hand side has two columns
if eltya <: Real
@test norm((lud.'\b) - full(d.')\b, 1) < ε*κd*n*2 # Two because the right hand side has two columns
end
if eltya <: Complex
@test norm((lud'\b) - full(d')\b, 1) < ε*κd*n*2 # Two because the right hand side has two columns
end
f = zeros(eltyb, n+1)
@test_throws DimensionMismatch lud\f

debug && println("Thin LU")
lua = @inferred lufact(a[:,1:n1])
@test_approx_eq lua[:L]*lua[:U] lua[:P]*a[:,1:n1]
Expand Down

0 comments on commit c9f0efd

Please sign in to comment.