Skip to content

Commit 8cc4216

Browse files
authored
define lu!(F::LU,A) for general matrices (#1131)
alternative to #1125 . Defines `lu!(F::LU{<:Any,<:AbstractMatrix},A)` and `lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}},A)`. closes #1125
2 parents a51f415 + b9d2146 commit 8cc4216

File tree

2 files changed

+61
-2
lines changed

2 files changed

+61
-2
lines changed

src/lu.jl

Lines changed: 60 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -102,6 +102,65 @@ function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupi
102102
end
103103
lu!(A.data, pivot; check, allowsingular)
104104
end
105+
106+
#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
107+
"""
108+
lu!(F::LU, pivot = RowMaximum(); check = true, allowsingular = false) -> LU
109+
110+
`lu!` is the same as [`lu`](@ref), but saves space by overwriting the
111+
input `F`, instead of creating a copy.
112+
113+
!!! compat "Julia 1.12"
114+
Reusing dense `LU` factorizations in `lu!` require Julia 1.12 or later.
115+
116+
# Examples
117+
```jldoctest
118+
julia> A = [4. 3.; 6. 3.]
119+
2×2 Matrix{Float64}:
120+
4.0 3.0
121+
6.0 3.0
122+
123+
julia> F = lu(A)
124+
LU{Float64, Matrix{Float64}, Vector{Int64}}
125+
L factor:
126+
2×2 Matrix{Float64}:
127+
1.0 0.0
128+
0.666667 1.0
129+
U factor:
130+
2×2 Matrix{Float64}:
131+
6.0 3.0
132+
0.0 1.0
133+
134+
julia> B = [8 3; 12 3]
135+
2×2 Matrix{Int64}:
136+
8 3
137+
12 3
138+
139+
julia> F2 = lu!(F,B)
140+
LU{Float64, Matrix{Float64}, Vector{Int64}}
141+
L factor:
142+
2×2 Matrix{Float64}:
143+
1.0 0.0
144+
0.666667 1.0
145+
U factor:
146+
2×2 Matrix{Float64}:
147+
12.0 3.0
148+
0.0 1.0
149+
```
150+
"""
151+
function lu!(F::LU{<:Any,<:AbstractMatrix}, A; check::Bool = true, allowsingular::Bool = false)
152+
copyto!(F.factors, A)
153+
return generic_lufact!(F.factors, lupivottype(eltype(A)), F.ipiv; check, allowsingular)
154+
end
155+
156+
#lu!(F::LU,A) should be dispatched on the type of matrix stored in the LU factorization.
157+
function lu!(F::LU{<:Any,<:StridedMatrix{<:BlasFloat}}, A; check::Bool = true, allowsingular::Bool = false)
158+
copyto!(F.factors, A)
159+
lpt = LAPACK.getrf!(F.factors, F.ipiv; check)
160+
check && _check_lu_success(lpt[3], allowsingular)
161+
return LU{eltype(lpt[1]),typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3])
162+
end
163+
105164
# for backward compatibility
106165
# TODO: remove towards Julia v2
107166
@deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, RowMaximum(); check=check)
@@ -149,7 +208,7 @@ Stacktrace:
149208
"""
150209
lu!(A::AbstractMatrix, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(eltype(A));
151210
check::Bool = true, allowsingular::Bool = false) = generic_lufact!(A, pivot; check, allowsingular)
152-
function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T);
211+
function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T), ipiv::AbstractVector{BlasInt} = Vector{BlasInt}(undef,min(size(A)...));
153212
check::Bool = true, allowsingular::Bool = false) where {T}
154213
check && LAPACK.chkfinite(A)
155214
# Extract values
@@ -158,7 +217,6 @@ function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,R
158217

159218
# Initialize variables
160219
info = 0
161-
ipiv = Vector{BlasInt}(undef, minmn)
162220
@inbounds begin
163221
for k = 1:minmn
164222
# find index max

test/lu.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@ dimg = randn(n)/2
6767
@test (l*u)[invperm(p),:] a
6868
@test a * inv(lua) Matrix(I, n, n)
6969
@test copy(lua) == lua
70+
@test lu!(copy(lua),a) == lua
7071
if eltya <: BlasFloat
7172
# test conversion of LU factorization's numerical type
7273
bft = eltya <: Real ? LinearAlgebra.LU{BigFloat} : LinearAlgebra.LU{Complex{BigFloat}}

0 commit comments

Comments
 (0)