Skip to content

Commit

Permalink
Use pivoting as the default in LU regardless of the element type. (#3…
Browse files Browse the repository at this point in the history
…2989)

For types that weren't subtypes of AbstractFloat, we used to try
to LU factorize without pivoting and only use pivoting when it failed.
This caused large numerical errors when computing the LU for element
types which promoted to float like numbers such as most integers.
The behavior was never documented and is error prone. Hence, this
PR removes the behavior.
  • Loading branch information
andreasnoack committed Aug 23, 2019
1 parent aee3fc2 commit 5af3c2a
Show file tree
Hide file tree
Showing 2 changed files with 14 additions and 28 deletions.
27 changes: 2 additions & 25 deletions stdlib/LinearAlgebra/src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -179,13 +179,6 @@ function generic_lufact!(A::StridedMatrix{T}, ::Val{Pivot} = Val(true);
return LU{T,typeof(A)}(A, ipiv, convert(BlasInt, info))
end

# floating point types doesn't have to be promoted for LU, but should default to pivoting
function lu(A::Union{AbstractMatrix{T}, AbstractMatrix{Complex{T}}},
pivot::Union{Val{false}, Val{true}} = Val(true);
check::Bool = true) where {T<:AbstractFloat}
lu!(copy(A), pivot; check = check)
end

function lutype(T::Type)
# In generic_lufact!, the elements of the lower part of the matrix are
# obtained using the division of two matrix elements. Hence their type can
Expand Down Expand Up @@ -274,26 +267,10 @@ julia> l == F.L && u == F.U && p == F.p
true
```
"""
function lu(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}};
function lu(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}}=Val(true);
check::Bool = true) where T
S = lutype(T)
AA = similar(A, S)
copyto!(AA, A)
lu!(AA, pivot; check = check)
end
# We can't assume an ordered field so we first try without pivoting
function lu(A::AbstractMatrix{T}; check::Bool = true) where T
S = lutype(T)
AA = similar(A, S)
copyto!(AA, A)
F = lu!(AA, Val(false); check = false)
if issuccess(F)
return F
else
AA = similar(A, S)
copyto!(AA, A)
return lu!(AA, Val(true); check = check)
end
lu!(copy_oftype(A, S), pivot; check = check)
end

lu(S::LU) = S
Expand Down
15 changes: 12 additions & 3 deletions stdlib/LinearAlgebra/test/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -301,10 +301,19 @@ include("trickyarithmetic.jl")
@testset "lu with type whose sum is another type" begin
A = TrickyArithmetic.A[1 2; 3 4]
ElT = TrickyArithmetic.D{TrickyArithmetic.C,TrickyArithmetic.C}
B = lu(A)
B = lu(A, Val(false))
@test B isa LinearAlgebra.LU{ElT,Matrix{ElT}}
C = lu(A, Val(false))
@test C isa LinearAlgebra.LU{ElT,Matrix{ElT}}
end

@testset "Issue #30917. Determinant of integer matrix" begin
@test det([1 1 0 0 1 0 0 0
1 0 1 0 0 1 0 0
1 0 0 1 0 0 1 0
0 1 1 1 0 0 0 0
0 1 0 0 0 0 1 1
0 0 1 0 1 0 0 1
0 0 0 1 1 1 0 0
0 0 0 0 1 1 0 1]) 6
end

end # module TestLU

2 comments on commit 5af3c2a

@Keno
Copy link
Member

@Keno Keno commented on 5af3c2a Aug 23, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Broke the doctests? https://build.julialang.org/#/builders/148/builds/1354/steps/3/logs/stdio. (I'm guessing since the doctest error has lu in the output).

@fredrikekre
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please sign in to comment.