Adjust to !(AbstractQ <: AbstractMatrix) (#122)

* Adjust to `!(AbstractQ <: AbstractMatrix)`

* be more restrictive

Author: dkarrasch

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.6.13"
+version = "0.6.14"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/InfiniteLinearAlgebra.jl

@@ -5,7 +5,8 @@ using BlockArrays, BlockBandedMatrices, BandedMatrices, LazyArrays, LazyBandedMa
 
 import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, size, axes, length,
             AbstractMatrix, AbstractArray, Matrix, Array, Vector, AbstractVector, Slice,
-            show, getproperty, copy, copyto!, map, require_one_based_indexing, similar, inv
+            show, getproperty, copy, copyto!, map, require_one_based_indexing, similar, inv,
+            oneto, unitrange
 import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
 import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
@@ -17,7 +18,7 @@ import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
 import InfiniteArrays: OneToInf, InfUnitRange, Infinity, PosInfinity, InfiniteCardinal, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
 import LinearAlgebra: matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose, AdjOrTrans
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, LazyArrayStyle,
-                    resizedata!, MemoryLayout, most,
+                    resizedata!, MemoryLayout,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
                     applylayout, ApplyLayout, PaddedLayout, CachedLayout, cacheddata, zero!, MulAddStyle,
                     LazyArray, LazyMatrix, LazyVector, paddeddata, arguments
@@ -40,17 +41,13 @@ import DSP: conv
 import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
 
-if VERSION < v"1.6-"
-    oneto(n) = Base.OneTo(n)
-else
-    import Base: oneto, unitrange
-end
-
 if VERSION ≥ v"1.7-"
-    LinearAlgebra._cut_B(x::AbstractVector, r::InfUnitRange) = x
-    LinearAlgebra._cut_B(X::AbstractMatrix, r::InfUnitRange) = X
+    LinearAlgebra._cut_B(x::AbstractVector, ::InfUnitRange) = x
+    LinearAlgebra._cut_B(X::AbstractMatrix, ::InfUnitRange) = X
 end
 
+const AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint
+
 # BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
 
 function ArrayLayouts._power_by_squaring(_, ::NTuple{2,InfiniteCardinal{0}}, A::AbstractMatrix{T}, p::Integer) where T

src/banded/hessenbergq.jl

@@ -1,33 +1,39 @@
 isorthogonal(::AbstractQ) = true
 isorthogonal(q) = q'q ≈ I
 
+convert_eltype(Q::AbstractMatrix, ::Type{T}) where {T} = convert(AbstractMatrix{T}, Q)
+if !(AbstractQ <: AbstractMatrix)
+    convert_eltype(Q::AbstractQ, ::Type{T}) where {T} = convert(AbstractQ{T}, Q)
+end
+
 """
     QLHessenberg(factors, q)
 
 represents a Hessenberg QL factorization where factors contains L in its
 lower triangular components and q is a vector of 2x2 orthogonal transformations
 whose product gives Q.
 """
-struct QLHessenberg{T,S<:AbstractMatrix{T},QT<:AbstractVector{<:AbstractMatrix{T}}} <: Factorization{T}
+struct QLHessenberg{T,S<:AbstractMatrix{T},QT<:AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}} <: Factorization{T}
     factors::S
     q::QT
 
-    function QLHessenberg{T,S,QT}(factors, q) where {T,S<:AbstractMatrix{T},QT<:AbstractVector{<:AbstractMatrix{T}}}
+    function QLHessenberg{T,S,QT}(factors, q) where {T,S<:AbstractMatrix{T},QT<:AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}}
         require_one_based_indexing(factors)
         new{T,S,QT}(factors, q)
     end
 end
 
-QLHessenberg(factors::AbstractMatrix{T}, q::AbstractVector{<:AbstractMatrix{T}}) where {T} = QLHessenberg{T,typeof(factors),typeof(q)}(factors, q)
+QLHessenberg(factors::AbstractMatrix{T}, q::AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}) where {T} =
+    QLHessenberg{T,typeof(factors),typeof(q)}(factors, q)
 QLHessenberg{T}(factors::AbstractMatrix, q::AbstractVector) where {T} =
-    QLHessenberg(convert(AbstractMatrix{T}, factors), convert.(AbstractMatrix{T}, q))
+    QLHessenberg(convert(AbstractMatrix{T}, factors), convert_eltype.(q, T))
 
 # iteration for destructuring into components
 Base.iterate(S::QLHessenberg) = (S.Q, Val(:L))
 Base.iterate(S::QLHessenberg, ::Val{:L}) = (S.L, Val(:done))
 Base.iterate(S::QLHessenberg, ::Val{:done}) = nothing
 
-QLHessenberg{T}(A::QLHessenberg) where {T} = QLHessenberg(convert(AbstractMatrix{T}, A.factors), convert.(AbstractMatrix{T}, A.q))
+QLHessenberg{T}(A::QLHessenberg) where {T} = QLHessenberg(convert(AbstractMatrix{T}, A.factors), convert_eltype.(A.q, T))
 Factorization{T}(A::QLHessenberg{T}) where {T} = A
 Factorization{T}(A::QLHessenberg) where {T} = QLHessenberg{T}(A)
 AbstractMatrix(F::QLHessenberg) = F.Q * F.L
@@ -69,16 +75,16 @@ abstract type AbstractHessenbergQ{T} <: LayoutQ{T} end
 
 for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
     @eval begin
-        struct $Typ{T, QT<:AbstractVector{<:AbstractMatrix{T}}} <: AbstractHessenbergQ{T}
+        struct $Typ{T, QT<:AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}} <: AbstractHessenbergQ{T}
             q::QT
-            function $Typ{T,QT}(q::QT) where {T,QT<:AbstractVector{<:AbstractMatrix{T}}}
+            function $Typ{T,QT}(q::QT) where {T,QT<:AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}}
                 all(isorthogonal.(q)) || throw(ArgumentError("input must be orthogonal"))
                 all(size.(q) .== Ref((2,2))) ||  throw(ArgumentError("input must be 2x2"))
                 new{T,QT}(q)
             end
         end
 
-        $Typ(q::AbstractVector{<:AbstractMatrix{T}}) where T =
+        $Typ(q::AbstractVector{<:Union{AbstractMatrix{T},AbstractQ{T}}}) where {T} =
             $Typ{T,typeof(q)}(q)
     end
 end
@@ -141,13 +147,12 @@ end
 # getindex
 ####
 
-function getindex(Q::UpperHessenbergQ, i::Int, j::Int)
-    y = zeros(eltype(Q), size(Q, 2))
-    y[j] = 1
-    lmul!(Q, y)[i]
-end
+getindex(Q::UpperHessenbergQ, I::AbstractVector{Int}, J::AbstractVector{Int}) =
+    hcat((Q[:,j][I] for j in J)...)
 
 getindex(Q::LowerHessenbergQ, i::Int, j::Int) = (Q')[j,i]'
+getindex(Q::LowerHessenbergQ, I::AbstractVector{Int}, J::AbstractVector{Int}) =
+    [Q[i,j] for i in I, j in J]
 
 
 ###

src/banded/infbanded.jl

@@ -15,7 +15,7 @@ getindex(B::InfBandCartesianIndices, k::Int) = B.b ≥ 0 ? CartesianIndex(k, k+B
 
 Base.checkindex(::Type{Bool}, ::NTuple{2,OneToInf{Int}}, ::InfBandCartesianIndices) = true
 BandedMatrices.band_to_indices(_, ::NTuple{2,OneToInf{Int}}, b) = (InfBandCartesianIndices(b),)
-Base.BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,Tuple{InfBandCartesianIndices}}}) = LazyArrayStyle{1}()
+BroadcastStyle(::Type{<:SubArray{<:Any,1,<:Any,Tuple{InfBandCartesianIndices}}}) = LazyArrayStyle{1}()
 
 _inf_banded_sub_materialize(_, V) = V
 function _inf_banded_sub_materialize(::BandedColumns, V)
@@ -487,7 +487,7 @@ for Typ in (:(LinearAlgebra.Tridiagonal{<:Any,<:InfFill}),
             :(LazyBandedMatrices.SymTridiagonal{<:Any,<:InfFill,<:InfFill}))
     @eval begin
         MemoryLayout(::Type{<:$Typ}) = TridiagonalToeplitzLayout()
-        Base.BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}()
+        BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}()
     end
 end
 
@@ -497,7 +497,7 @@ for Typ in (:(LinearAlgebra.Bidiagonal{<:Any,<:InfFill}),
             :(LazyBandedMatrices.Bidiagonal{<:Any,<:InfFill,<:InfFill}))
     @eval begin
         MemoryLayout(::Type{<:$Typ}) = BidiagonalToeplitzLayout()
-        Base.BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}()
+        BroadcastStyle(::Type{<:$Typ}) = LazyArrayStyle{2}()
     end
 end
 

src/banded/infqltoeplitz.jl

@@ -88,10 +88,18 @@ struct ProductQ{T,QQ<:Tuple} <: LayoutQ{T}
     Qs::QQ
 end
 
-ArrayLayouts.@layoutmatrix ProductQ
-ArrayLayouts.@_layoutlmul ProductQ
+if VERSION < v"1.8-"
+    ArrayLayouts.@layoutmatrix ProductQ
+    ArrayLayouts.@_layoutlmul ProductQ
+else # ArrayLayouts.@layoutmatrix ProductQ without ArrayLayouts.@layoutgetindex
+    ArrayLayouts.@layoutldiv ProductQ
+    ArrayLayouts.@layoutmul ProductQ
+    ArrayLayouts.@layoutlmul ProductQ
+    ArrayLayouts.@layoutfactorizations ProductQ
+    ArrayLayouts.@_layoutlmul ProductQ
+end
 
-ProductQ(Qs::AbstractMatrix...) = ProductQ{mapreduce(eltype, promote_type, Qs),typeof(Qs)}(Qs)
+ProductQ(Qs::Union{AbstractMatrix,AbstractQ}...) = ProductQ{mapreduce(eltype, promote_type, Qs),typeof(Qs)}(Qs)
 
 adjoint(Q::ProductQ) = ProductQ(reverse(map(adjoint, Q.Qs))...)
 
@@ -105,14 +113,21 @@ function lmul!(Q::ProductQ, v::AbstractVecOrMat)
     v
 end
 
-# Avoid ambiguities
-getindex(Q::ProductQ, i::Int, j::Int) = Q[:, j][i]
+if VERSION < v"1.8-"
+    # Avoid ambiguities
+    getindex(Q::ProductQ, i::Int, j::Int) = Q[:, j][i]
 
-function getindex(Q::ProductQ, ::Colon, j::Int)
-    y = zeros(eltype(Q), size(Q, 2))
-    y[j] = 1
-    lmul!(Q, y)
+    function getindex(Q::ProductQ, ::Colon, j::Int)
+        y = zeros(eltype(Q), size(Q, 2))
+        y[j] = 1
+        lmul!(Q, y)
+    end
 end
+if VERSION >= v"1.10-"
+    getindex(Q::ProductQ, I::AbstractVector{Int}, J::AbstractVector{Int}) =
+        hcat((Q[:,j][I] for j in J)...)
+end
+
 getindex(Q::ProductQ{<:Any,<:Tuple{Vararg{LowerHessenbergQ}}}, i::Int, j::Int) = (Q')[j, i]'
 
 function _productq_mul(A::ProductQ{T}, x::AbstractVector{S}) where {T,S}

src/blockbanded/blockbanded.jl

@@ -1,4 +1,4 @@
-const OneToInfCumsum = InfiniteArrays.RangeCumsum{Int,OneToInf{Int}}
+const OneToInfCumsum = RangeCumsum{Int,OneToInf{Int}}
 
 BlockArrays.sortedunion(::AbstractVector{<:PosInfinity}, ::AbstractVector{<:PosInfinity}) = [∞]
 function BlockArrays.sortedunion(::AbstractVector{<:PosInfinity}, b)

src/infql.jl

@@ -145,7 +145,7 @@ end
 function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:BandedColumns},<:PaddedLayout})
     adjA,B = M.A,M.B
     require_one_based_indexing(B)
-    A = adjA.parent
+    A = parent(adjA)
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
     if mA != mB
@@ -179,9 +179,9 @@ function _lmul_cache(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
 end
 
 (*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::AbstractVector) where {T} = _lmul_cache(A, x)
-(*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector) where {T} = _lmul_cache(A, x)
+(*)(A::AdjointQtype{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::AbstractVector) where {T} = _lmul_cache(A, x)
 (*)(A::QLPackedQ{T,<:InfBandedMatrix}, x::LazyVector) where {T} = _lmul_cache(A, x)
-(*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::LazyVector) where {T} = _lmul_cache(A, x)
+(*)(A::AdjointQtype{T,<:QLPackedQ{T,<:InfBandedMatrix}}, x::LazyVector) where {T} = _lmul_cache(A, x)
 
 
 function blocktailiterate(c,a,b, d=c, e=a)
@@ -230,9 +230,9 @@ ql(A::Adjoint{T,BlockTriPertToeplitz{T}}) where T = ql(BlockTridiagonal(A))
 
 const InfBlockBandedMatrix{T} = BlockSkylineMatrix{T,<:Vcat{T,1,<:Tuple{Vector{T},<:BlockArray{T,1,<:Fill{<:Any,1,Tuple{OneToInf{Int64}}}}}}}
 
-function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:InfBlockBandedMatrix}}, B::AbstractVector)
+function lmul!(adjA::AdjointQtype{<:Any,<:QLPackedQ{<:Any,<:InfBlockBandedMatrix}}, B::AbstractVector)
     require_one_based_indexing(B)
-    A = adjA.parent
+    A = parent(adjA)
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
     if mA != mB
@@ -270,7 +270,7 @@ function (*)(A::QLPackedQ{T,<:InfBlockBandedMatrix}, x::AbstractVector{S}) where
     lmul!(A, cache(convert(AbstractVector{TS},x)))
 end
 
-function (*)(A::Adjoint{T,<:QLPackedQ{T,<:InfBlockBandedMatrix}}, x::AbstractVector{S}) where {T,S}
+function (*)(A::AdjointQtype{T,<:QLPackedQ{T,<:InfBlockBandedMatrix}}, x::AbstractVector{S}) where {T,S}
     TS = promote_op(matprod, T, S)
     lmul!(A, cache(convert(AbstractVector{TS},x)))
 end
@@ -308,7 +308,7 @@ _data_tail(::CachedLayout, a) = cacheddata(a), getindex_value(a.array)
 function _data_tail(::ApplyLayout{typeof(vcat)}, a)
     args = arguments(vcat, a)
     dat,tl = _data_tail(last(args))
-    vcat(most(args)..., dat), tl
+    vcat(Base.front(args)..., dat), tl
 end
 _data_tail(a) = _data_tail(MemoryLayout(a), a)
 

src/infqr.jl

@@ -158,6 +158,13 @@ _factorize(::AbstractBandedLayout, ::NTuple{2,OneToInf{Int}}, A) = qr(A)
 partialqr!(F::QR, n) = partialqr!(F.factors, n)
 partialqr!(F::AdaptiveQRFactors, n) = partialqr!(F.data, n)
 
+#########
+# getindex
+#########
+
+getindex(Q::QRPackedQ{<:Any,<:AdaptiveQRFactors,<:AdaptiveQRTau}, I::AbstractVector{Int}, J::AbstractVector{Int64}) =
+    hcat((Q[:,j][I] for j in J)...)
+
 #########
 # lmul!
 #########
@@ -210,7 +217,7 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout},<:Pad
     COLGROWTH = 1000 # rate to grow columns
 
     require_one_based_indexing(B)
-    A = adjA.parent
+    A = parent(adjA)
     mA, nA = size(A.factors)
     mB = length(B)
     if mA != mB
@@ -274,7 +281,7 @@ end
 
 function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:AbstractBlockBandedLayout}},<:PaddedLayout}; tolerance=1E-30)
     adjA,B_in = M.A,M.B
-    A = adjA.parent
+    A = parent(adjA)
     T = eltype(M)
     COLGROWTH = 300 # rate to grow columns
     ax1,ax2 = axes(A.factors.data.data)
@@ -309,15 +316,15 @@ function materialize!(M::MatLmulVec{<:AdjQRPackedQLayout{<:AdaptiveLayout{<:Abst
 end
 
 
-function _lmul_copymutable(A::AbstractMatrix{T}, x::AbstractVector{S}; kwds...) where {T,S}
+function _lmul_copymutable(A::Union{AbstractMatrix{T},AbstractQ{T}}, x::AbstractVector{S}; kwds...) where {T,S}
     TS = promote_op(matprod, T, S)
     lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)); kwds...)
 end
 
 (*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::AbstractVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
-(*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
+(*)(A::AdjointQtype{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
 (*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::LayoutVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
-(*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::LayoutVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
+(*)(A::AdjointQtype{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::LayoutVector; kwds...) where {T} = _lmul_copymutable(A, x; kwds...)
 
 function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}})
     n = B.datasize[1]

test/runtests.jl

@@ -150,7 +150,7 @@ include("test_infbanded.jl")
             N = 1000
             v = view(n, Block.(Base.OneTo(N)))
             @test view(v, Block(2)) ≡ Fill(2, 2)
-            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
+            @test axes(v) isa Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64,Base.OneTo{Int64}}}}
             @test @allocated(axes(v)) ≤ 40
 
             dest = PseudoBlockArray{Float64}(undef, axes(v))
@@ -159,7 +159,7 @@ include("test_infbanded.jl")
 
             v = view(k, Block.(Base.OneTo(N)))
             @test view(v, Block(2)) ≡ Base.OneTo(2)
-            @test axes(v) isa Tuple{BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,Base.OneTo{Int64}}}}
+            @test axes(v) isa Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64,Base.OneTo{Int64}}}}
             @test @allocated(axes(v)) ≤ 40
             @test copyto!(dest, v) == v
 
@@ -173,7 +173,7 @@ include("test_infbanded.jl")
             @test v[Block(1)] == 1:2
             @test v[Block(1)] ≡ k[Block(2)] ≡ Base.OneTo(2)
 
-            @test axes(n, 1) isa BlockedUnitRange{InfiniteArrays.RangeCumsum{Int64,OneToInf{Int64}}}
+            @test axes(n, 1) isa BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64,OneToInf{Int64}}}
         end
 
         @testset "BlockHcat copyto!" begin