diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl index c90b621..04aebe1 100644 --- a/src/MatrixFactorizations.jl +++ b/src/MatrixFactorizations.jl @@ -1,6 +1,6 @@ module MatrixFactorizations using Base, LinearAlgebra, ArrayLayouts -import Base: axes, axes1, getproperty, iterate, tail, oneto +import Base: axes, axes1, getproperty, iterate, tail, oneto, BroadcastStyle import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!, copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym, det, logdet, logabsdet, isposdef @@ -50,12 +50,18 @@ const TransposeFact = isdefined(LinearAlgebra, :TransposeFactorization) ? Linear abstract type LayoutQ{T} <: AbstractQ{T} end +# Support orthogonal matrices as an abstract matrix +abstract type LayoutQMatrix{T} <: LayoutMatrix{T} end + +const LayoutQTypes{T} = Union{LayoutQ{T}, LayoutQMatrix{T}} + + @_layoutlmul LayoutQ @_layoutlmul AdjointQtype{<:Any,<:LayoutQ} @_layoutrmul LayoutQ @_layoutrmul AdjointQtype{<:Any,<:LayoutQ} -LinearAlgebra.copymutable(Q::LayoutQ) = copymutable_size(size(Q), Q) +LinearAlgebra.copymutable(Q::LayoutQTypes) = copymutable_size(size(Q), Q) copymutable_size(sz, Q) = lmul!(Q, Matrix{eltype(Q)}(I, sz)) (*)(Q::LayoutQ, b::AbstractVector) = _mul(Q, b) @@ -106,16 +112,17 @@ end *(A::AbstractTriangular, B::LayoutQ) = mul(A, B) *(A::AbstractTriangular, B::AdjointQtype{<:Any,<:LayoutQ}) = mul(A, B) -axes(Q::LayoutQ, dim::Integer) = axes(getfield(Q, :factors), dim == 2 ? 1 : dim) -axes(Q::LayoutQ) = axes(Q, 1), axes(Q, 2) -copy(Q::LayoutQ) = Q +axes(Q::LayoutQTypes, dim::Integer) = axes(getfield(Q, :factors), dim == 2 ? 1 : dim) +axes(Q::LayoutQTypes) = axes(Q, 1), axes(Q, 2) +copy(Q::LayoutQTypes) = Q Base.@propagate_inbounds getindex(A::LayoutQ, I...) = layout_getindex(A, I...) +Base.@propagate_inbounds getindex(A::LayoutQMatrix, k::Int, j::Int) = AbstractQ(A)[k, j] # by default, fall back to AbstractQ methods layout_getindex(A::LayoutQ, I...) = Base.invoke(Base.getindex, Tuple{AbstractQ, typeof.(I)...}, A, I...) -size(Q::LayoutQ, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim) -size(Q::LayoutQ) = size(Q, 1), size(Q, 2) +size(Q::LayoutQTypes, dim::Integer) = size(getfield(Q, :factors), dim == 2 ? 1 : dim) +size(Q::LayoutQTypes) = size(Q, 1), size(Q, 2) include("ul.jl") include("qr.jl") diff --git a/src/qr.jl b/src/qr.jl index 397c9af..650c8cc 100644 --- a/src/qr.jl +++ b/src/qr.jl @@ -145,6 +145,8 @@ ldiv!(F::QR, B::AbstractVecOrMat) = ArrayLayouts.ldiv!(F, B) ldiv!(F::QR, B::LayoutVector) = ArrayLayouts.ldiv!(F, B) ldiv!(F::QR, B::LayoutMatrix) = ArrayLayouts.ldiv!(F, B) +size(F::QR, dim::Integer) = size(getfield(F, :factors), dim) +size(F::QR) = size(getfield(F, :factors)) """ QRPackedQ <: LinearAlgebra.AbstractQ @@ -163,28 +165,51 @@ end -QRPackedQ(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = QRPackedQ{T,typeof(factors),typeof(τ)}(factors, τ) -function QRPackedQ{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} - QRPackedQ(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ)) +""" + QRPackedQMatrix <: LayoutQMatrix + +The orthogonal/unitary ``Q`` matrix of a QR factorization stored in [`QR`](@ref), +conforming to the `AbstractMatrix` interface. +""" +struct QRPackedQMatrix{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} <: LayoutQMatrix{T} + factors::S + τ::Tau + + function QRPackedQMatrix{T,S,Tau}(factors, τ) where {T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} + require_one_based_indexing(factors) + new{T,S,Tau}(factors, τ) + end +end + + +for QTyp in (:QRPackedQ, :QRPackedQMatrix) + @eval begin + $QTyp(factors::AbstractMatrix{T}, τ::AbstractVector{T}) where {T} = $QTyp{T,typeof(factors),typeof(τ)}(factors, τ) + function $QTyp{T}(factors::AbstractMatrix, τ::AbstractVector) where {T} + $QTyp(convert(AbstractMatrix{T}, factors), convert(AbstractVector{T}, τ)) + end + + $QTyp{T}(Q::$QTyp) where {T} = $QTyp(convert(AbstractMatrix{T}, Q.factors), convert(AbstractVector{T}, Q.τ)) + end end -QRPackedQ{T}(Q::QRPackedQ) where {T} = QRPackedQ(convert(AbstractMatrix{T}, Q.factors), convert(AbstractVector{T}, Q.τ)) QRPackedQ(Q::LinearAlgebra.QRPackedQ) = QRPackedQ(Q.factors, Q.τ) +const QRPackedQTypes{T,S<:AbstractMatrix{T},Tau<:AbstractVector{T}} = Union{QRPackedQ{T,S,Tau}, QRPackedQMatrix{T,S,Tau}} + +Matrix{T}(Q::QRPackedQTypes{S}) where {T,S} = convert(Matrix{T}, lmul!(Q, Matrix{S}(I, size(Q, 1), min(size(Q.factors)...)))) +Matrix(Q::QRPackedQTypes{S}) where {S} = Matrix{S}(Q) + AbstractQ{T}(Q::QRPackedQ{T}) where {T} = Q AbstractQ{T}(Q::QRPackedQ) where {T} = QRPackedQ{T}(Q) +AbstractQ{T}(Q::QRPackedQMatrix{T}) where {T} = QRPackedQ(Q.factors, Q.τ) +AbstractQ{T}(Q::QRPackedQMatrix) where {T} = QRPackedQ{T}(Q.factors, Q.τ) +AbstractQ(Q::QRPackedQMatrix) = QRPackedQ(Q.factors, Q.τ) convert(::Type{AbstractQ{T}}, Q::QRPackedQ) where {T} = QRPackedQ{T}(Q) - -Matrix{T}(Q::QRPackedQ{S}) where {T,S} = - convert(Matrix{T}, lmul!(Q, Matrix{S}(I, size(Q, 1), min(size(Q.factors)...)))) -Matrix(Q::QRPackedQ{S}) where {S} = Matrix{S}(Q) - AbstractMatrix{T}(Q::QRPackedQ) where {T} = Matrix{T}(Q) -size(F::QR, dim::Integer) = size(getfield(F, :factors), dim) -size(F::QR) = size(getfield(F, :factors)) -MemoryLayout(::Type{<:QRPackedQ{<:Any,S,T}}) where {S,T} = +MemoryLayout(::Type{<:QRPackedQTypes{<:Any,S,T}}) where {S,T} = QRPackedQLayout{typeof(MemoryLayout(S)),typeof(MemoryLayout(T))}() MemoryLayout(::Type{<:QR{<:Any,S,T}}) where {S,T} = @@ -227,3 +252,9 @@ function (\)(A::QR{T}, BIn::VecOrMat{Complex{T}}) where T<:BlasReal XX = reshape(collect(reinterpret(Complex{T}, copy(transpose(reshape(X, div(length(X), 2), 2))))), _ret_size(A, BIn)) return _cut_B(XX, 1:n) end + + + +# support lazy broadcasting + +BroadcastStyle(::Type{<:QRPackedQMatrix{<:Any,F}}) where {F} = BroadcastStyle(F) # TODO: broken for banded? \ No newline at end of file