|
| 1 | +""" |
| 2 | + LULinv <: Factorization |
| 3 | +
|
| 4 | +Matrix factorization type of the `LUL⁻¹` factorization of a square matrix `A`. This |
| 5 | +is the return type of [`lulinv`](@ref), the corresponding matrix factorization function. |
| 6 | +
|
| 7 | +The individual components of the factorization `F::LULinv` can be accessed via [`getproperty`](@ref): |
| 8 | +
|
| 9 | +| Component | Description | |
| 10 | +|:----------|:--------------------------------------------| |
| 11 | +| `F.L` | `L` (unit lower triangular) part of `LUL⁻¹` | |
| 12 | +| `F.U` | `U` (upper triangular) part of `LUL⁻¹` | |
| 13 | +
|
| 14 | +Iterating the factorization produces the components `F.L` and `F.U`. |
| 15 | +
|
| 16 | +# Examples |
| 17 | +```jldoctest |
| 18 | +julia> A = [4 3; 6 3] |
| 19 | +2×2 Array{Int64,2}: |
| 20 | + 4 3 |
| 21 | + 6 3 |
| 22 | +
|
| 23 | +julia> F = lulinv(A) |
| 24 | +LULinv{Float64, Matrix{Float64}} |
| 25 | +L factor: |
| 26 | +2×2 Matrix{Float64}: |
| 27 | + 1.0 0.0 |
| 28 | + -1.59067 1.0 |
| 29 | +U factor: |
| 30 | +2×2 Matrix{Float64}: |
| 31 | + -0.772002 3.0 |
| 32 | + 0.0 7.772 |
| 33 | +
|
| 34 | +julia> F.L * F.U / F.L ≈ A |
| 35 | +true |
| 36 | +
|
| 37 | +julia> l, u = lulinv(A); # destructuring via iteration |
| 38 | +
|
| 39 | +julia> l == F.L && u == F.U |
| 40 | +true |
| 41 | +
|
| 42 | +julia> A = [-150 334 778; -89 195 464; 5 -10 -27] |
| 43 | +3×3 Matrix{Int64}: |
| 44 | + -150 334 778 |
| 45 | + -89 195 464 |
| 46 | + 5 -10 -27 |
| 47 | +
|
| 48 | +julia> F = lulinv(A, [17, -2, 3//1]) # can input rational eigenvalues directly |
| 49 | +LULinv{Rational{Int64}, Matrix{Rational{Int64}}} |
| 50 | +L factor: |
| 51 | +3×3 Matrix{Rational{Int64}}: |
| 52 | + 1 0 0 |
| 53 | + 1//2 1 0 |
| 54 | + 0 -2//5 1 |
| 55 | +U factor: |
| 56 | +3×3 Matrix{Rational{Int64}}: |
| 57 | + 17 114//5 778 |
| 58 | + 0 -2 75 |
| 59 | + 0 0 3 |
| 60 | +
|
| 61 | +julia> F.L * F.U / F.L == A |
| 62 | +true |
| 63 | +``` |
| 64 | +""" |
| 65 | +struct LULinv{T, S <: AbstractMatrix{T}} <: Factorization{T} |
| 66 | + factors::S |
| 67 | + function LULinv{T, S}(factors) where {T, S <: AbstractMatrix{T}} |
| 68 | + require_one_based_indexing(factors) |
| 69 | + new{T, S}(factors) |
| 70 | + end |
| 71 | +end |
| 72 | + |
| 73 | + |
| 74 | +LULinv(factors::AbstractMatrix{T}) where T = LULinv{T, typeof(factors)}(factors) |
| 75 | +LULinv{T}(factors::AbstractMatrix) where T = LULinv(convert(AbstractMatrix{T}, factors)) |
| 76 | +LULinv{T}(F::LULinv) where T = LULinv{T}(F.factors) |
| 77 | + |
| 78 | +iterate(F::LULinv) = (F.L, Val(:U)) |
| 79 | +iterate(F::LULinv, ::Val{:U}) = (F.U, Val(:done)) |
| 80 | +iterate(F::LULinv, ::Val{:done}) = nothing |
| 81 | + |
| 82 | + |
| 83 | +function lulinvtype(T::Type) |
| 84 | + # In generic_ulfact!, the elements of the lower part of the matrix are |
| 85 | + # obtained using the division of two matrix elements. Hence their type can |
| 86 | + # be different (e.g. the division of two types with the same unit is a type |
| 87 | + # without unit). |
| 88 | + # The elements of the upper part are obtained by U - L * U / L |
| 89 | + # where U is an upper part element and L is a lower part element. |
| 90 | + # Therefore, the types LT, UT should be invariant under the map: |
| 91 | + # (LT, UT) -> begin |
| 92 | + # L = oneunit(UT) / oneunit(UT) |
| 93 | + # U = oneunit(UT) - L * oneunit(UT) / L |
| 94 | + # typeof(L), typeof(U) |
| 95 | + # end |
| 96 | + # The following should handle most cases |
| 97 | + UT = typeof(oneunit(T) - (oneunit(T) / (oneunit(T) + zero(T)) * oneunit(T) * (oneunit(T) + zero(T)) / oneunit(T))) |
| 98 | + LT = typeof(oneunit(UT) / oneunit(UT)) |
| 99 | + S = promote_type(T, LT, UT) |
| 100 | +end |
| 101 | + |
| 102 | + |
| 103 | +size(A::LULinv) = size(getfield(A, :factors)) |
| 104 | +size(A::LULinv, i) = size(getfield(A, :factors), i) |
| 105 | + |
| 106 | +function getL(F::LULinv{T}) where T |
| 107 | + n = size(F.factors, 1) |
| 108 | + L = tril!(getindex(getfield(F, :factors), 1:n, 1:n)) |
| 109 | + for i in 1:n L[i, i] = one(T) end |
| 110 | + return L |
| 111 | +end |
| 112 | + |
| 113 | +function getU(F::LULinv) |
| 114 | + n = size(F.factors, 1) |
| 115 | + triu!(getindex(getfield(F, :factors), 1:n, 1:n)) |
| 116 | +end |
| 117 | + |
| 118 | +function getproperty(F::LULinv{T, <: AbstractMatrix}, d::Symbol) where T |
| 119 | + if d === :L |
| 120 | + return getL(F) |
| 121 | + elseif d === :U |
| 122 | + return getU(F) |
| 123 | + else |
| 124 | + getfield(F, d) |
| 125 | + end |
| 126 | +end |
| 127 | + |
| 128 | +propertynames(F::LULinv, private::Bool=false) = |
| 129 | + (:L, :U, (private ? fieldnames(typeof(F)) : ())...) |
| 130 | + |
| 131 | +function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LULinv) |
| 132 | + summary(io, F); println(io) |
| 133 | + println(io, "L factor:") |
| 134 | + show(io, mime, F.L) |
| 135 | + println(io, "\nU factor:") |
| 136 | + show(io, mime, F.U) |
| 137 | +end |
| 138 | + |
| 139 | + |
| 140 | +lulinv(A::AbstractMatrix; kwds...) = lulinv(A, eigvals(A); kwds...) |
| 141 | +function lulinv(A::AbstractMatrix{T}, λ::AbstractVector{T}; kwds...) where T |
| 142 | + S = lulinvtype(T) |
| 143 | + lulinv!(copy_oftype(A, S), copy_oftype(λ, S); kwds...) |
| 144 | +end |
| 145 | +function lulinv(A::AbstractMatrix{T1}, λ::AbstractVector{T2}; kwds...) where {T1, T2} |
| 146 | + T = promote_type(T1, T2) |
| 147 | + S = lulinvtype(T) |
| 148 | + lulinv!(copy_oftype(A, S), copy_oftype(λ, S); kwds...) |
| 149 | +end |
| 150 | + |
| 151 | +function lulinv!(A::Matrix{T}, λ::Vector{T}; rtol::Real = size(A, 1)*eps(real(float(oneunit(T))))) where T |
| 152 | + n = checksquare(A) |
| 153 | + n == length(λ) || throw(ArgumentError("Eigenvalue count does not match matrix dimensions.")) |
| 154 | + v = zeros(T, n) |
| 155 | + for i in 1:n-1 |
| 156 | + # We must find an eigenvector with nonzero "first" entry. A failed UL factorization of A-λI reveals this vector provided L₁₁ == 0. |
| 157 | + for j in 1:length(λ) |
| 158 | + F = ul!(view(A, i:n, i:n) - λ[j]*I; check=false) |
| 159 | + nrm = norm(F.L) |
| 160 | + if norm(F.L[1]) ≤ rtol*nrm # v[i] is a free parameter; we set it to 1. |
| 161 | + fill!(v, zero(T)) |
| 162 | + v[i] = one(T) |
| 163 | + # Next, we must scan the remaining free parameters, set them to 0, so that we find a nonsingular lower-triangular linear system for the nontrivial remaining part of the eigenvector. |
| 164 | + idx = Int[] |
| 165 | + for k in 2:n+1-i |
| 166 | + if norm(F.L[k, k]) > rtol*nrm |
| 167 | + push!(idx, k) |
| 168 | + end |
| 169 | + end |
| 170 | + if !isempty(idx) |
| 171 | + v[idx.+(i-1)] .= -F.L[idx, 1] |
| 172 | + ldiv!(LowerTriangular(view(F.L, idx, idx)), view(v, idx.+(i-1))) |
| 173 | + end |
| 174 | + deleteat!(λ, j) |
| 175 | + break |
| 176 | + end |
| 177 | + end |
| 178 | + for k in 1:n |
| 179 | + for j in i+1:n |
| 180 | + A[k, i] += A[k, j]*v[j] |
| 181 | + end |
| 182 | + end |
| 183 | + for j in i:n |
| 184 | + for k in i+1:n |
| 185 | + A[k, j] -= A[i, j]*v[k] |
| 186 | + end |
| 187 | + end |
| 188 | + for k in i+1:n |
| 189 | + A[k, i] = v[k] |
| 190 | + end |
| 191 | + end |
| 192 | + return LULinv(A) |
| 193 | +end |
| 194 | + |
| 195 | +function ldiv!(F::LULinv, B::AbstractVecOrMat) |
| 196 | + L = UnitLowerTriangular(F.factors) |
| 197 | + return lmul!(L, ldiv!(UpperTriangular(F.factors), ldiv!(L, B))) |
| 198 | +end |
| 199 | + |
| 200 | +function rdiv!(B::AbstractVecOrMat, F::LULinv) |
| 201 | + L = UnitLowerTriangular(F.factors) |
| 202 | + return rdiv!(rdiv!(rmul!(B, L), UpperTriangular(F.factors)), L) |
| 203 | +end |
| 204 | + |
| 205 | +det(F::LULinv) = det(UpperTriangular(F.factors)) |
| 206 | +logabsdet(F::LULinv) = logabsdet(UpperTriangular(F.factors)) |
0 commit comments