From 817c894a8436fff264604f31c8d52519daa26b6d Mon Sep 17 00:00:00 2001 From: Benjamin Cohen-Stead Date: Fri, 8 Nov 2024 15:41:41 -0500 Subject: [PATCH] v1.5.0 using Diagonal matrix --- Project.toml | 2 +- docs/src/developer_api.md | 7 -- src/LDR.jl | 3 +- src/developer_functions.jl | 162 ------------------------------------ src/exported_functions.jl | 105 ++++++++--------------- src/overloaded_functions.jl | 57 ++++++++----- 6 files changed, 77 insertions(+), 259 deletions(-) diff --git a/Project.toml b/Project.toml index b41e6e2..97efc9d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StableLinearAlgebra" uuid = "7d06c537-f8ff-4c22-91e1-ce4088e3cfd7" authors = ["Benjamin Cohen-Stead "] -version = "1.4.2" +version = "1.5.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/src/developer_api.md b/docs/src/developer_api.md index b13b5ce..4774703 100644 --- a/docs/src/developer_api.md +++ b/docs/src/developer_api.md @@ -1,13 +1,6 @@ # Developer API ```@docs -StableLinearAlgebra.det_D -StableLinearAlgebra.mul_D! -StableLinearAlgebra.div_D! -StableLinearAlgebra.lmul_D! -StableLinearAlgebra.rmul_D! -StableLinearAlgebra.ldiv_D! -StableLinearAlgebra.rdiv_D! StableLinearAlgebra.mul_P! StableLinearAlgebra.inv_P! StableLinearAlgebra.perm_sign diff --git a/src/LDR.jl b/src/LDR.jl index faa12d0..f6eabb3 100644 --- a/src/LDR.jl +++ b/src/LDR.jl @@ -175,7 +175,8 @@ function ldr!(F::LDR, ws::LDRWorkspace{T}) where{T} end # calculate R = D⁻¹⋅R - ldiv_D!(d, M) + D = Diagonal(d) + ldiv!(D, M) # calculate R⋅Pᵀ mul_P!(R, M, qr_ws.jpvt) diff --git a/src/developer_functions.jl b/src/developer_functions.jl index fbfdc34..1b22019 100644 --- a/src/developer_functions.jl +++ b/src/developer_functions.jl @@ -1,165 +1,3 @@ -@doc raw""" - det_D(d::AbstractVector{T}) where {T} - -Given a diagonal matrix ``D`` represented by the vector `d`, -return ``\textrm{sign}(\det D)`` and ``\log(|\det A|).`` -""" -function det_D(d::AbstractVector{T}) where {T} - - logdetD = zero(real(T)) - sgndetD = oneunit(T) - for i in eachindex(d) - logdetD += log(abs(d[i])) - sgndetD *= sign(d[i]) - end - - return logdetD, sgndetD -end - -@doc raw""" - mul_D!(A, d, B) - -Calculate the matrix product ``A = D \cdot B,`` where ``D`` is a diagonal matrix -represented by the vector `d`. -""" -function mul_D!(A::AbstractMatrix, d::AbstractVector, B::AbstractMatrix) - - @inbounds @fastmath for c in axes(B,2) - for r in eachindex(d) - A[r,c] = d[r] * B[r,c] - end - end - - return nothing -end - - -@doc raw""" - mul_D!(A, B, d) - -Calculate the matrix product ``A = B \cdot D,`` where ``D`` is a diagonal matrix -represented by the vector `d`. -""" -function mul_D!(A::AbstractMatrix, B::AbstractMatrix, d::AbstractVector) - - @inbounds @fastmath for c in eachindex(d) - for r in axes(A,1) - A[r,c] = B[r,c] * d[c] - end - end - - return nothing -end - - -@doc raw""" - div_D!(A, d, B) - -Calculate the matrix product ``A = D^{-1} \cdot B,`` where ``D`` is a diagonal matrix -represented by the vector `d`. -""" -function div_D!(A::AbstractMatrix, d::AbstractVector, B::AbstractMatrix) - - @inbounds @fastmath for c in axes(B,2) - for r in eachindex(d) - A[r,c] = B[r,c] / d[r] - end - end - - return nothing -end - - -@doc raw""" - div_D!(A, B, d) - -Calculate the matrix product ``A = B \cdot D^{-1},`` where ``D`` is a diagonal matrix -represented by the vector `d`. -""" -function div_D!(A::AbstractMatrix, B::AbstractMatrix, d::AbstractVector) - - @inbounds @fastmath for c in eachindex(d) - for r in axes(B,1) - A[r,c] = B[r,c] / d[c] - end - end - - return nothing -end - - -@doc raw""" - lmul_D!(d, M) - -In-place calculation of the matrix product ``M = D \cdot M,`` where ``D`` is a diagonal -matrix represented by the vector `d`. -""" -function lmul_D!(d::AbstractVector, M::AbstractMatrix) - - @inbounds @fastmath for c in axes(M,2) - for r in eachindex(d) - M[r,c] *= d[r] - end - end - - return nothing -end - - -@doc raw""" - rmul_D!(M, d) - -In-place calculation of the matrix product ``M = M \cdot D,`` where ``D`` is a diagonal -matrix represented by the vector `d`. -""" -function rmul_D!(M::AbstractMatrix, d::AbstractVector) - - @inbounds @fastmath for c in eachindex(d) - for r in axes(M,1) - M[r,c] *= d[c] - end - end - - return nothing -end - - -@doc raw""" - ldiv_D!(d, M) - -In-place calculation of the matrix product ``M = D^{-1} \cdot M,`` where ``D`` is a diagonal -matrix represented by the vector `d`. -""" -function ldiv_D!(d::AbstractVector, M::AbstractMatrix) - - @inbounds @fastmath for c in axes(M,2) - for r in eachindex(d) - M[r,c] /= d[r] - end - end - - return nothing -end - - -@doc raw""" - rdiv_D!(M, d) - -In-place calculation of the matrix product ``M = M \cdot D^{-1},`` where ``D`` is a diagonal -matrix represented by the vector `d`. -""" -function rdiv_D!(M::AbstractMatrix, d::AbstractVector) - - @inbounds @fastmath for c in eachindex(d) - for r in axes(M,1) - M[r,c] /= d[c] - end - end - - return nothing -end - - @doc raw""" mul_P!(A, p, B) diff --git a/src/exported_functions.jl b/src/exported_functions.jl index 7ab9035..61c3f5d 100644 --- a/src/exported_functions.jl +++ b/src/exported_functions.jl @@ -1,50 +1,3 @@ -# @doc raw""" -# adjoint_inv_ldr!(A⁻ᵀ::LDR{T}, A::AbstractMatrix{T}, ws::LDRWorkspace{T}) where {T} - -# Given a matrix ``A,`` calculate the [`LDR`](@ref) factorization ``(A^{-1})^{\dagger},`` -# storing the result in `A⁻ᵀ`. - -# # Algorithm - -# By calculating the pivoted QR factorization of `A`, we can calculate the [`LDR`](@ref) -# facorization `(A^{-1})^{\dagger} = L_0 D_0 R_0` using the procedure - -# ```math -# \begin{align*} -# (A^{-1})^{\dagger}:= & ([QRP^{T}]^{-1})^{\dagger}\\ -# = & ([Q\,\textrm{|diag}(R)|\,|\textrm{diag}(R)|^{-1}\,RP^{T}]^{-1})^{\dagger}\\ -# = & [\overset{R_{0}^{\dagger}}{\overbrace{PR^{-1}|\textrm{diag}(R)}|}\,\overset{D_{0}}{\overbrace{|\textrm{diag}(R)|^{-1}}}\,\overset{L_{0}^{\dagger}}{\overbrace{Q^{\dagger}}}]^{\dagger}\\ -# = & [R_{0}^{\dagger}D_{0}L_{0}^{\dagger}]^{\dagger}\\ -# = & L_{0}D_{0}R_{0}. -# \end{align*} -# ``` -# """ -# function adjoint_inv_ldr!(A⁻ᵀ::LDR{T}, A::AbstractMatrix{T}, ws::LDRWorkspace{T}) where {T} - -# # calculate A = Q⋅R⋅Pᵀ -# copyto!(A⁻ᵀ.L, A) -# LAPACK.geqp3!(A⁻ᵀ.L, ws.qr_ws) -# copyto!(A⁻ᵀ.R, A⁻ᵀ.L) -# triu!(A⁻ᵀ.R) # R (set lower triangular to 0) -# pᵀ = ws.qr_ws.jpvt -# p = ws.lu_ws.ipiv -# inv_P!(p, pᵀ) # P -# LAPACK.orgqr!(Aᵀ.L, ws.qr_ws) # L₀ = Q -# @fastmath @inbounds for i in eachindex(A⁻ᵀ.d) -# A⁻ᵀ.d[i] = inv(abs(A⁻ᵀ.R[i,i])) # D₀ = |diag(R)|⁻¹ -# end -# R′ = A⁻ᵀ.R -# lmul_D!(A⁻ᵀ.d, R′) # R′ := |diag(R)|⁻¹⋅R -# R″ = R′ -# LinearAlgebra.inv!(UpperTriangular(R″)) # R″ = (R′)⁻¹ = R⁻¹⋅|diag(R)| -# R₀ᵀ = ws.M -# mul_P!(R₀ᵀ, p, R″) # R₀ᵀ = P⋅R″ = P⋅R⁻¹⋅diag(R) -# adjoint(A⁻ᵀ.R, R₀ᵀ) # R₀ = [P⋅R⁻¹⋅diag(R)]ᵀ - -# return nothing -# end - - @doc raw""" inv_IpA!(G::AbstractMatrix{T}, A::LDR{T,E}, ws::LDRWorkspace{T,E})::Tuple{E,T} where {T,E} @@ -85,18 +38,20 @@ function inv_IpA!(G::AbstractMatrix{T}, A::LDR{T,E}, ws::LDRWorkspace{T,E})::Tup @. d₋ = min(dₐ, 1) # calculate Lₐ⋅D₋ - mul_D!(ws.M, Lₐ, d₋) + D₋ = Diagonal(d₋) + mul!(ws.M, Lₐ, D₋) # calculate D₊ = max(Dₐ, 1) d₊ = ws.v @. d₊ = max(dₐ, 1) # calculate sign(det(D₊)) and log(|det(D₊)|) - logdetD₊, sgndetD₊ = det_D(d₊) + D₊ = Diagonal(d₊) + logdetD₊, sgndetD₊ = logabsdet(D₊) # calculate Rₐ⁻¹⋅D₊⁻¹ Rₐ⁻¹D₊ = Rₐ⁻¹ - rdiv_D!(Rₐ⁻¹D₊, d₊) + rdiv!(Rₐ⁻¹D₊, D₊) # calculate M = Rₐ⁻¹⋅D₊⁻¹ + Lₐ⋅D₋ axpy!(1.0, Rₐ⁻¹D₊, ws.M) @@ -160,40 +115,44 @@ function inv_IpUV!(G::AbstractMatrix{T}, U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorksp # calcuate Dᵥ₊ = max(Dᵥ, 1) dᵥ₊ = ws.v @. dᵥ₊ = max(dᵥ, 1) + Dᵥ₊ = Diagonal(dᵥ₊) # calculate sign(det(Dᵥ₊)) and log(|det(Dᵥ₊)|) - logdetDᵥ₊, sgndetDᵥ₊ = det_D(dᵥ₊) + logdetDᵥ₊, sgndetDᵥ₊ = logabsdet(Dᵥ₊) # calculate Rᵥ⁻¹⋅Dᵥ₊⁻¹ - rdiv_D!(Rᵥ⁻¹, dᵥ₊) + rdiv!(Rᵥ⁻¹, Dᵥ₊) Rᵥ⁻¹Dᵥ₊⁻¹ = Rᵥ⁻¹ # calcuate Dᵤ₊ = max(Dᵤ, 1) dᵤ₊ = ws.v @. dᵤ₊ = max(dᵤ, 1) + Dᵤ₊ = Diagonal(dᵤ₊) # calculate sign(det(Dᵥ₊)) and log(|det(Dᵥ₊)|) - logdetDᵤ₊, sgndetDᵤ₊ = det_D(dᵤ₊) + logdetDᵤ₊, sgndetDᵤ₊ = logabsdet(Dᵤ₊) # calcualte Dᵤ₊⁻¹⋅Lᵤᵀ adjoint!(ws.M, Lᵤ) - ldiv_D!(dᵤ₊, ws.M) + ldiv!(Dᵤ₊, ws.M) Dᵤ₊⁻¹Lᵤᵀ = ws.M # calculate Dᵤ₋ = min(Dᵤ, 1) dᵤ₋ = ws.v @. dᵤ₋ = min(dᵤ, 1) + Dᵤ₋ = Diagonal(dᵤ₋) # calculate Dᵤ₋⋅Rᵤ⋅Lᵥ mul!(G, Rᵤ, Lᵥ) # Rᵤ⋅Lᵥ - lmul_D!(dᵤ₋, G) # Dᵤ₋⋅Rᵤ⋅Lᵥ + lmul!(Dᵤ₋, G) # Dᵤ₋⋅Rᵤ⋅Lᵥ # calculate Dᵥ₋ = min(Dᵥ, 1) dᵥ₋ = ws.v @. dᵥ₋ = min(dᵥ, 1) + Dᵥ₋ = Diagonal(dᵥ₋) # caluclate Dᵤ₋⋅Rᵤ⋅Lᵥ⋅Dᵥ₋ - rmul_D!(G, dᵥ₋) + rmul!(G, Dᵥ₋) # caluclate Dᵤ₊⁻¹⋅Lᵤᵀ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹ mul!(ws.M″, Dᵤ₊⁻¹Lᵤᵀ, Rᵥ⁻¹Dᵥ₊⁻¹) @@ -270,41 +229,45 @@ function inv_UpV!(G::AbstractMatrix{T}, U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorkspa # calcuate Dᵥ₊ = max(Dᵥ, 1) dᵥ₊ = ws.v @. dᵥ₊ = max(dᵥ, 1) + Dᵥ₊ = Diagonal(dᵥ₊) # calculate sign(det(Dᵥ₊)) and log(|det(Dᵥ₊)|) - logdetDᵥ₊, sgndetDᵥ₊ = det_D(dᵥ₊) + logdetDᵥ₊, sgndetDᵥ₊ = logabsdet(Dᵥ₊) # calculate Rᵥ⁻¹⋅Dᵥ₊⁻¹ - rdiv_D!(Rᵥ⁻¹, dᵥ₊) + rdiv!(Rᵥ⁻¹, Dᵥ₊) Rᵥ⁻¹Dᵥ₊⁻¹ = Rᵥ⁻¹ # calcuate Dᵤ₊ = max(Dᵤ, 1) dᵤ₊ = ws.v @. dᵤ₊ = max(dᵤ, 1) + Dᵤ₊ = Diagonal(dᵤ₊) - # calculate sign(det(Dᵥ₊)) and log(|det(Dᵥ₊)|) - logdetDᵤ₊, sgndetDᵤ₊ = det_D(dᵤ₊) + # calculate sign(det(Dᵤ₊)) and log(|det(Dᵤ₊)|) + logdetDᵤ₊, sgndetDᵤ₊ = logabsdet(Dᵤ₊) # calcualte Dᵤ₊⁻¹⋅Lᵤᵀ adjoint!(ws.M, Lᵤ) - ldiv_D!(dᵤ₊, ws.M) + ldiv!(Dᵤ₊, ws.M) Dᵤ₊⁻¹Lᵤᵀ = ws.M # calculate Dᵤ₋ = min(Dᵤ, 1) dᵤ₋ = ws.v @. dᵤ₋ = min(dᵤ, 1) + Dᵤ₋ = Diagonal(dᵤ₋) # calculate Dᵤ₋⋅Rᵤ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹ mul!(G, Rᵤ, Rᵥ⁻¹Dᵥ₊⁻¹) # Rᵤ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹ - lmul_D!(dᵤ₋, G) # Dᵤ₋⋅[Rᵤ⋅Rᵥ⁻¹⋅Dᵥ₊] + lmul!(Dᵤ₋, G) # Dᵤ₋⋅[Rᵤ⋅Rᵥ⁻¹⋅Dᵥ₊] # calculate Dᵥ₋ = min(Dᵥ, 1) dᵥ₋ = ws.v @. dᵥ₋ = min(dᵥ, 1) + Dᵥ₋ = Diagonal(dᵥ₋) # calculate Dᵤ₊⁻¹⋅Lᵤᵀ⋅Lᵥ⋅Dᵥ₋ mul!(ws.M″, Dᵤ₊⁻¹Lᵤᵀ, Lᵥ) - rmul_D!(ws.M″, dᵥ₋) + rmul!(ws.M″, Dᵥ₋) # calculate M = Dᵤ₋⋅Rᵤ⋅Rᵥ⁻¹⋅Dᵥ₊ + Dᵤ₊⁻¹⋅Lᵤᵀ⋅Lᵥ⋅Dᵥ₋ M = G @@ -374,14 +337,15 @@ function inv_invUpV!(G::AbstractMatrix{T}, U::LDR{T,E}, V::LDR{T,E}, ws::LDRWork # calculate Dᵤ₋ = min(Dᵤ, 1) dᵤ₋ = ws.v @. dᵤ₋ = min(dᵤ, 1) + Dᵤ₋ = Diagonal(dᵤ₋) # calculate sign(det(Dᵤ₋)) and log(|det(Dᵤ₋)|) - logdetDᵤ₋, sgndetDᵤ₋ = det_D(dᵤ₋) + logdetDᵤ₋, sgndetDᵤ₋ = logabsdet(Dᵤ₋) # calculate Dᵤ₋⋅Rᵤ Dᵤ₋Rᵤ = ws.M′ copyto!(Dᵤ₋Rᵤ, Rᵤ) - lmul_D!(dᵤ₋, Dᵤ₋Rᵤ) + lmul!(Dᵤ₋, Dᵤ₋Rᵤ) # calculate Rᵥ⁻¹, sign(det(Rᵥ⁻¹)) and log(|det(Rᵥ⁻¹)|) Rᵥ⁻¹ = ws.M″ @@ -391,30 +355,33 @@ function inv_invUpV!(G::AbstractMatrix{T}, U::LDR{T,E}, V::LDR{T,E}, ws::LDRWork # calculate Dᵥ₊ = max(Dᵥ, 1) dᵥ₊ = ws.v @. dᵥ₊ = max(dᵥ, 1) + Dᵥ₊ = Diagonal(dᵥ₊) # calculate sign(det(Dᵥ₊)) and log(|det(Dᵥ₊)|) - logdetDᵥ₊, sgndetDᵥ₊ = det_D(dᵥ₊) + logdetDᵥ₊, sgndetDᵥ₊ = logabsdet(Dᵥ₊) # calculate Rᵥ⁻¹⋅Dᵥ₊⁻¹ Rᵥ⁻¹Dᵥ₊⁻¹ = Rᵥ⁻¹ - rdiv_D!(Rᵥ⁻¹Dᵥ₊⁻¹, dᵥ₊) + rdiv!(Rᵥ⁻¹Dᵥ₊⁻¹, Dᵥ₊) # calculate Dᵥ₋ = min(Dᵥ, 1) dᵥ₋ = ws.v @. dᵥ₋ = min(dᵥ, 1) + Dᵥ₋ = Diagonal(dᵥ₋) # calculate Dᵤ₋⋅Rᵤ⋅Lᵥ⋅Dᵥ₋ mul!(G, Dᵤ₋Rᵤ, Lᵥ) # Dᵤ₋⋅Rᵤ⋅Lᵥ - rmul_D!(G, dᵥ₋) # [Dᵤ₋⋅Rᵤ⋅Lᵥ]⋅Dᵥ₋ + rmul!(G, Dᵥ₋) # [Dᵤ₋⋅Rᵤ⋅Lᵥ]⋅Dᵥ₋ # calculate Dᵤ₊ = max(Dᵤ, 1) dᵤ₊ = ws.v @. dᵤ₊ = max(dᵤ, 1) + Dᵤ₊ = Diagonal(dᵤ₊) # calculate Dᵤ₊⁻¹⋅Lᵤᵀ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹ Lᵤᵀ = adjoint(Lᵤ) mul!(ws.M, Lᵤᵀ, Rᵥ⁻¹Dᵥ₊⁻¹) # Lᵤᵀ⋅[Rᵥ⁻¹⋅Dᵥ₊⁻¹] - ldiv_D!(dᵤ₊, ws.M) # Dᵤ₊⁻¹⋅[Lᵤᵀ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹] + ldiv!(Dᵤ₊, ws.M) # Dᵤ₊⁻¹⋅[Lᵤᵀ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹] # calculate Dᵤ₊⁻¹⋅Lᵤᵀ⋅Rᵥ⁻¹⋅Dᵥ₊⁻¹ + Dᵤ₋⋅Rᵤ⋅Lᵥ⋅Dᵥ₋ axpy!(1.0, ws.M, G) diff --git a/src/overloaded_functions.jl b/src/overloaded_functions.jl index 5557098..fa2f16c 100644 --- a/src/overloaded_functions.jl +++ b/src/overloaded_functions.jl @@ -37,7 +37,8 @@ function copyto!(U::AbstractMatrix{T}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where (; M) = ws copyto!(M, R) # R - lmul_D!(d, M) # D⋅R + D = Diagonal(d) + lmul!(D, M) # D⋅R mul!(U, L, M) # U = L⋅D⋅R return nothing @@ -77,7 +78,8 @@ function adjoint!(Aᵀ::AbstractMatrix{T}, A::LDR{T,E}, ws::LDRWorkspace{T,E}) w Rᵀ = ws.M′ adjoint!(Rᵀ, R) adjoint!(ws.M, L) # Lᵀ - lmul_D!(d, ws.M) # D⋅Lᵀ + D = Diagonal(d) + lmul!(D, ws.M) # D⋅Lᵀ mul!(Aᵀ, Rᵀ, ws.M) # Rᵀ⋅D⋅Lᵀ return nothing @@ -96,7 +98,8 @@ function lmul!(U::LDR{T,E}, V::AbstractMatrix{T}, ws::LDRWorkspace{T,E}) where { # calculate V := Lᵤ⋅Dᵤ⋅Rᵤ⋅V mul!(ws.M, U.R, V) # Rᵤ⋅V - lmul_D!(U.d, ws.M) # Dᵤ⋅Rᵤ⋅V + Dᵤ = Diagonal(U.d) + lmul!(Dᵤ, ws.M) # Dᵤ⋅Rᵤ⋅V mul!(V, U.L, ws.M) # V := Lᵤ⋅Dᵤ⋅Rᵤ⋅V return nothing @@ -128,7 +131,8 @@ function lmul!(U::AbstractMatrix{T}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where { # calculate product U⋅Lᵥ⋅Dᵥ mul!(ws.M, U, V.L) # U⋅Lᵥ - mul_D!(V.L, ws.M, V.d) # U⋅Lᵥ⋅Dᵥ + Dᵥ = Diagonal(V.d) + mul!(V.L, ws.M, Dᵥ) # U⋅Lᵥ⋅Dᵥ # calcualte [L₀⋅D₀⋅R₀] = U⋅Lᵥ⋅Dᵥ ldr!(V, ws) @@ -169,8 +173,10 @@ function lmul!(U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where {T,E} mul!(ws.M, U.R, V.L) # calculate Dᵤ⋅M⋅Dᵥ - rmul_D!(ws.M, V.d) # M⋅Dᵥ - mul_D!(V.L, U.d, ws.M) # Dᵤ⋅M⋅Dᵥ + Dᵥ = Diagonal(V.d) + Dᵤ = Diagonal(U.d) + rmul!(ws.M, Dᵥ) # M⋅Dᵥ + mul!(V.L, Dᵤ, ws.M) # Dᵤ⋅M⋅Dᵥ # calculate [L₀⋅D₀⋅R₀] = Dᵤ⋅M⋅Dᵥ ldr!(V, ws) @@ -200,7 +206,8 @@ function rmul!(U::AbstractMatrix{T}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where { # calculate U := U⋅Lᵥ⋅Dᵥ⋅Rᵥ mul!(ws.M, U, V.L) # U⋅Lᵥ - rmul_D!(ws.M, V.d) # U⋅Lᵥ⋅Dᵥ + Dᵥ = Diagonal(V.d) + rmul!(ws.M, Dᵥ) # U⋅Lᵥ⋅Dᵥ mul!(U, ws.M, V.R) # U := U⋅Lᵥ⋅Dᵥ⋅Rᵥ return nothing @@ -232,7 +239,8 @@ function rmul!(U::LDR{T,E}, V::AbstractMatrix{T}, ws::LDRWorkspace{T,E}) where { # calculate Dᵤ⋅Rᵤ⋅V mul!(U.L, U.R, V) - lmul_D!(U.d, U.L) + Dᵤ = Diagonal(U.d) + lmul!(Dᵤ, U.L) # calculate [L₀⋅D₀⋅R₀] = Dᵤ⋅Rᵤ⋅V ldr!(U, ws) @@ -273,8 +281,10 @@ function rmul!(U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where {T,E} mul!(ws.M, U.R, V.L) # calculate Dᵤ⋅Rᵤ⋅Lᵥ⋅Dᵥ - rmul_D!(ws.M, V.d) - mul_D!(U.L, U.d, ws.M) + Dᵥ = Diagonal(V.d) + Dᵤ = Diagonal(U.d) + rmul!(ws.M, Dᵥ) + mul!(U.L, Dᵤ, ws.M) # calculate [L₀⋅D₀⋅R₀] = Dᵤ⋅Rᵤ⋅Lᵥ⋅Dᵥ ldr!(U, ws) @@ -381,7 +391,8 @@ function ldiv!(U::LDR{T,E}, V::AbstractMatrix{T}, ws::LDRWorkspace{T,E}) where { Lᵤ = ws.M copyto!(Lᵤ, U.L) ldiv_lu!(Lᵤ, V, ws.lu_ws) # Lᵤ⁻¹⋅V - ldiv_D!(U.d, V) # Dᵤ⁻¹⋅Lᵤ⁻¹⋅V + Dᵤ = Diagonal(U.d) + ldiv!(Dᵤ, V) # Dᵤ⁻¹⋅Lᵤ⁻¹⋅V Rᵤ = ws.M copyto!(Rᵤ, U.R) ldiv_lu!(Rᵤ, V, ws.lu_ws) # V := Rᵤ⁻¹⋅Dᵤ⁻¹⋅Lᵤ⁻¹⋅V @@ -435,8 +446,10 @@ function ldiv!(U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where {T,E} copyto!(Rᵥ, V.R) # calculate Rᵤ⁻¹⋅Dᵤ⁻¹⋅Lᵤᵀ⋅Lᵥ⋅Dᵥ - rmul_D!(V.L, V.d) # [Lᵤᵀ⋅Lᵥ]⋅Dᵥ - ldiv_D!(U.d, V.L) # Dᵤ⁻¹⋅[Lᵤᵀ⋅Lᵥ⋅Dᵥ] + Dᵥ = Diagonal(V.d) + Dᵤ = Diagonal(U.d) + rmul!(V.L, Dᵥ) # [Lᵤᵀ⋅Lᵥ]⋅Dᵥ + ldiv!(Dᵤ, V.L) # Dᵤ⁻¹⋅[Lᵤᵀ⋅Lᵥ⋅Dᵥ] Rᵤ = ws.M copyto!(Rᵤ, U.R) ldiv_lu!(Rᵤ, V.L, ws.lu_ws) # Rᵤ⁻¹⋅[Dᵤ⁻¹⋅Lᵤᵀ⋅Lᵥ⋅Dᵥ] @@ -493,7 +506,8 @@ function ldiv!(U::AbstractMatrix{T}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where { copyto!(Rᵥ, V.R) # calculate U⁻¹⋅Lᵥ⋅Dᵥ - rmul_D!(V.L, V.d) # Lᵥ⋅Dᵥ + Dᵥ = Diagonal(V.d) + rmul!(V.L, Dᵥ) # Lᵥ⋅Dᵥ copyto!(ws.M, U) ldiv_lu!(ws.M, V.L, ws.lu_ws) # U⁻¹⋅Lᵥ⋅Dᵥ @@ -539,7 +553,8 @@ function rdiv!(U::AbstractMatrix{T}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where { Lᵥ = ws.M′ copyto!(Lᵥ, V.L) ldiv_lu!(Lᵥ, ws.M, ws.lu_ws) # Lᵥ⁻¹ - ldiv_D!(V.d, ws.M) # Dᵥ⁻¹⋅Lᵥ⁻¹ + Dᵥ = Diagonal(V.d) + ldiv!(Dᵥ, ws.M) # Dᵥ⁻¹⋅Lᵥ⁻¹ Rᵥ = ws.M′ copyto!(Rᵥ, V.R) ldiv_lu!(Rᵥ, ws.M, ws.lu_ws) # Rᵥ⁻¹⋅Dᵥ⁻¹⋅Lᵥ⁻¹ @@ -600,8 +615,10 @@ function rdiv!(U::LDR{T,E}, V::LDR{T,E}, ws::LDRWorkspace{T,E}) where {T,E} copyto!(Lᵤ, U.L) # calculate Dᵤ⋅M⋅Dᵥ⁻¹ - div_D!(U.L, ws.M, V.d) - lmul_D!(U.d, U.L) + Dᵥ = Diagonal(V.d) + Dᵤ = Diagonal(U.d) + mul!(U.L, Dᵤ, ws.M) + rdiv!(U.L, Dᵥ) # calculate [L₀⋅D₀⋅R₀] = Dᵤ⋅M⋅Dᵥ⁻¹ ldr!(U, ws) @@ -660,10 +677,11 @@ function rdiv!(U::LDR{T,E}, V::AbstractMatrix{T}, ws::LDRWorkspace{T,E}) where { copyto!(Lᵤ, U.L) # calculate Dᵤ⋅Rᵤ⋅V⁻¹ + Dᵤ = Diagonal(U.d) copyto!(ws.M, V) inv_lu!(ws.M, ws.lu_ws) # V⁻¹ mul!(U.L, U.R, ws.M) # Rᵤ⋅V⁻¹ - lmul_D!(U.d, U.L) # Dᵤ⋅Rᵤ⋅V⁻¹ + lmul!(Dᵤ, U.L) # Dᵤ⋅Rᵤ⋅V⁻¹ # calcualte [L₀⋅D₀⋅R₀] = Dᵤ⋅Rᵤ⋅V⁻¹ ldr!(U, ws) @@ -704,7 +722,8 @@ function logabsdet(A::LDR{T,E}, ws::LDRWorkspace{T,E}) where {T,E} copyto!(ws.M, A.L) logdetL, sgndetL = det_lu!(ws.M, ws.lu_ws) - logdetD, sgndetD = det_D(A.d) + D = Diagonal(A.d) + logdetD, sgndetD = logabsdet(D) copyto!(ws.M, A.R) logdetR, sgndetR = det_lu!(ws.M, ws.lu_ws) logdetA = logdetL + logdetD + logdetR