Line | Exclusive | Inclusive | Code |
---|---|---|---|
1 | # This file is a part of Julia. License is MIT: https://julialang.org/license | ||
2 | |||
3 | #################### | ||
4 | # LU Factorization # | ||
5 | #################### | ||
6 | """ | ||
7 | LU <: Factorization | ||
8 | |||
9 | Matrix factorization type of the `LU` factorization of a square matrix `A`. This | ||
10 | is the return type of [`lu`](@ref), the corresponding matrix factorization function. | ||
11 | |||
12 | The individual components of the factorization `F::LU` can be accessed via [`getproperty`](@ref): | ||
13 | |||
14 | | Component | Description | | ||
15 | |:----------|:-----------------------------------------| | ||
16 | | `F.L` | `L` (unit lower triangular) part of `LU` | | ||
17 | | `F.U` | `U` (upper triangular) part of `LU` | | ||
18 | | `F.p` | (right) permutation `Vector` | | ||
19 | | `F.P` | (right) permutation `Matrix` | | ||
20 | |||
21 | Iterating the factorization produces the components `F.L`, `F.U`, and `F.p`. | ||
22 | |||
23 | # Examples | ||
24 | ```jldoctest | ||
25 | julia> A = [4 3; 6 3] | ||
26 | 2×2 Matrix{Int64}: | ||
27 | 4 3 | ||
28 | 6 3 | ||
29 | |||
30 | julia> F = lu(A) | ||
31 | LU{Float64, Matrix{Float64}, Vector{Int64}} | ||
32 | L factor: | ||
33 | 2×2 Matrix{Float64}: | ||
34 | 1.0 0.0 | ||
35 | 0.666667 1.0 | ||
36 | U factor: | ||
37 | 2×2 Matrix{Float64}: | ||
38 | 6.0 3.0 | ||
39 | 0.0 1.0 | ||
40 | |||
41 | julia> F.L * F.U == A[F.p, :] | ||
42 | true | ||
43 | |||
44 | julia> l, u, p = lu(A); # destructuring via iteration | ||
45 | |||
46 | julia> l == F.L && u == F.U && p == F.p | ||
47 | true | ||
48 | ``` | ||
49 | """ | ||
50 | struct LU{T,S<:AbstractMatrix{T},P<:AbstractVector{<:Integer}} <: Factorization{T} | ||
51 | factors::S | ||
52 | ipiv::P | ||
53 | info::BlasInt | ||
54 | |||
55 | function LU{T,S,P}(factors, ipiv, info) where {T, S<:AbstractMatrix{T}, P<:AbstractVector{<:Integer}} | ||
56 | require_one_based_indexing(factors) | ||
57 | new{T,S,P}(factors, ipiv, info) | ||
58 | end | ||
59 | end | ||
60 | LU(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, info::BlasInt) where {T} = | ||
61 | LU{T,typeof(factors),typeof(ipiv)}(factors, ipiv, info) | ||
62 | LU{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T} = | ||
63 | LU(convert(AbstractMatrix{T}, factors), ipiv, BlasInt(info)) | ||
64 | # backwards-compatible constructors (remove with Julia 2.0) | ||
65 | @deprecate(LU{T,S}(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, | ||
66 | info::BlasInt) where {T,S}, | ||
67 | LU{T,S,typeof(ipiv)}(factors, ipiv, info), false) | ||
68 | |||
69 | # iteration for destructuring into components | ||
70 | Base.iterate(S::LU) = (S.L, Val(:U)) | ||
71 | Base.iterate(S::LU, ::Val{:U}) = (S.U, Val(:p)) | ||
72 | Base.iterate(S::LU, ::Val{:p}) = (S.p, Val(:done)) | ||
73 | Base.iterate(S::LU, ::Val{:done}) = nothing | ||
74 | |||
75 | # LU prefers transpose over adjoint in the real case, override the generic fallback | ||
76 | adjoint(F::LU{<:Real}) = TransposeFactorization(F) | ||
77 | transpose(F::LU{<:Real}) = TransposeFactorization(F) | ||
78 | |||
79 | # the following method is meant to catch calls to lu!(A::LAPACKArray) without a pivoting stategy | ||
80 | lu!(A::StridedMatrix{<:BlasFloat}; check::Bool = true) = lu!(A, RowMaximum(); check=check) | ||
81 | function lu!(A::StridedMatrix{T}, ::RowMaximum; check::Bool = true) where {T<:BlasFloat} | ||
82 | lpt = LAPACK.getrf!(A; check) | ||
83 | check && checknonsingular(lpt[3]) | ||
84 | return LU{T,typeof(lpt[1]),typeof(lpt[2])}(lpt[1], lpt[2], lpt[3]) | ||
85 | end | ||
86 | function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T); check::Bool = true) where {T} | ||
87 | copytri!(A.data, A.uplo, isa(A, Hermitian)) | ||
88 | lu!(A.data, pivot; check = check) | ||
89 | end | ||
90 | # for backward compatibility | ||
91 | # TODO: remove towards Julia v2 | ||
92 | @deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{true}; check::Bool = true) lu!(A, RowMaximum(); check=check) | ||
93 | @deprecate lu!(A::Union{StridedMatrix,HermOrSym,Tridiagonal}, ::Val{false}; check::Bool = true) lu!(A, NoPivot(); check=check) | ||
94 | |||
95 | """ | ||
96 | lu!(A, pivot = RowMaximum(); check = true) -> LU | ||
97 | |||
98 | `lu!` is the same as [`lu`](@ref), but saves space by overwriting the | ||
99 | input `A`, instead of creating a copy. An [`InexactError`](@ref) | ||
100 | exception is thrown if the factorization produces a number not representable by the | ||
101 | element type of `A`, e.g. for integer types. | ||
102 | |||
103 | # Examples | ||
104 | ```jldoctest | ||
105 | julia> A = [4. 3.; 6. 3.] | ||
106 | 2×2 Matrix{Float64}: | ||
107 | 4.0 3.0 | ||
108 | 6.0 3.0 | ||
109 | |||
110 | julia> F = lu!(A) | ||
111 | LU{Float64, Matrix{Float64}, Vector{Int64}} | ||
112 | L factor: | ||
113 | 2×2 Matrix{Float64}: | ||
114 | 1.0 0.0 | ||
115 | 0.666667 1.0 | ||
116 | U factor: | ||
117 | 2×2 Matrix{Float64}: | ||
118 | 6.0 3.0 | ||
119 | 0.0 1.0 | ||
120 | |||
121 | julia> iA = [4 3; 6 3] | ||
122 | 2×2 Matrix{Int64}: | ||
123 | 4 3 | ||
124 | 6 3 | ||
125 | |||
126 | julia> lu!(iA) | ||
127 | ERROR: InexactError: Int64(0.6666666666666666) | ||
128 | Stacktrace: | ||
129 | [...] | ||
130 | ``` | ||
131 | """ | ||
132 | lu!(A::AbstractMatrix, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(eltype(A)); check::Bool = true) = | ||
133 | generic_lufact!(A, pivot; check = check) | ||
134 | function generic_lufact!(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T); | ||
135 | check::Bool = true) where {T} | ||
136 | check && LAPACK.chkfinite(A) | ||
137 | # Extract values | ||
138 | m, n = size(A) | ||
139 | minmn = min(m,n) | ||
140 | |||
141 | # Initialize variables | ||
142 | info = 0 | ||
143 | ipiv = Vector{BlasInt}(undef, minmn) | ||
144 | @inbounds begin | ||
145 | for k = 1:minmn | ||
146 | # find index max | ||
147 | kp = k | ||
148 | if pivot === RowMaximum() && k < m | ||
149 | amax = abs(A[k, k]) | ||
150 | for i = k+1:m | ||
151 | absi = abs(A[i,k]) | ||
152 | if absi > amax | ||
153 | kp = i | ||
154 | amax = absi | ||
155 | end | ||
156 | end | ||
157 | elseif pivot === RowNonZero() | ||
158 | for i = k:m | ||
159 | if !iszero(A[i,k]) | ||
160 | kp = i | ||
161 | break | ||
162 | end | ||
163 | end | ||
164 | end | ||
165 | ipiv[k] = kp | ||
166 | if !iszero(A[kp,k]) | ||
167 | if k != kp | ||
168 | # Interchange | ||
169 | for i = 1:n | ||
170 | tmp = A[k,i] | ||
171 | A[k,i] = A[kp,i] | ||
172 | A[kp,i] = tmp | ||
173 | end | ||
174 | end | ||
175 | # Scale first column | ||
176 | Akkinv = inv(A[k,k]) | ||
177 | for i = k+1:m | ||
178 | A[i,k] *= Akkinv | ||
179 | end | ||
180 | elseif info == 0 | ||
181 | info = k | ||
182 | end | ||
183 | # Update the rest | ||
184 | for j = k+1:n | ||
185 | for i = k+1:m | ||
186 | A[i,j] -= A[i,k]*A[k,j] | ||
187 | end | ||
188 | end | ||
189 | end | ||
190 | end | ||
191 | check && checknonsingular(info, pivot) | ||
192 | return LU{T,typeof(A),typeof(ipiv)}(A, ipiv, convert(BlasInt, info)) | ||
193 | end | ||
194 | |||
195 | function lutype(T::Type) | ||
196 | # In generic_lufact!, the elements of the lower part of the matrix are | ||
197 | # obtained using the division of two matrix elements. Hence their type can | ||
198 | # be different (e.g. the division of two types with the same unit is a type | ||
199 | # without unit). | ||
200 | # The elements of the upper part are obtained by U - U * L | ||
201 | # where U is an upper part element and L is a lower part element. | ||
202 | # Therefore, the types LT, UT should be invariant under the map: | ||
203 | # (LT, UT) -> begin | ||
204 | # L = oneunit(UT) / oneunit(UT) | ||
205 | # U = oneunit(UT) - oneunit(UT) * L | ||
206 | # typeof(L), typeof(U) | ||
207 | # end | ||
208 | # The following should handle most cases | ||
209 | UT = typeof(oneunit(T) - oneunit(T) * (oneunit(T) / (oneunit(T) + zero(T)))) | ||
210 | LT = typeof(oneunit(UT) / oneunit(UT)) | ||
211 | S = promote_type(T, LT, UT) | ||
212 | end | ||
213 | |||
214 | lupivottype(::Type{T}) where {T} = RowMaximum() | ||
215 | |||
216 | # for all other types we must promote to a type which is stable under division | ||
217 | """ | ||
218 | lu(A, pivot = RowMaximum(); check = true) -> F::LU | ||
219 | |||
220 | Compute the LU factorization of `A`. | ||
221 | |||
222 | When `check = true`, an error is thrown if the decomposition fails. | ||
223 | When `check = false`, responsibility for checking the decomposition's | ||
224 | validity (via [`issuccess`](@ref)) lies with the user. | ||
225 | |||
226 | In most cases, if `A` is a subtype `S` of `AbstractMatrix{T}` with an element | ||
227 | type `T` supporting `+`, `-`, `*` and `/`, the return type is `LU{T,S{T}}`. | ||
228 | |||
229 | In general, LU factorization involves a permutation of the rows of the matrix | ||
230 | (corresponding to the `F.p` output described below), known as "pivoting" (because it | ||
231 | corresponds to choosing which row contains the "pivot", the diagonal entry of `F.U`). | ||
232 | One of the following pivoting strategies can be selected via the optional `pivot` argument: | ||
233 | |||
234 | * `RowMaximum()` (default): the standard pivoting strategy; the pivot corresponds | ||
235 | to the element of maximum absolute value among the remaining, to be factorized rows. | ||
236 | This pivoting strategy requires the element type to also support [`abs`](@ref) and | ||
237 | [`<`](@ref). (This is generally the only numerically stable option for floating-point | ||
238 | matrices.) | ||
239 | * `RowNonZero()`: the pivot corresponds to the first non-zero element among the remaining, | ||
240 | to be factorized rows. (This corresponds to the typical choice in hand calculations, and | ||
241 | is also useful for more general algebraic number types that support [`iszero`](@ref) but | ||
242 | not `abs` or `<`.) | ||
243 | * `NoPivot()`: pivoting turned off (may fail if a zero entry is encountered). | ||
244 | |||
245 | The individual components of the factorization `F` can be accessed via [`getproperty`](@ref): | ||
246 | |||
247 | | Component | Description | | ||
248 | |:----------|:------------------------------------| | ||
249 | | `F.L` | `L` (lower triangular) part of `LU` | | ||
250 | | `F.U` | `U` (upper triangular) part of `LU` | | ||
251 | | `F.p` | (right) permutation `Vector` | | ||
252 | | `F.P` | (right) permutation `Matrix` | | ||
253 | |||
254 | Iterating the factorization produces the components `F.L`, `F.U`, and `F.p`. | ||
255 | |||
256 | The relationship between `F` and `A` is | ||
257 | |||
258 | `F.L*F.U == A[F.p, :]` | ||
259 | |||
260 | `F` further supports the following functions: | ||
261 | |||
262 | | Supported function | `LU` | `LU{T,Tridiagonal{T}}` | | ||
263 | |:---------------------------------|:-----|:-----------------------| | ||
264 | | [`/`](@ref) | ✓ | | | ||
265 | | [`\\`](@ref) | ✓ | ✓ | | ||
266 | | [`inv`](@ref) | ✓ | ✓ | | ||
267 | | [`det`](@ref) | ✓ | ✓ | | ||
268 | | [`logdet`](@ref) | ✓ | ✓ | | ||
269 | | [`logabsdet`](@ref) | ✓ | ✓ | | ||
270 | | [`size`](@ref) | ✓ | ✓ | | ||
271 | |||
272 | # Examples | ||
273 | ```jldoctest | ||
274 | julia> A = [4 3; 6 3] | ||
275 | 2×2 Matrix{Int64}: | ||
276 | 4 3 | ||
277 | 6 3 | ||
278 | |||
279 | julia> F = lu(A) | ||
280 | LU{Float64, Matrix{Float64}, Vector{Int64}} | ||
281 | L factor: | ||
282 | 2×2 Matrix{Float64}: | ||
283 | 1.0 0.0 | ||
284 | 0.666667 1.0 | ||
285 | U factor: | ||
286 | 2×2 Matrix{Float64}: | ||
287 | 6.0 3.0 | ||
288 | 0.0 1.0 | ||
289 | |||
290 | julia> F.L * F.U == A[F.p, :] | ||
291 | true | ||
292 | |||
293 | julia> l, u, p = lu(A); # destructuring via iteration | ||
294 | |||
295 | julia> l == F.L && u == F.U && p == F.p | ||
296 | true | ||
297 | ``` | ||
298 | """ | ||
299 | function lu(A::AbstractMatrix{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T); check::Bool = true) where {T} | ||
300 | lu!(_lucopy(A, lutype(T)), pivot; check = check) | ||
301 | end | ||
302 | # TODO: remove for Julia v2.0 | ||
303 | @deprecate lu(A::AbstractMatrix, ::Val{true}; check::Bool = true) lu(A, RowMaximum(); check=check) | ||
304 | @deprecate lu(A::AbstractMatrix, ::Val{false}; check::Bool = true) lu(A, NoPivot(); check=check) | ||
305 | |||
306 | _lucopy(A::AbstractMatrix, T) = copy_similar(A, T) | ||
307 | _lucopy(A::HermOrSym, T) = copymutable_oftype(A, T) | ||
308 | _lucopy(A::Tridiagonal, T) = copymutable_oftype(A, T) | ||
309 | |||
310 | lu(S::LU) = S | ||
311 | function lu(x::Number; check::Bool=true) | ||
312 | info = x == 0 ? one(BlasInt) : zero(BlasInt) | ||
313 | check && checknonsingular(info) | ||
314 | return LU(fill(x, 1, 1), BlasInt[1], info) | ||
315 | end | ||
316 | |||
317 | function LU{T}(F::LU) where T | ||
318 | M = convert(AbstractMatrix{T}, F.factors) | ||
319 | LU{T,typeof(M),typeof(F.ipiv)}(M, F.ipiv, F.info) | ||
320 | end | ||
321 | LU{T,S,P}(F::LU) where {T,S,P} = LU{T,S,P}(convert(S, F.factors), convert(P, F.ipiv), F.info) | ||
322 | Factorization{T}(F::LU{T}) where {T} = F | ||
323 | Factorization{T}(F::LU) where {T} = LU{T}(F) | ||
324 | |||
325 | copy(A::LU{T,S,P}) where {T,S,P} = LU{T,S,P}(copy(A.factors), copy(A.ipiv), A.info) | ||
326 | |||
327 | size(A::LU) = size(getfield(A, :factors)) | ||
328 | size(A::LU, i::Integer) = size(getfield(A, :factors), i) | ||
329 | |||
330 | function ipiv2perm(v::AbstractVector{T}, maxi::Integer) where T | ||
331 | require_one_based_indexing(v) | ||
332 | p = T[1:maxi;] | ||
333 | @inbounds for i in 1:length(v) | ||
334 | p[i], p[v[i]] = p[v[i]], p[i] | ||
335 | end | ||
336 | return p | ||
337 | end | ||
338 | |||
339 | function getproperty(F::LU{T}, d::Symbol) where T | ||
340 | m, n = size(F) | ||
341 | if d === :L | ||
342 | L = tril!(getfield(F, :factors)[1:m, 1:min(m,n)]) | ||
343 | for i = 1:min(m,n); L[i,i] = one(T); end | ||
344 | return L | ||
345 | elseif d === :U | ||
346 | return triu!(getfield(F, :factors)[1:min(m,n), 1:n]) | ||
347 | elseif d === :p | ||
348 | return ipiv2perm(getfield(F, :ipiv), m) | ||
349 | elseif d === :P | ||
350 | return Matrix{T}(I, m, m)[:,invperm(F.p)] | ||
351 | else | ||
352 | getfield(F, d) | ||
353 | end | ||
354 | end | ||
355 | |||
356 | Base.propertynames(F::LU, private::Bool=false) = | ||
357 | (:L, :U, :p, :P, (private ? fieldnames(typeof(F)) : ())...) | ||
358 | |||
359 | issuccess(F::LU) = F.info == 0 | ||
360 | |||
361 | function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LU) | ||
362 | if issuccess(F) | ||
363 | summary(io, F); println(io) | ||
364 | println(io, "L factor:") | ||
365 | show(io, mime, F.L) | ||
366 | println(io, "\nU factor:") | ||
367 | show(io, mime, F.U) | ||
368 | else | ||
369 | print(io, "Failed factorization of type $(typeof(F))") | ||
370 | end | ||
371 | end | ||
372 | |||
373 | _apply_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, 1 : length(A.ipiv), B) | ||
374 | _apply_inverse_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, length(A.ipiv) : -1 : 1, B) | ||
375 | |||
376 | function _ipiv_rows!(A::LU, order::OrdinalRange, B::AbstractVecOrMat) | ||
377 | for i = order | ||
378 | if i != A.ipiv[i] | ||
379 | _swap_rows!(B, i, A.ipiv[i]) | ||
380 | end | ||
381 | end | ||
382 | B | ||
383 | end | ||
384 | |||
385 | function _swap_rows!(B::AbstractVector, i::Integer, j::Integer) | ||
386 | B[i], B[j] = B[j], B[i] | ||
387 | B | ||
388 | end | ||
389 | |||
390 | function _swap_rows!(B::AbstractMatrix, i::Integer, j::Integer) | ||
391 | for col = 1 : size(B, 2) | ||
392 | B[i,col], B[j,col] = B[j,col], B[i,col] | ||
393 | end | ||
394 | B | ||
395 | end | ||
396 | |||
397 | _apply_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, 1 : length(A.ipiv), B) | ||
398 | _apply_inverse_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, length(A.ipiv) : -1 : 1, B) | ||
399 | |||
400 | function _ipiv_cols!(A::LU, order::OrdinalRange, B::AbstractVecOrMat) | ||
401 | for i = order | ||
402 | if i != A.ipiv[i] | ||
403 | _swap_cols!(B, i, A.ipiv[i]) | ||
404 | end | ||
405 | end | ||
406 | B | ||
407 | end | ||
408 | |||
409 | function _swap_cols!(B::AbstractVector, i::Integer, j::Integer) | ||
410 | _swap_rows!(B, i, j) | ||
411 | end | ||
412 | |||
413 | function _swap_cols!(B::AbstractMatrix, i::Integer, j::Integer) | ||
414 | for row = 1 : size(B, 1) | ||
415 | B[row,i], B[row,j] = B[row,j], B[row,i] | ||
416 | end | ||
417 | B | ||
418 | end | ||
419 | |||
420 | function rdiv!(A::AbstractVecOrMat, B::LU) | ||
421 | rdiv!(rdiv!(A, UpperTriangular(B.factors)), UnitLowerTriangular(B.factors)) | ||
422 | _apply_inverse_ipiv_cols!(B, A) | ||
423 | end | ||
424 | |||
425 | 9 (3 %) |
9 (100 %)
samples spent calling
getrs!
ldiv!(A::LU{T,<:StridedMatrix}, B::StridedVecOrMat{T}) where {T<:BlasFloat} =
|
|
426 | LAPACK.getrs!('N', A.factors, A.ipiv, B) | ||
427 | |||
428 | function ldiv!(A::LU, B::AbstractVecOrMat) | ||
429 | _apply_ipiv_rows!(A, B) | ||
430 | ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B)) | ||
431 | end | ||
432 | |||
433 | ldiv!(transA::TransposeFactorization{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} = | ||
434 | (A = transA.parent; LAPACK.getrs!('T', A.factors, A.ipiv, B)) | ||
435 | |||
436 | function ldiv!(transA::TransposeFactorization{<:Any,<:LU}, B::AbstractVecOrMat) | ||
437 | A = transA.parent | ||
438 | ldiv!(transpose(UnitLowerTriangular(A.factors)), ldiv!(transpose(UpperTriangular(A.factors)), B)) | ||
439 | _apply_inverse_ipiv_rows!(A, B) | ||
440 | end | ||
441 | |||
442 | ldiv!(adjA::AdjointFactorization{T,<:LU{T,<:StridedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} = | ||
443 | (A = adjA.parent; LAPACK.getrs!('C', A.factors, A.ipiv, B)) | ||
444 | |||
445 | function ldiv!(adjA::AdjointFactorization{<:Any,<:LU}, B::AbstractVecOrMat) | ||
446 | A = adjA.parent | ||
447 | ldiv!(adjoint(UnitLowerTriangular(A.factors)), ldiv!(adjoint(UpperTriangular(A.factors)), B)) | ||
448 | _apply_inverse_ipiv_rows!(A, B) | ||
449 | end | ||
450 | |||
451 | (\)(A::AdjointFactorization{T,<:LU{T,<:StridedMatrix}}, B::Adjoint{T,<:StridedVecOrMat{T}}) where {T<:BlasComplex} = | ||
452 | LAPACK.getrs!('C', A.parent.factors, A.parent.ipiv, copy(B)) | ||
453 | (\)(A::TransposeFactorization{T,<:LU{T,<:StridedMatrix}}, B::Transpose{T,<:StridedVecOrMat{T}}) where {T<:BlasFloat} = | ||
454 | LAPACK.getrs!('T', A.parent.factors, A.parent.ipiv, copy(B)) | ||
455 | |||
456 | function det(F::LU{T}) where T | ||
457 | n = checksquare(F) | ||
458 | issuccess(F) || return zero(T) | ||
459 | P = one(T) | ||
460 | c = 0 | ||
461 | @inbounds for i = 1:n | ||
462 | P *= F.factors[i,i] | ||
463 | if F.ipiv[i] != i | ||
464 | c += 1 | ||
465 | end | ||
466 | end | ||
467 | s = (isodd(c) ? -one(T) : one(T)) | ||
468 | return P * s | ||
469 | end | ||
470 | |||
471 | function logabsdet(F::LU{T}) where T # return log(abs(det)) and sign(det) | ||
472 | n = checksquare(F) | ||
473 | issuccess(F) || return log(zero(real(T))), log(one(T)) | ||
474 | c = 0 | ||
475 | P = one(T) | ||
476 | abs_det = zero(real(T)) | ||
477 | @inbounds for i = 1:n | ||
478 | dg_ii = F.factors[i,i] | ||
479 | P *= sign(dg_ii) | ||
480 | if F.ipiv[i] != i | ||
481 | c += 1 | ||
482 | end | ||
483 | abs_det += log(abs(dg_ii)) | ||
484 | end | ||
485 | s = ifelse(isodd(c), -one(real(T)), one(real(T))) * P | ||
486 | abs_det, s | ||
487 | end | ||
488 | |||
489 | inv!(A::LU{<:BlasFloat,<:StridedMatrix}) = | ||
490 | LAPACK.getri!(A.factors, A.ipiv) | ||
491 | inv!(A::LU{T,<:StridedMatrix}) where {T} = | ||
492 | ldiv!(A.factors, copy(A), Matrix{T}(I, size(A, 1), size(A, 1))) | ||
493 | inv(A::LU{<:BlasFloat,<:StridedMatrix}) = inv!(copy(A)) | ||
494 | |||
495 | # Tridiagonal | ||
496 | |||
497 | # See dgttrf.f | ||
498 | function lu!(A::Tridiagonal{T,V}, pivot::Union{RowMaximum,NoPivot} = RowMaximum(); check::Bool = true) where {T,V} | ||
499 | # Extract values | ||
500 | n = size(A, 1) | ||
501 | |||
502 | # Initialize variables | ||
503 | info = 0 | ||
504 | ipiv = Vector{BlasInt}(undef, n) | ||
505 | dl = A.dl | ||
506 | d = A.d | ||
507 | du = A.du | ||
508 | if dl === du | ||
509 | throw(ArgumentError("off-diagonals of `A` must not alias")) | ||
510 | end | ||
511 | # Check if Tridiagonal matrix already has du2 for pivoting | ||
512 | has_du2_defined = isdefined(A, :du2) && length(A.du2) == max(0, n-2) | ||
513 | if has_du2_defined | ||
514 | du2 = A.du2::V | ||
515 | else | ||
516 | du2 = similar(d, max(0, n-2))::V | ||
517 | end | ||
518 | fill!(du2, 0) | ||
519 | |||
520 | @inbounds begin | ||
521 | for i = 1:n | ||
522 | ipiv[i] = i | ||
523 | end | ||
524 | for i = 1:n-2 | ||
525 | # pivot or not? | ||
526 | if pivot === NoPivot() || abs(d[i]) >= abs(dl[i]) | ||
527 | # No interchange | ||
528 | if d[i] != 0 | ||
529 | fact = dl[i]/d[i] | ||
530 | dl[i] = fact | ||
531 | d[i+1] -= fact*du[i] | ||
532 | du2[i] = 0 | ||
533 | end | ||
534 | else | ||
535 | # Interchange | ||
536 | fact = d[i]/dl[i] | ||
537 | d[i] = dl[i] | ||
538 | dl[i] = fact | ||
539 | tmp = du[i] | ||
540 | du[i] = d[i+1] | ||
541 | d[i+1] = tmp - fact*d[i+1] | ||
542 | du2[i] = du[i+1] | ||
543 | du[i+1] = -fact*du[i+1] | ||
544 | ipiv[i] = i+1 | ||
545 | end | ||
546 | end | ||
547 | if n > 1 | ||
548 | i = n-1 | ||
549 | if pivot === NoPivot() || abs(d[i]) >= abs(dl[i]) | ||
550 | if d[i] != 0 | ||
551 | fact = dl[i]/d[i] | ||
552 | dl[i] = fact | ||
553 | d[i+1] -= fact*du[i] | ||
554 | end | ||
555 | else | ||
556 | fact = d[i]/dl[i] | ||
557 | d[i] = dl[i] | ||
558 | dl[i] = fact | ||
559 | tmp = du[i] | ||
560 | du[i] = d[i+1] | ||
561 | d[i+1] = tmp - fact*d[i+1] | ||
562 | ipiv[i] = i+1 | ||
563 | end | ||
564 | end | ||
565 | # check for a zero on the diagonal of U | ||
566 | for i = 1:n | ||
567 | if d[i] == 0 | ||
568 | info = i | ||
569 | break | ||
570 | end | ||
571 | end | ||
572 | end | ||
573 | B = has_du2_defined ? A : Tridiagonal{T,V}(dl, d, du, du2) | ||
574 | check && checknonsingular(info, pivot) | ||
575 | return LU{T,Tridiagonal{T,V},typeof(ipiv)}(B, ipiv, convert(BlasInt, info)) | ||
576 | end | ||
577 | |||
578 | factorize(A::Tridiagonal) = lu(A) | ||
579 | |||
580 | function getproperty(F::LU{T,Tridiagonal{T,V}}, d::Symbol) where {T,V} | ||
581 | m, n = size(F) | ||
582 | if d === :L | ||
583 | dl = getfield(getfield(F, :factors), :dl) | ||
584 | L = Array(Bidiagonal(fill!(similar(dl, n), one(T)), dl, d)) | ||
585 | for i = 2:n | ||
586 | tmp = L[getfield(F, :ipiv)[i], 1:i - 1] | ||
587 | L[getfield(F, :ipiv)[i], 1:i - 1] = L[i, 1:i - 1] | ||
588 | L[i, 1:i - 1] = tmp | ||
589 | end | ||
590 | return L | ||
591 | elseif d === :U | ||
592 | U = Array(Bidiagonal(getfield(getfield(F, :factors), :d), getfield(getfield(F, :factors), :du), d)) | ||
593 | for i = 1:n - 2 | ||
594 | U[i,i + 2] = getfield(getfield(F, :factors), :du2)[i] | ||
595 | end | ||
596 | return U | ||
597 | elseif d === :p | ||
598 | return ipiv2perm(getfield(F, :ipiv), m) | ||
599 | elseif d === :P | ||
600 | return Matrix{T}(I, m, m)[:,invperm(F.p)] | ||
601 | end | ||
602 | return getfield(F, d) | ||
603 | end | ||
604 | |||
605 | # See dgtts2.f | ||
606 | function ldiv!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V} | ||
607 | require_one_based_indexing(B) | ||
608 | n = size(A,1) | ||
609 | if n != size(B,1) | ||
610 | throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows")) | ||
611 | end | ||
612 | nrhs = size(B,2) | ||
613 | dl = A.factors.dl | ||
614 | d = A.factors.d | ||
615 | du = A.factors.du | ||
616 | du2 = A.factors.du2 | ||
617 | ipiv = A.ipiv | ||
618 | @inbounds begin | ||
619 | for j = 1:nrhs | ||
620 | for i = 1:n-1 | ||
621 | ip = ipiv[i] | ||
622 | tmp = B[i+1-ip+i,j] - dl[i]*B[ip,j] | ||
623 | B[i,j] = B[ip,j] | ||
624 | B[i+1,j] = tmp | ||
625 | end | ||
626 | B[n,j] /= d[n] | ||
627 | if n > 1 | ||
628 | B[n-1,j] = (B[n-1,j] - du[n-1]*B[n,j])/d[n-1] | ||
629 | end | ||
630 | for i = n-2:-1:1 | ||
631 | B[i,j] = (B[i,j] - du[i]*B[i+1,j] - du2[i]*B[i+2,j])/d[i] | ||
632 | end | ||
633 | end | ||
634 | end | ||
635 | return B | ||
636 | end | ||
637 | |||
638 | function ldiv!(transA::TransposeFactorization{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} | ||
639 | require_one_based_indexing(B) | ||
640 | A = transA.parent | ||
641 | n = size(A,1) | ||
642 | if n != size(B,1) | ||
643 | throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows")) | ||
644 | end | ||
645 | nrhs = size(B,2) | ||
646 | dl = A.factors.dl | ||
647 | d = A.factors.d | ||
648 | du = A.factors.du | ||
649 | du2 = A.factors.du2 | ||
650 | ipiv = A.ipiv | ||
651 | @inbounds begin | ||
652 | for j = 1:nrhs | ||
653 | B[1,j] /= d[1] | ||
654 | if n > 1 | ||
655 | B[2,j] = (B[2,j] - du[1]*B[1,j])/d[2] | ||
656 | end | ||
657 | for i = 3:n | ||
658 | B[i,j] = (B[i,j] - du[i-1]*B[i-1,j] - du2[i-2]*B[i-2,j])/d[i] | ||
659 | end | ||
660 | for i = n-1:-1:1 | ||
661 | if ipiv[i] == i | ||
662 | B[i,j] = B[i,j] - dl[i]*B[i+1,j] | ||
663 | else | ||
664 | tmp = B[i+1,j] | ||
665 | B[i+1,j] = B[i,j] - dl[i]*tmp | ||
666 | B[i,j] = tmp | ||
667 | end | ||
668 | end | ||
669 | end | ||
670 | end | ||
671 | return B | ||
672 | end | ||
673 | |||
674 | # Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B) | ||
675 | function ldiv!(adjA::AdjointFactorization{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} | ||
676 | require_one_based_indexing(B) | ||
677 | A = adjA.parent | ||
678 | n = size(A,1) | ||
679 | if n != size(B,1) | ||
680 | throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows")) | ||
681 | end | ||
682 | nrhs = size(B,2) | ||
683 | dl = A.factors.dl | ||
684 | d = A.factors.d | ||
685 | du = A.factors.du | ||
686 | du2 = A.factors.du2 | ||
687 | ipiv = A.ipiv | ||
688 | @inbounds begin | ||
689 | for j = 1:nrhs | ||
690 | B[1,j] /= conj(d[1]) | ||
691 | if n > 1 | ||
692 | B[2,j] = (B[2,j] - conj(du[1])*B[1,j])/conj(d[2]) | ||
693 | end | ||
694 | for i = 3:n | ||
695 | B[i,j] = (B[i,j] - conj(du[i-1])*B[i-1,j] - conj(du2[i-2])*B[i-2,j])/conj(d[i]) | ||
696 | end | ||
697 | for i = n-1:-1:1 | ||
698 | if ipiv[i] == i | ||
699 | B[i,j] = B[i,j] - conj(dl[i])*B[i+1,j] | ||
700 | else | ||
701 | tmp = B[i+1,j] | ||
702 | B[i+1,j] = B[i,j] - conj(dl[i])*tmp | ||
703 | B[i,j] = tmp | ||
704 | end | ||
705 | end | ||
706 | end | ||
707 | end | ||
708 | return B | ||
709 | end | ||
710 | |||
711 | rdiv!(B::AbstractMatrix, A::LU) = transpose(ldiv!(transpose(A), transpose(B))) | ||
712 | |||
713 | # Conversions | ||
714 | AbstractMatrix(F::LU) = (F.L * F.U)[invperm(F.p),:] | ||
715 | AbstractArray(F::LU) = AbstractMatrix(F) | ||
716 | Matrix(F::LU) = Array(AbstractArray(F)) | ||
717 | Array(F::LU) = Matrix(F) | ||
718 | |||
719 | function Tridiagonal(F::LU{T,Tridiagonal{T,V}}) where {T,V} | ||
720 | n = size(F, 1) | ||
721 | |||
722 | dl = copy(F.factors.dl) | ||
723 | d = copy(F.factors.d) | ||
724 | du = copy(F.factors.du) | ||
725 | du2 = copy(F.factors.du2) | ||
726 | |||
727 | for i = n - 1:-1:1 | ||
728 | li = dl[i] | ||
729 | dl[i] = li*d[i] | ||
730 | d[i + 1] += li*du[i] | ||
731 | if i < n - 1 | ||
732 | du[i + 1] += li*du2[i] | ||
733 | end | ||
734 | |||
735 | if F.ipiv[i] != i | ||
736 | tmp = dl[i] | ||
737 | dl[i] = d[i] | ||
738 | d[i] = tmp | ||
739 | |||
740 | tmp = d[i + 1] | ||
741 | d[i + 1] = du[i] | ||
742 | du[i] = tmp | ||
743 | |||
744 | if i < n - 1 | ||
745 | tmp = du[i + 1] | ||
746 | du[i + 1] = du2[i] | ||
747 | du2[i] = tmp | ||
748 | end | ||
749 | end | ||
750 | end | ||
751 | return Tridiagonal(dl, d, du) | ||
752 | end | ||
753 | AbstractMatrix(F::LU{T,Tridiagonal{T,V}}) where {T,V} = Tridiagonal(F) | ||
754 | AbstractArray(F::LU{T,Tridiagonal{T,V}}) where {T,V} = AbstractMatrix(F) | ||
755 | Matrix(F::LU{T,Tridiagonal{T,V}}) where {T,V} = Array(AbstractArray(F)) | ||
756 | Array(F::LU{T,Tridiagonal{T,V}}) where {T,V} = Matrix(F) |