Dl/blockbandedqr (#29)

* Block banded QR

* Block adaptive QR

* Update test_infqr.jl

* Block-+ adaptive QR

* AbstractQ

* Update Project.toml

* Update runtests.jl

* Tests pass

* Update Project.toml

* Update Project.toml

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,16 +16,16 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.2.4"
-BandedMatrices = "0.15.3"
-BlockArrays = "0.12"
-BlockBandedMatrices = "0.8"
-FillArrays = "0.8.6"
+ArrayLayouts = "0.3"
+BandedMatrices = "0.15.7"
+BlockArrays = "0.12.6"
+BlockBandedMatrices = "0.8.4"
+FillArrays = "0.8.8"
 InfiniteArrays = "0.7"
-LazyArrays = "0.16.5"
-LazyBandedMatrices = "0.2.6"
-MatrixFactorizations = "0.3.1"
-SemiseparableMatrices = "0.0.1"
+LazyArrays = "0.16.8"
+LazyBandedMatrices = "0.2.8"
+MatrixFactorizations = "0.4.1"
+SemiseparableMatrices = "0.1"
 julia = "1.3"
 
 [extras]

src/InfiniteLinearAlgebra.jl

@@ -7,7 +7,8 @@ import Base: +, -, *, /, \, ^, OneTo, getindex, promote_op, _unsafe_getindex, pr
             show, getproperty, copy, map, require_one_based_indexing
 import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
-import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout, _qr
+import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, sublayout, _qr, 
+                        MatLmulVec, MatLmulMat, AbstractQLayout, materialize!
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast
 import FillArrays: AbstractFill, getindex_value      
@@ -16,18 +17,20 @@ import LinearAlgebra: lmul!,  ldiv!, matprod, qr, AbstractTriangular, AbstractQ,
 import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
                     resizedata!, MemoryLayout, mulapplystyle, LmulStyle, RmulStyle,
                     factorize, sub_materialize, LazyLayout, LazyArrayStyle, layout_getindex,
-                    @lazymul, applylayout, ApplyLayout, PaddedLayout, materialize!, zero!, MulAddStyle,
-                    LazyArray, LazyMatrix, LazyVector
-import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ
+                    @lazymul, applylayout, ApplyLayout, PaddedLayout, zero!, MulAddStyle,
+                    LazyArray, LazyMatrix, LazyVector, paddeddata
+import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, getL, getR, reflector!, reflectorApply!, QL, QR, QRPackedQ,
+                            QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
 
-import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange
+import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport
 
 import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
 
-import LazyBandedMatrices: MulBandedLayout, BroadcastBandedLayout
+import LazyBandedMatrices: MulBandedLayout, BroadcastBandedLayout, _krontrav_axes
 
 import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
-        BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal
+        BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal,
+        AbstractBlockBandedLayout, _blockbanded_qr!, BlockBandedLayout
 
 import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
 
@@ -78,5 +81,13 @@ broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{<:Fill,1,Ty
 
 BlockArrays._length(::BlockedUnitRange, ::OneToInf) = ∞
 BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
+
+###
+# KronTrav
+###
+
+_krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N = 
+     @. blockedrange(OneTo(length(A)))
+
 
 end # module

src/banded/hessenbergq.jl

@@ -1,7 +1,13 @@
 isorthogonal(::AbstractQ) = true
 isorthogonal(q) = q'q ≈ I
 
+"""
+    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}
     factors::S
     q::QT
@@ -43,8 +49,8 @@ end
 end
 @inline getQ(F::QLHessenberg, _) = LowerHessenbergQ(F.q)
 
-getL(F::QLHessenberg) = getL(F, axes(F.factors))
-getQ(F::QLHessenberg) = getQ(F, axes(F.factors))
+getL(F::QLHessenberg) = getL(F, size(F.factors))
+getQ(F::QLHessenberg) = getQ(F, size(F.factors))
 
 function getproperty(F::QLHessenberg, d::Symbol)
     if d == :L
@@ -59,7 +65,7 @@ end
 Base.propertynames(F::QLHessenberg, private::Bool=false) =
     (:L, :Q, (private ? fieldnames(typeof(F)) : ())...)
 
-abstract type AbstractHessenbergQ{T} <: AbstractQ{T} end
+abstract type AbstractHessenbergQ{T} <: LayoutQ{T} end
 
 for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
     @eval begin
@@ -78,6 +84,7 @@ for Typ in (:UpperHessenbergQ, :LowerHessenbergQ)
 end
 
 size(Q::AbstractHessenbergQ, k::Integer) = length(Q.q)+1
+axes(Q::AbstractHessenbergQ, k::Integer) = Base.OneTo(length(Q.q)+1)
 size(F::QLHessenberg, dim::Integer) = size(getfield(F, :factors), dim)
 size(F::QLHessenberg) = size(getfield(F, :factors))
 
@@ -90,17 +97,14 @@ adjoint(Q::LowerHessenbergQ) = UpperHessenbergQ(adjoint.(Q.q))
 check_mul_axes(A::AbstractHessenbergQ, B, C...) =
     axes(A,2) == axes(B,1) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and first axis of B, $(axes(B,1)) must match"))
 
-@lazymul AbstractHessenbergQ
+struct HessenbergQLayout{UPLO} <: AbstractQLayout end
 
-# ambiguities
-for Arr in (:AbstractBandedMatrix, :StridedVector, :StridedMatrix, :(BandedMatrix{<:Any,<:Any,<:OneToInf}))
-    @eval begin
-        *(A::AbstractHessenbergQ, B::$Arr) = apply(*, A, B)
-        *(A::$Arr, B::AbstractHessenbergQ) = apply(*, A, B)
-    end
-end
+MemoryLayout(::Type{<:UpperHessenbergQ}) = HessenbergQLayout{'U'}()
+MemoryLayout(::Type{<:LowerHessenbergQ}) = HessenbergQLayout{'L'}()
 
-function lmul!(Q::LowerHessenbergQ{T}, x::AbstractVector) where T
+function materialize!(L::MatLmulVec{<:HessenbergQLayout{'L'}})
+    Q, x = L.A,L.B
+    T = eltype(Q)
     t = Array{T}(undef, 2)
     nz = nzzeros(x,1)
     for n = 1:length(Q.q)
@@ -112,7 +116,9 @@ function lmul!(Q::LowerHessenbergQ{T}, x::AbstractVector) where T
     x
 end
 
-function lmul!(Q::UpperHessenbergQ{T}, x::AbstractVector) where T
+function materialize!(L::MatLmulVec{<:HessenbergQLayout{'U'}})
+    Q, x = L.A,L.B
+    T = eltype(Q)
     t = Array{T}(undef, 2)
     for n = min(length(Q.q),nzzeros(x,1)):-1:1
         v = view(x, n:n+1)
@@ -122,9 +128,10 @@ function lmul!(Q::UpperHessenbergQ{T}, x::AbstractVector) where T
     x
 end
 
-function lmul!(Q::AbstractHessenbergQ, X::AbstractMatrix)
+function materialize!(L::MatLmulMat{<:HessenbergQLayout})
+    Q,X = L.A,L.B
     for j in axes(X,2)
-        lmul!(Q, view(X,:,j))
+        ArrayLayouts.lmul!(Q, view(X,:,j))
     end
     X
 end

src/banded/infbanded.jl

@@ -374,5 +374,5 @@ mulapplystyle(::BandedColumns{FillLayout}, ::PertToeplitzLayout) = LazyArrayAppl
 mulapplystyle(::PertToeplitzLayout, ::BandedColumns{FillLayout}) = LazyArrayApplyStyle()
 mulapplystyle(::BandedColumns{FillLayout}, ::BandedToeplitzLayout) = LazyArrayApplyStyle()
 mulapplystyle(::BandedToeplitzLayout, ::BandedColumns{FillLayout}) = LazyArrayApplyStyle()
-mulapplystyle(::QLayout, ::BandedToeplitzLayout) = LazyArrayApplyStyle()
-mulapplystyle(::QLayout, ::PertToeplitzLayout) = LazyArrayApplyStyle()
+mulapplystyle(::AbstractQLayout, ::BandedToeplitzLayout) = LazyArrayApplyStyle()
+mulapplystyle(::AbstractQLayout, ::PertToeplitzLayout) = LazyArrayApplyStyle()

src/blockbanded/blockbanded.jl

@@ -16,7 +16,7 @@ blocklayout(::LazyLayout) = BlockLayout{LazyLayout}()
 # for LazyLay in (:(BlockLayout{LazyLayout}), :(TriangularLayout{UPLO,UNIT,BlockLayout{LazyLayout}} where {UPLO,UNIT}))
 #     @eval begin
 #         combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()
-#         mulapplystyle(::QLayout, ::$LazyLay) = LazyArrayApplyStyle()
+#         mulapplystyle(::AbstractQLayout, ::$LazyLay) = LazyArrayApplyStyle()
 #     end
 # end
 

src/infql.jl

@@ -97,10 +97,10 @@ function ql_hessenberg!(B::InfBandedMatrix{TT}; kwds...) where TT
 end
 
 getindex(Q::QLPackedQ{T,<:InfBandedMatrix{T}}, i::Integer, j::Integer) where T =
-    (Q'*Vcat(Zeros{T}(i-1), one(T), Zeros{T}(∞)))[j]'
+    (Q'*[Zeros{T}(i-1); one(T); Zeros{T}(∞)])[j]'
 
-getL(Q::QL, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = LowerTriangular(Q.factors)
-getL(Q::QLHessenberg, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = LowerTriangular(Q.factors)
+getL(Q::QL, ::NTuple{2,Infinity}) where T = LowerTriangular(Q.factors)
+getL(Q::QLHessenberg, ::NTuple{2,Infinity}) where T = LowerTriangular(Q.factors)
 
 # number of structural non-zeros in axis k
 nzzeros(A::AbstractArray, k) = size(A,k)
@@ -112,7 +112,8 @@ function nzzeros(B::AbstractMatrix, k)
     k == 1 ? size(B,2) + l : size(B,1) + u
 end
 
-function lmul!(A::QLPackedQ{<:Any,<:InfBandedMatrix}, B::AbstractVecOrMat)
+function materialize!(M::Lmul{<:QLPackedQLayout{<:BandedColumns},<:PaddedLayout})
+    A,B = M.A,M.B
     require_one_based_indexing(B)
     mA, nA = size(A.factors)
     mB, nB = size(B,1), size(B,2)
@@ -145,7 +146,8 @@ function lmul!(A::QLPackedQ{<:Any,<:InfBandedMatrix}, B::AbstractVecOrMat)
     B
 end
 
-function lmul!(adjA::Adjoint{<:Any,<:QLPackedQ{<:Any,<:InfBandedMatrix}}, B::AbstractVecOrMat)
+function materialize!(M::Lmul{<:AdjQLPackedQLayout{<:BandedColumns},<:PaddedLayout})
+    adjA,B = M.A,M.B
     require_one_based_indexing(B)
     A = adjA.parent
     mA, nA = size(A.factors)