Skip to content

Commit

Permalink
Avoid slow fallback Matrix constructor when converting Q from QR, (#2…
Browse files Browse the repository at this point in the history
…9865)

LQ and Hessenberg to Matrix

Fixes #29846
  • Loading branch information
andreasnoack authored Nov 2, 2018
1 parent 05fb47d commit 6a3a504
Show file tree
Hide file tree
Showing 3 changed files with 8 additions and 34 deletions.
34 changes: 2 additions & 32 deletions stdlib/LinearAlgebra/src/hessenberg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ true
hessenberg(A::StridedMatrix{T}) where T =
hessenberg!(copy_oftype(A, eigtype(T)))

struct HessenbergQ{T,S<:AbstractMatrix} <: AbstractMatrix{T}
struct HessenbergQ{T,S<:AbstractMatrix} <: AbstractQ{T}
factors::S
τ::Vector{T}
function HessenbergQ{T,S}(factors, τ) where {T,S<:AbstractMatrix}
Expand All @@ -77,8 +77,6 @@ struct HessenbergQ{T,S<:AbstractMatrix} <: AbstractMatrix{T}
end
HessenbergQ(factors::AbstractMatrix{T}, τ::Vector{T}) where {T} = HessenbergQ{T,typeof(factors)}(factors, τ)
HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.τ)
size(A::HessenbergQ, d) = size(A.factors, d)
size(A::HessenbergQ) = size(A.factors)

function getproperty(F::Hessenberg, d::Symbol)
d == :Q && return HessenbergQ(F)
Expand All @@ -89,17 +87,8 @@ end
Base.propertynames(F::Hessenberg, private::Bool=false) =
(:Q, :H, (private ? fieldnames(typeof(F)) : ())...)

function getindex(A::HessenbergQ, i::Integer, j::Integer)
x = zeros(eltype(A), size(A, 1))
x[i] = 1
y = zeros(eltype(A), size(A, 2))
y[j] = 1
return dot(x, lmul!(A, y))
end

## reconstruct the original matrix
Matrix(A::HessenbergQ{<:BlasFloat}) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.τ)
Array(A::HessenbergQ) = Matrix(A)
Matrix{T}(Q::HessenbergQ) where {T} = convert(Matrix{T}, LAPACK.orghr!(1, size(Q.factors, 1), copy(Q.factors), Q.τ))
AbstractMatrix(F::Hessenberg) = (fq = Array(F.Q); (fq * F.H) * fq')
AbstractArray(F::Hessenberg) = AbstractMatrix(F)
Matrix(F::Hessenberg) = Array(AbstractArray(F))
Expand All @@ -113,22 +102,3 @@ lmul!(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{T}) where {T<:Bl
(Q = adjQ.parent; LAPACK.ormhr!('L', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))
rmul!(X::StridedMatrix{T}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T<:BlasFloat} =
(Q = adjQ.parent; LAPACK.ormhr!('R', ifelse(T<:Real, 'T', 'C'), 1, size(Q.factors, 1), Q.factors, Q.τ, X))

function (*)(Q::HessenbergQ{T}, X::StridedVecOrMat{S}) where {T,S}
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
return lmul!(Q, copy_oftype(X, TT))
end
function (*)(X::StridedVecOrMat{S}, Q::HessenbergQ{T}) where {T,S}
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
return rmul!(copy_oftype(X, TT), Q)
end
function *(adjQ::Adjoint{<:Any,<:HessenbergQ{T}}, X::StridedVecOrMat{S}) where {T,S}
Q = adjQ.parent
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
return lmul!(adjoint(Q), copy_oftype(X, TT))
end
function *(X::StridedVecOrMat{S}, adjQ::Adjoint{<:Any,<:HessenbergQ{T}}) where {T,S}
Q = adjQ.parent
TT = typeof(zero(T)*zero(S) + zero(T)*zero(S))
return rmul!(copy_oftype(X, TT), adjoint(Q))
end
4 changes: 3 additions & 1 deletion stdlib/LinearAlgebra/src/lq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -113,7 +113,9 @@ end

LQPackedQ{T}(Q::LQPackedQ) where {T} = LQPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(Vector{T}, Q.τ))
AbstractMatrix{T}(Q::LQPackedQ) where {T} = LQPackedQ{T}(Q)
Matrix(A::LQPackedQ) = LAPACK.orglq!(copy(A.factors),A.τ)
Matrix{T}(A::LQPackedQ) where {T} = convert(Matrix{T}, LAPACK.orglq!(copy(A.factors),A.τ))
Matrix(A::LQPackedQ{T}) where {T} = Matrix{T}(A)
Array{T}(A::LQPackedQ{T}) where {T} = Matrix{T}(A)
Array(A::LQPackedQ) = Matrix(A)

size(F::LQ, dim::Integer) = size(getfield(F, :factors), dim)
Expand Down
4 changes: 3 additions & 1 deletion stdlib/LinearAlgebra/src/qr.jl
Original file line number Diff line number Diff line change
Expand Up @@ -500,7 +500,9 @@ AbstractMatrix{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q)
QRCompactWYQ{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ(convert(AbstractMatrix{S}, Q.factors), convert(AbstractMatrix{S}, Q.T))
AbstractMatrix{S}(Q::QRCompactWYQ{S}) where {S} = Q
AbstractMatrix{S}(Q::QRCompactWYQ) where {S} = QRCompactWYQ{S}(Q)
Matrix(A::AbstractQ{T}) where {T} = lmul!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...)))
Matrix{T}(A::AbstractQ) where {T} = lmul!(A, Matrix{T}(I, size(A.factors, 1), min(size(A.factors)...)))
Matrix(A::AbstractQ{T}) where {T} = Matrix{T}(A)
Array{T}(A::AbstractQ) where {T} = Matrix{T}(A)
Array(A::AbstractQ) = Matrix(A)

size(A::Union{QR,QRCompactWY,QRPivoted}, dim::Integer) = size(getfield(A, :factors), dim)
Expand Down

0 comments on commit 6a3a504

Please sign in to comment.