Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

define lu!(F::LU,A) for general matrices #1131

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 60 additions & 2 deletions src/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,65 @@ function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupi
end
lu!(A.data, pivot; check, allowsingular)
end

#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
"""
lu!(F::LU, pivot = RowMaximum(); check = true, allowsingular = false) -> LU

`lu!` is the same as [`lu`](@ref), but saves space by overwriting the
input `F`, instead of creating a copy.

!!! compat "Julia 1.12"
Reusing dense `LU` factorizations in `lu!` require Julia 1.12 or later.

# Examples
```jldoctest
julia> A = [4. 3.; 6. 3.]
2×2 Matrix{Float64}:
4.0 3.0
6.0 3.0

julia> F = lu(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
2×2 Matrix{Float64}:
1.0 0.0
0.666667 1.0
U factor:
2×2 Matrix{Float64}:
6.0 3.0
0.0 1.0

julia> B = [8 3; 12 3]
2×2 Matrix{Int64}:
8 3
12 3

julia> F2 = lu!(F,A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
2×2 Matrix{Float64}:
1.0 0.0
0.666667 1.0
U factor:
2×2 Matrix{Float64}:
12.0 3.0
0.0 1.0
```
"""
function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
longemen3000 marked this conversation as resolved.
Show resolved Hide resolved
copyto!(F.factors, A)
return generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
end

#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
copyto!(F.factors, A)
lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
check && _check_lu_success(lpt[3], allowsingular)
return LU{eltype(lpt[1]),typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
end

# for backward compatibility
# TODO: remove towards Julia v2
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, RowMaximum(); check=check)
Expand Down Expand Up @@ -149,7 +208,7 @@ Stacktrace:
"""
lu!(A::AbstractMatrix, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(eltype(A));
check::Bool = true, allowsingular::Bool = false) = generic_lufact!(A, pivot; check, allowsingular)
function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T);
function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T), ipiv::AbstractVector{BlasInt} = Vector{BlasInt}(undef,min(size(A)...));
check::Bool = true, allowsingular::Bool = false) where {T}
check && LAPACK.chkfinite(A)
# Extract values
Expand All @@ -158,7 +217,6 @@ function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,R

# Initialize variables
info = 0
ipiv = Vector{BlasInt}(undef, minmn)
@inbounds begin
for k = 1:minmn
# find index max
Expand Down
1 change: 1 addition & 0 deletions test/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ dimg = randn(n)/2
@test (l*u)[invperm(p),:] ≈ a
@test a * inv(lua) ≈ Matrix(I, n, n)
@test copy(lua) == lua
@test lu!(copy(lua),a) == lua
if eltya <: BlasFloat
# test conversion of LU factorization's numerical type
bft = eltya <: Real ? LinearAlgebra.LU{BigFloat} : LinearAlgebra.LU{Complex{BigFloat}}
Expand Down
Loading