diff --git a/src/lu.jl b/src/lu.jl index 0837ac0..b2d8380 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -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) + 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) @@ -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 @@ -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 diff --git a/test/lu.jl b/test/lu.jl index 56a402d..20358f3 100644 --- a/test/lu.jl +++ b/test/lu.jl @@ -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}}