Generalize QR (#128)

* Support Julia v1.3

* Skip Float16 det test, add tests for #118

* v0.10.1, add missing similar

* Start to Generalize QR to support subviews

* FillArrays 0.7

* updates for latest LazyArrays

* Add special copy for BandedStyle broadcasting, use lazy_getindex for more types

* Fix materialize!, BandedMatrices 0.11

Project.toml

@@ -1,6 +1,6 @@
 name = "BandedMatrices"
 uuid = "aae01518-5342-5314-be14-df237901396f"
-version = "0.10.1"
+version = "0.11"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
@@ -11,8 +11,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-FillArrays = "0.6.4"
-LazyArrays = "0.10"
+FillArrays = "0.6.4, 0.7"
+LazyArrays = "0.11"
 MatrixFactorizations = "0.1"
 julia = "1"
 

src/BandedMatrices.jl

@@ -31,19 +31,22 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Br
 
 import LazyArrays: MemoryLayout, @lazymul, @lazylmul, @lazyldiv,
                     AbstractStridedLayout, AbstractColumnMajor, AbstractRowMajor,
+                    LazyLayout, UnknownLayout,
                     transposelayout, triangulardata,
                     ConjLayout, conjlayout, SymmetricLayout, symmetriclayout, symmetricdata,
                     triangularlayout, MatLdivVec, TriangularLayout,
-                    AbstractBandedLayout, DiagonalLayout,
+                    AbstractBandedLayout, DiagonalLayout, AbstractFillLayout,
                     HermitianLayout, hermitianlayout, hermitiandata,
                     MulAdd, materialize!, BlasMatMulMatAdd, BlasMatMulVecAdd, BlasMatLmulVec, BlasMatLdivVec,
-                    VcatLayout, ZerosLayout,
+                    ZerosLayout,
                     AbstractColumnMajor, MulLayout, colsupport, rowsupport,
                     DenseColumnMajor, DenseRowMajor, ApplyArrayBroadcastStyle,
-                    mulapplystyle, AbstractMulAddStyle, symmetricuplo, MatMulMatAdd, MatMulVecAdd,
-                    _fill_lmul!, applybroadcaststyle, subarraylayout, sub_materialize, lazy_getindex
+                    symmetricuplo, MatMulMatAdd, MatMulVecAdd,
+                    _fill_lmul!, applybroadcaststyle, subarraylayout, sub_materialize, lazy_getindex,
+                    resizedata!, CachedMatrix, ApplyLayout,
+                    combine_mul_styles, LazyArrayApplyStyle, QLayout, mulapplystyle
 
-import FillArrays: AbstractFill
+import FillArrays: AbstractFill, getindex_value
 
 export BandedMatrix,
        bandrange,

src/banded/BandedLU.jl

@@ -91,7 +91,8 @@ end
 lu!(A::AbstractBandedMatrix, pivot::Union{Val{false}, Val{true}} = Val(true); check::Bool = true) =
     banded_lufact!(A, pivot; check = check)
 
-function lu(A::Union{AbstractBandedMatrix{T}, AbstractBandedMatrix{Complex{T}}},
+function lu(A::Union{AbstractBandedMatrix{T}, AbstractBandedMatrix{Complex{T}}, Adjoint{T,<:AbstractBandedMatrix{T}}, Adjoint{Complex{T},<:AbstractBandedMatrix{Complex{T}}},
+                    Transpose{T,<:AbstractBandedMatrix{T}}, Transpose{Complex{T},<:AbstractBandedMatrix{Complex{T}}}},
     pivot::Union{Val{false}, Val{true}} = Val(true);
     check::Bool = true) where {T<:Real}
     l,u = bandwidths(A)

src/banded/BandedMatrix.jl

@@ -29,7 +29,18 @@ end
 
 _BandedMatrix(data::AbstractMatrix, m::Integer, l, u) = _BandedMatrix(data, Base.OneTo(m), l, u)
 
-MemoryLayout(A::Type{BandedMatrix{T,Cont,Axes}}) where {T,Cont,Axes} = BandedColumns{typeof(MemoryLayout(Cont))}()
+bandedcolumns(_) = BandedColumns{UnknownLayout}()
+bandedcolumns(::ML) where ML<:AbstractStridedLayout = BandedColumns{ML}()
+bandedcolumns(::ML) where ML<:AbstractFillLayout = BandedColumns{ML}()
+bandedcolumns(::ML) where ML<:LazyLayout = BandedColumns{ML}()
+bandedcolumns(::ML) where ML<:ApplyLayout = BandedColumns{LazyLayout}()
+
+combine_mul_styles(::BandedColumns{LazyLayout}) = LazyArrayApplyStyle()
+combine_mul_styles(::BandedRows{LazyLayout}) = LazyArrayApplyStyle()
+
+mulapplystyle(::QLayout, ::BandedColumns{LazyLayout}) = LazyArrayApplyStyle()
+
+MemoryLayout(A::Type{BandedMatrix{T,Cont,Axes}}) where {T,Cont,Axes} = bandedcolumns(MemoryLayout(Cont))
 
 
 ## Constructors
@@ -257,12 +268,11 @@ function similar(bm::BandedMatrix, T::Type=eltype(bm),
     _BandedMatrix(data, n, l, u)
 end
 
-similar(bm::AbstractBandedMatrix, n::Integer, m::Integer) = similar(bm, eltype(bm), m, n)
+similar(bm::AbstractBandedMatrix, n::Integer, m::Integer) = similar(bm, eltype(bm), n, m)
 similar(bm::AbstractBandedMatrix, n::Integer, m::Integer, l::Integer, u::Integer) =
     similar(bm, eltype(bm), m, n, l, u)
 similar(bm::AbstractBandedMatrix, nm::Tuple{<:Integer,<:Integer}) = similar(bm, nm...)
 
-
 ## Abstract Array Interface
 
 axes(A::BandedMatrix) = (A.raxis, axes(A.data,2))
@@ -664,9 +674,6 @@ end
 
 ## BandedSubBandedMatrix routines
 
-
-# getindex(A::AbstractBandedMatrix, I...) = _materialize(view(A, I...))
-
 # gives the band which is diagonal for the parent
 bandshift(a::AbstractRange, b::AbstractRange) = first(a)-first(b)
 bandshift(::Slice{OneTo{Int}}, b::AbstractRange) = 1-first(b)
@@ -682,7 +689,7 @@ const BandedSubBandedMatrix{T, C, R, I1<:AbstractUnitRange, I2<:AbstractUnitRang
     SubArray{T,2,BandedMatrix{T, C, R},Tuple{I1,I2},t}
 
 isbanded(::BandedSubBandedMatrix) = true
-MemoryLayout(::Type{BandedSubBandedMatrix{T,C,R,I1,I2,t}}) where {T,C,R,I1,I2,t} = 
+MemoryLayout(::Type{BandedSubBandedMatrix{T,C,R,I1,I2,t}}) where {T,C,R,I1<:AbstractUnitRange,I2<:AbstractUnitRange,t} = 
     BandedColumns{typeof(MemoryLayout(SubArray{T,2,C,Tuple{Slice{OneTo{Int}},I2},t}))}()
 BroadcastStyle(::Type{<:BandedSubBandedMatrix}) = BandedStyle()
 

src/banded/bandedqr.jl

@@ -5,12 +5,11 @@ qr(A::Tridiagonal{T}) where T = qr!(BandedMatrix{float(T)}(A, (1,2)))
 qr!(A::BandedMatrix) = banded_qr!(A)
 
 
-function banded_qr!(R::BandedMatrix{T}) where T
+function banded_qr!(R::AbstractMatrix{T}, τ) where T
     D = bandeddata(R)
     l,u = bandwidths(R)
     ν = l+u+1
     m,n=size(R)
-    τ = zeros(T, min(m,n))
 
     for k = 1:min(m - 1 + !(T<:Real), n)
         x = view(D,u+1:min(ν,m-k+u+1), k)
@@ -23,6 +22,8 @@ function banded_qr!(R::BandedMatrix{T}) where T
     QR(R, τ)
 end
 
+banded_qr!(R::AbstractMatrix{T}) where T = banded_qr!(R, zeros(T, min(size(R)...)))
+
 function lmul!(A::QRPackedQ{<:Any,<:BandedMatrix}, B::AbstractVecOrMat)
     require_one_based_indexing(B)
     mA, nA = size(A.factors)

src/banded/linalg.jl

@@ -5,7 +5,7 @@
 function materialize!(L::Ldiv{<:BandedColumnMajor})
     A, B = L.A, L.B
     checksquare(A)
-    ldiv!(factorize(A), dest)
+    ldiv!(factorize(A), B)
 end
 
 function copyto!(dest::AbstractVecOrMat, L::Ldiv{<:BandedRowMajor})

src/generic/broadcast.jl

@@ -353,6 +353,7 @@ function _banded_broadcast!(dest::AbstractMatrix, f, (x, src)::Tuple{Number,Abst
     dest
 end
 
+
 function copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix}})
     (A,) = bc.args
     _banded_broadcast!(dest, bc.f, A, MemoryLayout(typeof(dest)), MemoryLayout(typeof(A)))
@@ -374,6 +375,49 @@ function copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any,
     _banded_broadcast!(dest, bc.f, bc.args, MemoryLayout(typeof(dest)), MemoryLayout.(typeof.(bc.args)))
 end
 
+# override copy in case data has special broadcast
+_default_banded_broadcast(bc::Broadcasted{Style}) where Style = Base.invoke(copy, Tuple{Broadcasted{Style}}, bc)
+
+_banded_broadcast(f, args::Tuple, _...) = _default_banded_broadcast(broadcasted(f, args...))
+_banded_broadcast(f, arg, _...) = _default_banded_broadcast(broadcasted(f, arg))
+
+
+function _banded_broadcast(f, A::AbstractMatrix{T}, ::BandedColumns) where T
+    iszero(f(zero(T))) || return _default_banded_broadcast(broadcasted(f, A))
+    _BandedMatrix(f.(bandeddata(A)), axes(A,1), bandwidths(A)...)
+end
+function _banded_broadcast(f, (src,x)::Tuple{AbstractMatrix{T},Number}, ::BandedColumns) where T
+    iszero(f(zero(T),x)) || return _default_banded_broadcast(broadcasted(f, src,x))
+    _BandedMatrix(f.(bandeddata(src),x), axes(src,1), bandwidths(src)...)
+end
+function _banded_broadcast(f, (x, src)::Tuple{Number,AbstractMatrix{T}}, ::BandedColumns) where T
+    iszero(f(x, zero(T))) || return _default_banded_broadcast(broadcasted(f, x,src))
+    _BandedMatrix(f.(x, bandeddata(src)), axes(src,1), bandwidths(src)...)
+end
+
+
+function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix}})
+    (A,) = bc.args
+    _banded_broadcast(bc.f, A, MemoryLayout(typeof(A)))
+end
+
+function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:Number}})
+    (A,x) = bc.args
+    _banded_broadcast(bc.f, (A, x), MemoryLayout(typeof(A)))
+end
+
+
+function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:Number,<:AbstractMatrix}})
+    (x,A) = bc.args
+    _banded_broadcast(bc.f, (x,A), MemoryLayout(typeof(A)))
+end
+
+
+function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:AbstractMatrix}})
+    _banded_broadcast(bc.f, bc.args, MemoryLayout.(typeof.(bc.args)))
+end
+
+
 _bandwidths(::Number) = (-720,-720)
 _bandwidths(A) = bandwidths(A)
 

src/generic/indexing.jl

@@ -30,6 +30,9 @@ sub_materialize(::AbstractBandedLayout, V) = BandedMatrix(V)
 @inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
 @inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::Colon) = lazy_getindex(A, kr, jr)
 @inline getindex(A::AbstractBandedMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
+@inline getindex(A::AbstractBandedMatrix, kr::AbstractVector, jr::AbstractVector) = lazy_getindex(A, kr, jr)
+@inline getindex(A::AbstractBandedMatrix, kr::Colon, jr::AbstractVector) = lazy_getindex(A, kr, jr)
+@inline getindex(A::AbstractBandedMatrix, kr::AbstractVector, jr::Colon) = lazy_getindex(A, kr, jr)
 
 @inline getindex(A::AbstractMatrix, b::Band) = lazy_getindex(A, b)
 @inline getindex(A::AbstractMatrix, kr::BandRangeType, j::Integer) = lazy_getindex(A, kr, j)

src/generic/matmul.jl

@@ -5,30 +5,19 @@ BroadcastStyle(::BandedStyle, M::ApplyArrayBroadcastStyle{2}) = M
 
 bandwidths(M::Mul) = prodbandwidths(M.args...)
 
-struct BandedMulAddStyle <: AbstractMulAddStyle end
-
-similar(M::Mul{BandedMulAddStyle}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
+bandwidths(M::MulAdd) = prodbandwidths(M.A,M.B)
+similar(M::MulAdd{<:DiagonalLayout,<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
     BandedMatrix{T}(undef, axes, bandwidths(M))
-
-mulapplystyle(A::AbstractBandedLayout, B::AbstractBandedLayout) = BandedMulAddStyle()
-
-### Symmetric
-
-for Lay in (:SymmetricLayout, :HermitianLayout)
-    @eval begin
-        mulapplystyle(::$Lay{<:AbstractBandedLayout}, ::$Lay{<:AbstractBandedLayout}) = BandedMulAddStyle()
-        mulapplystyle(::$Lay{<:AbstractBandedLayout}, ::AbstractBandedLayout) = BandedMulAddStyle()
-        mulapplystyle(::AbstractBandedLayout, ::$Lay{<:AbstractBandedLayout}) = BandedMulAddStyle()
-    end
-end
-
-### Triangular
-
-mulapplystyle(::TriangularLayout{uplo1,unit1,<:AbstractBandedLayout}, ::TriangularLayout{uplo2,unit2,<:AbstractBandedLayout}) where {uplo1,uplo2,unit1,unit2} = BandedMulAddStyle()
-mulapplystyle(::TriangularLayout{uplo,unit,<:AbstractBandedLayout}, ::AbstractBandedLayout) where {uplo,unit} = BandedMulAddStyle()
-mulapplystyle(::AbstractBandedLayout, ::TriangularLayout{uplo,unit,<:AbstractBandedLayout}) where {uplo,unit} = BandedMulAddStyle()
-
-
+similar(M::MulAdd{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
+    BandedMatrix{T}(undef, axes, bandwidths(M))
+similar(M::MulAdd{<:AbstractBandedLayout,<:DiagonalLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T = 
+    BandedMatrix{T}(undef, axes, bandwidths(M))
+similar(M::MulAdd{<:SymmetricLayout{<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
+    BandedMatrix{T}(undef, axes, bandwidths(M))    
+similar(M::MulAdd{<:HermitianLayout{<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
+    BandedMatrix{T}(undef, axes, bandwidths(M))    
+similar(M::MulAdd{<:TriangularLayout{uplo,trans,<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where {uplo,trans,T} =
+    BandedMatrix{T}(undef, axes, bandwidths(M))    
 ##
 # BLAS routines
 ##
@@ -288,3 +277,20 @@ end
 #     end
 #     C
 # end
+
+
+
+###
+# Special Fill Diagonal
+####
+
+function materialize!(M::MatMulMatAdd{<:DiagonalLayout{<:AbstractFillLayout},<:AbstractBandedLayout})
+    M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
+    M.C
+end
+
+function materialize!(M::MatMulMatAdd{<:AbstractBandedLayout,<:DiagonalLayout{<:AbstractFillLayout}})
+    M.C .= (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
+    M.C
+end
+

src/interfaceimpl.jl

@@ -37,9 +37,9 @@ inbands_setindex!(J::SymTridiagonal, v, k::Integer, j::Integer) =
 isbanded(::Tridiagonal) = true
 bandwidths(::Tridiagonal) = (1,1)
 
-isbanded(K::Kron{<:Any,2}) = all(isbanded, K.arrays)
+isbanded(K::Kron{<:Any,2}) = all(isbanded, K.args)
 function bandwidths(K::Kron{<:Any,2})
-    A,B = K.arrays
+    A,B = K.args
     (size(B,1)*bandwidth(A,1) + max(0,size(B,1)-size(B,2))*size(A,1)   + bandwidth(B,1),
         size(B,2)*bandwidth(A,2) + max(0,size(B,2)-size(B,1))*size(A,2) + bandwidth(B,2))
 end
@@ -69,23 +69,22 @@ end
 # needed for ∞-dimensional banded linear algebra
 ###
 
-struct PaddedMulAddStyle <: AbstractMulAddStyle end
-mulapplystyle(::AbstractBandedLayout, ::VcatLayout{<:Tuple{<:Any,ZerosLayout}}) = PaddedMulAddStyle()
-
-function similar(M::Mul{PaddedMulAddStyle}, ::Type{T}, axes) where T
-    A,x = M.args
-    xf,_ = x.arrays
+function similar(M::MulAdd{<:AbstractBandedLayout,<:ApplyLayout{typeof(vcat),<:Tuple{<:Any,ZerosLayout}}}, ::Type{T}, axes) where T
+    A,x = M.A,M.B
+    xf,_ = x.args
     n = max(0,min(length(xf) + bandwidth(A,1),length(M)))
     Vcat(Vector{T}(undef, n), Zeros{T}(size(A,1)-n))
 end
 
-function materialze!(M::MatMulVecAdd{<:AbstractBandedLayout,<:VcatLayout{<:Tuple{<:Any,ZerosLayout}},<:VcatLayout{<:Tuple{<:Any,ZerosLayout}}})
-    α,A,x,β,y = M.α,Μ.A,Μ.B,Μ.β,M.C
+function materialize!(M::MatMulVecAdd{<:AbstractBandedLayout,
+                                    <:ApplyLayout{typeof(vcat),<:Tuple{<:Any,ZerosLayout}},
+                                     <:ApplyLayout{typeof(vcat),<:Tuple{<:Any,ZerosLayout}}})
+    α,A,x,β,y = M.α,M.A,M.B,M.β,M.C
     length(y) == size(A,1) || throw(DimensionMismatch())
     length(x) == size(A,2) || throw(DimensionMismatch())
 
-    ỹ,_ = y.arrays
-    x̃,_ = x.arrays
+    ỹ,_ = y.args
+    x̃,_ = x.args
 
     length(ỹ) ≥ min(length(M),length(x̃)+bandwidth(A,1)) ||
         throw(InexactError("Cannot assign non-zero entries to Zero"))
@@ -101,9 +100,9 @@ end
 # MulMatrix
 ###
 
-bandwidths(M::MulMatrix) = bandwidths(M.applied)
+bandwidths(M::MulMatrix) = bandwidths(Applied(M))
 isbanded(M::Mul) = all(isbanded, M.args)
-isbanded(M::MulMatrix) = isbanded(M.applied)
+isbanded(M::MulMatrix) = isbanded(Applied(M))
 
 const MulBandedLayout = MulLayout{<:Tuple{Vararg{<:AbstractBandedLayout}}}
 
@@ -125,11 +124,11 @@ subarraylayout(M::MulBandedLayout, ::Type{<:Tuple{Vararg{AbstractUnitRange}}}) =
 # Concat banded matrix
 ######
 
-bandwidths(M::Hcat) = (bandwidth(M.arrays[1],1),sum(size.(M.arrays[1:end-1],2)) + bandwidth(M.arrays[end],2))
-isbanded(M::Hcat) = all(isbanded, M.arrays)
+bandwidths(M::Hcat) = (bandwidth(M.args[1],1),sum(size.(M.args[1:end-1],2)) + bandwidth(M.args[end],2))
+isbanded(M::Hcat) = all(isbanded, M.args)
 
-bandwidths(M::Vcat) = (sum(size.(M.arrays[1:end-1],1)) + bandwidth(M.arrays[end],1), bandwidth(M.arrays[1],2))
-isbanded(M::Vcat) = all(isbanded, M.arrays)
+bandwidths(M::Vcat) = (sum(size.(M.args[1:end-1],1)) + bandwidth(M.args[end],1), bandwidth(M.args[1],2))
+isbanded(M::Vcat) = all(isbanded, M.args)
 
 
 const HcatBandedMatrix{T,N} = Hcat{T,NTuple{N,BandedMatrix{T,Matrix{T},OneTo{Int}}}}
@@ -149,3 +148,30 @@ hcat(A::BandedMatrix, B::AbstractMatrix...) = Matrix(Hcat(A, B...))
 vcat(A::BandedMatrix...) = BandedMatrix(Vcat(A...))    
 vcat(A::BandedMatrix, B::AbstractMatrix...) = Matrix(Vcat(A, B...))    
 
+
+
+#######
+# CachedArray
+#######
+
+bandwidths(B::CachedMatrix) = bandwidths(B.data)
+isbanded(B::CachedMatrix) = isbanded(B.data)
+
+function resizedata!(B::CachedMatrix{T,BandedMatrix{T,Matrix{T},OneTo{Int}}}, n::Integer, m::Integer) where T<:Number
+    @boundscheck checkbounds(Bool, B, n, m) || throw(ArgumentError("Cannot resize beyound size of operator"))
+
+    # increase size of array if necessary
+    olddata = B.data
+    ν,μ = size(olddata)
+    n,m = max(ν,n), max(μ,m)
+    if (ν,μ) ≠ (n,m)
+        B.data = BandedMatrix{T}(undef, (n,m), bandwidths(olddata))
+        B.data.data[:,1:μ] .= olddata.data
+        bd = bandeddata(B.array)
+        B.data.data[:,μ+1:end] .= view(bd,:,μ+1:m)
+    end
+
+    B.datasize = (n,m)
+
+    B
+end

test/test_banded.jl

@@ -343,6 +343,7 @@ Base.similar(::MyMatrix, ::Type{T}, m::Int, n::Int) where T = MyMatrix{T}(undef,
             @test @inferred(similar(banded,5,5)) isa BandedMatrix{Int32, Matrix{Int32}}
             @test @inferred(similar(banded,5,5,1,1)) isa BandedMatrix{Int32, Matrix{Int32}}
 
+
             banded = convert(BandedMatrix{<:, MyMatrix}, brand(Int32, 10, 12, 1, 2))
             @test banded isa BandedMatrix{Int32, MyMatrix{Int32}}
             @test @inferred(similar(banded)) isa BandedMatrix{Int32, MyMatrix{Int32}}