Dl/layoutmatrix (#66)

* Support new ArrayLayouts

* layoutmatrix macro

* Update BlockBandedMatrices.jl

* Update .travis.yml

* _qr

* Update Project.toml

* Update Project.toml

* Move code to BlockArrays

* update for new blockarrays

* updates

Author: dlfivefifty

.gitignore

@@ -1,2 +1,3 @@
 
 .DS_Store
+Manifest.toml

.travis.yml

@@ -5,7 +5,6 @@ os:
   - osx
 julia:
   - 1.2
-  - 1.3
   - 1.4
   - nightly
 matrix:

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.7.2"
+version = "0.8"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,9 +18,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "0.1"
-BandedMatrices = "0.14"
-BlockArrays = "0.11"
+ArrayLayouts = "0.2.3"
+BandedMatrices = "0.15.2"
+BlockArrays = "0.12"
 FillArrays = "0.8"
-MatrixFactorizations = "0.2, 0.3"
+MatrixFactorizations = "0.3.1"
 julia = "1.2"

examples/sharedarrays_setup.jl

@@ -2,7 +2,7 @@ using Pkg
 
 Pkg.activate(homedir() * "/Documents/Coding/gpublockbanded")
 using BandedMatrices, BlockBandedMatrices, SharedArrays, ArrayLayouts, BlockArrays, FillArrays
-import BlockBandedMatrices: _BandedBlockBandedMatrix, BandedBlockBandedSizes, BlockSizes, blockcolrange
+import BlockBandedMatrices: _BandedBlockBandedMatrix, BandedBlockBandedSizes, BlockSizes, blockcolsupport
 import BandedMatrices: AbstractBandedMatrix, bandwidths, BandedStyle
 import Base: getindex, size
 
@@ -28,7 +28,7 @@ Base.BroadcastStyle(::Type{<:Eye}) = BandedStyle()
 
 
 function blockcolbuild!(ret, A, JR)
-    for J in JR, K in blockcolrange(A,J)
+    for J in JR, K in blockcolsupport(A,J)
         view(ret,K,J) .= view(A,K,J)
     end
 end

src/AbstractBlockBandedMatrix.jl

@@ -6,10 +6,11 @@
 # in addition to the banded matrix interface
 ####
 
-abstract type AbstractBlockBandedLayout <: MemoryLayout end
+abstract type AbstractBlockBandedLayout <: AbstractBlockLayout end
+abstract type AbstractBandedBlockBandedLayout <: AbstractBlockBandedLayout end
 
-struct BandedBlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
-struct BandedBlockBandedRows{LAY} <: AbstractBlockBandedLayout end
+struct BandedBlockBandedColumns{LAY} <: AbstractBandedBlockBandedLayout end
+struct BandedBlockBandedRows{LAY} <: AbstractBandedBlockBandedLayout end
 struct BlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
 struct BlockBandedRows{LAY} <: AbstractBlockBandedLayout end
 
@@ -22,7 +23,6 @@ transposelayout(::BandedBlockBandedColumnMajor) = BandedBlockBandedRowMajor()
 transposelayout(::BandedBlockBandedRowMajor) = BandedBlockBandedColumnMajor()
 transposelayout(::BlockBandedColumnMajor) = BlockBandedRowMajor()
 transposelayout(::BlockBandedRowMajor) = BlockBandedColumnMajor()
-conjlayout(::Type{<:Complex}, M::AbstractBlockBandedLayout) = ConjLayout(M)
 
 
 
@@ -56,22 +56,21 @@ blockbandrange(A::AbstractMatrix) = -blockbandwidth(A,1):blockbandwidth(A,2)
 
 
 # start/stop indices of the i-th column/row, bounded by actual matrix size
-@inline blockcolstart(A::AbstractVecOrMat, i::Integer) = Block(max(i-colblockbandwidth(A,2)[i], 1))
-@inline  blockcolstop(A::AbstractVecOrMat, i::Integer) = Block(max(min(i+colblockbandwidth(A,1)[i], blocksize(A, 1)), 0))
-@inline blockrowstart(A::AbstractVecOrMat, i::Integer) = Block(max(i-rowblockbandwidth(A,1)[i], 1))
-@inline  blockrowstop(A::AbstractVecOrMat, i::Integer) = Block(max(min(i+rowblockbandwidth(A,2)[i], blocksize(A, 2)), 0))
-
-for Func in (:blockcolstart, :blockcolstop, :blockrowstart, :blockrowstop)
+@inline blockbanded_blockcolstart(A, i::Integer) = Block(max(i-colblockbandwidth(A,2)[i], 1))
+@inline blockbanded_blockcolstop(A, i::Integer) = Block(max(min(i+colblockbandwidth(A,1)[i], blocksize(A, 1)), 0))
+@inline blockbanded_blockrowstart(A, i::Integer) = Block(max(i-rowblockbandwidth(A,1)[i], 1))
+@inline blockbanded_blockrowstop(A, i::Integer) = Block(max(min(i+rowblockbandwidth(A,2)[i], blocksize(A, 2)), 0))
+for Func in (:blockbanded_blockcolstart, :blockbanded_blockcolstop, :blockbanded_blockrowstart, :blockbanded_blockrowstop)
     @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
 end
 
-@inline blockcolrange(A::AbstractVecOrMat, i) = blockcolstart(A,i):blockcolstop(A,i)
-@inline blockrowrange(A::AbstractVecOrMat, i) = blockrowstart(A,i):blockrowstop(A,i)
+@inline blockcolsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
+@inline blockrowsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)
 
 
 # length of i-the column/row
-@inline blockcollength(A::AbstractVecOrMat, i) = max(Int(blockcolstop(A, i)) - Int(blockcolstart(A, i)) + 1, 0)
-@inline blockrowlength(A::AbstractVecOrMat, i) = max(Int(blockrowstop(A, i)) - Int(blockrowstart(A, i)) + 1, 0)
+@inline blockcollength(A, i) = max(Int(blockcolstop(A, i)) - Int(blockcolstart(A, i)) + 1, 0)
+@inline blockrowlength(A, i) = max(Int(blockrowstop(A, i)) - Int(blockrowstart(A, i)) + 1, 0)
 
 # this gives the block bandwidth in each block column/row
 @inline colblockbandwidths(A::AbstractMatrix) = Fill.(blockbandwidths(A), blocksize(A,2))
@@ -107,17 +106,22 @@ function _firstblockdiagcol(A::AbstractMatrix, j::Int)
 end
 
 
-## BlockSlice1 is a conveneience for views
-
-const BlockSlice1 = BlockSlice{Block{1,Int}}
 
 
 ######################################
 #  RaggedMatrix interface
 ######################################
 
-@inline colstart(A::AbstractBlockBandedMatrix, i::Integer) =
-    isempty(axes(A,1)) ? 1 : first(axes(A,1)[blockcolstart(A,findblock(axes(A,2),i))])
+@inline function colstart(A::AbstractBlockBandedMatrix, i::Integer)
+    bs = blockcolstart(A,findblock(axes(A,2),i))
+    if isempty(axes(A,1)) 
+        1
+    elseif Int(bs) ≤ blocksize(A, 2)
+        first(axes(A,1)[bs])
+    else
+        size(A,1)+1
+    end
+end
 
 @inline function colstop(A::AbstractBlockBandedMatrix, i::Integer)
     CS = blockcolstop(A,findblock(axes(A,2),i))

src/BandedBlockBandedMatrix.jl

@@ -153,7 +153,7 @@ BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
 function BandedBlockBandedMatrix{T,Blocks,RR}(A::AbstractMatrix, axes::NTuple{2,AbstractUnitRange{Int}}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,Blocks,RR<:AbstractUnitRange{Int}}
     ret = BandedBlockBandedMatrix{T,Blocks,RR}(Zeros{T}(size(A)), axes, lu, λμ)
     L,M = λμ 
-    for J = Block.(1:blocksize(ret, 2)), K = blockcolrange(ret, Int(J))
+    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
         kr, jr = axes[1][K], axes[2][J]
 
         # We have the correct block - now we need to only add the entries from
@@ -202,7 +202,7 @@ function sparse(A::BandedBlockBandedMatrix)
     i = Vector{Int}()
     j = Vector{Int}()
     z = Vector{eltype(A)}()
-    for J = blockaxes(A,2), K = blockcolrange(A, J)
+    for J = blockaxes(A,2), K = blockcolsupport(A, J)
         B = view(A, K, J)
         ĩ = _banded_rowval(B)
         j̃ = _banded_colval(B)
@@ -221,7 +221,7 @@ end
 ################################
 
 bandedblockbandedcolumns(L::AbstractColumnMajor) = BandedBlockBandedColumnMajor()
-bandedblockbandedcolumns(_) = UnknownLayout()
+bandedblockbandedcolumns(_) = BandedBlockBandedColumns{UnknownLayout}()
 
 MemoryLayout(::Type{<:BandedBlockBandedMatrix{<:Any,BLOCKS}}) where BLOCKS = 
     bandedblockbandedcolumns(MemoryLayout(BLOCKS))

src/BlockBandedMatrices.jl

@@ -18,18 +18,21 @@ import LinearAlgebra: UniformScaling, isdiag, rmul!, lmul!, ldiv!, rdiv!,
 import LinearAlgebra.BLAS: BlasInt, BlasFloat, @blasfunc, libblas, BlasComplex, BlasReal
 import LinearAlgebra.LAPACK: chktrans, chkdiag, liblapack, chklapackerror, checksquare, chkstride1,
                     chkuplo
-import MatrixFactorizations: ql, ql!, QLPackedQ
+import MatrixFactorizations: ql, ql!, _ql, QLPackedQ
 import SparseArrays: sparse
 
-import ArrayLayouts: @lazymul, MatMulMatAdd, MatMulVecAdd, BlasMatLmulVec,
+import ArrayLayouts: BlasMatLmulVec,
                     triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
-                    triangulardata, sublayout, @lazyldiv, @lazylmul,
+                    triangulardata, sublayout, 
                     AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
-                    DiagonalLayout, materialize!, MulAdd, mul, colsupport, rowsupport
+                    DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
+                    _qr, _factorize, _copyto!
 
 import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
                         Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
-                        _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks
+                        _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
+                        blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
+                        AbstractBlockLayout
 
 import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                         inbands_setindex!, inbands_getindex, banded_setindex!,
@@ -47,7 +50,10 @@ export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockband
 
 const Block1 = Block{1,Int}
 const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
-const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}        
+const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}    
+
+blockcolrange(A...) = blockcolsupport(A...)
+blockrowrange(A...) = blockrowsupport(A...)
 
 include("AbstractBlockBandedMatrix.jl")
 include("broadcast.jl")

src/BlockSkylineMatrix.jl

@@ -160,7 +160,7 @@ BlockSkylineMatrix
 
 function BlockSkylineMatrix{T}(A::AbstractMatrix, block_sizes::BlockSkylineSizes) where T
     ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
-    for J = Block.(1:blocksize(ret, 2)), K = blockcolrange(ret, Int(J))
+    for J = blockaxes(ret,2), K = blockcolsupport(ret, Int(J))
         kr, jr = getindex.(block_sizes.axes, (K, J))
         view(ret, K, J) .= view(A, kr, jr)
     end
@@ -170,7 +170,7 @@ end
 function BlockSkylineMatrix{T}(A::AbstractBlockBandedMatrix, block_sizes::BlockSkylineSizes) where T
     ret = BlockSkylineMatrix(Zeros{T}(size(A)), block_sizes)
     blockisequal(axes(A), block_sizes.axes) || throw(ArgumentError())
-    for J = blockaxes(ret,2), K = blockcolrange(ret, Int(J))
+    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
         view(ret, K, J) .= view(A, K, J)
     end
     ret
@@ -234,7 +234,7 @@ function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
     block_sizes = BlockSkylineSizes(axes(A), colblockbandwidths(A)...)
 
     ret = BlockSkylineMatrix{eltype(A)}(undef, block_sizes)
-    for J = blockaxes(ret,2), K = blockcolrange(ret, J)
+    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
         view(ret, K, J) .= view(A, K, J)
     end
     ret

src/blockskylineqr.jl

@@ -57,11 +57,9 @@ function ql!(A::BlockBandedMatrix{T}) where T
     QL(A,τ.blocks)
 end
 
-qr(A::BlockBandedMatrix) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
-ql(A::BlockBandedMatrix) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))
-
-qr(A::BandedBlockBandedMatrix) = qr(BlockBandedMatrix(A))
-ql(A::BandedBlockBandedMatrix) = ql(BlockBandedMatrix(A))
+_qr(::AbstractBlockBandedLayout, _, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
+_ql(::AbstractBlockBandedLayout, _, A) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))
+_factorize(::AbstractBlockBandedLayout, _, A) = qr(A)
 
 function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
     Q = parent(adjQ)
@@ -123,5 +121,3 @@ function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractMatrix)
     lmul!(adjoint(A.Q), B)
     materialize!(Ldiv(LowerTriangular(A.factors), B))
 end
-
-\(A::AbstractBlockBandedMatrix, b::AbstractVecOrMat) = qr(A)\b  # use QR because LU would be a _mess_ to implement

src/broadcast.jl

@@ -79,7 +79,7 @@ function blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
 end
 
 
-copyto!(dest::AbstractMatrix, src::AbstractBlockBandedMatrix) = blockbanded_copyto!(dest, src)
+_copyto!(_, ::AbstractBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) = blockbanded_copyto!(dest, src)
 
 
 function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})
@@ -100,9 +100,8 @@ end
 function blockbanded_fill!(B::AbstractMatrix{T}, x) where T
     x == zero(T) || throw(BandError(B))
 
-    M,N = blocksize(B)
-    for J = 1:N, K = blockcolrange(B,J)
-        fill!(view(B,K,Block(J)), x)
+    for J = blockaxes(B,2), K = blockcolsupport(B,J)
+        fill!(view(B,K,J), x)
     end
     B
 end
@@ -111,19 +110,17 @@ end
 function blockbanded_rmul!(B::AbstractMatrix{T}, x::Number) where T
     x == zero(T) || throw(BandError(B))
 
-    M,N = blocksize(B)
-    for J = 1:N, K = blockcolrange(B,J)
-        rmul!(view(B,K,Block(J)), x)
+    for J = blockaxes(B,2), K = blockcolsupport(B,J)
+        rmul!(view(B,K,J), x)
     end
     B
 end
 
 function blockbanded_lmul!(x::Number, B::AbstractMatrix{T}) where T
     x == zero(T) || throw(BandError(B))
 
-    M,N = blocksize(B)
-    for J = 1:N, K = blockcolrange(B,J)
-        lmul!(x, view(B,K,Block(J)))
+    for J = blockaxes(B,2), K = blockcolsupport(B,J)
+        lmul!(x, view(B,K,J))
     end
     B
 end
@@ -205,7 +202,7 @@ similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(\), <:Tuple{<:
 function blockbanded_axpy!(a, X::AbstractMatrix, Y::AbstractMatrix)
     size(X) == size(Y) || throw(DimensionMismatch())
 
-    for J=Block(1):Block(blocksize(X,2)), K=blockcolrange(X,J)
+    for J = blockaxes(X,2), K = blockcolsupport(X,J)
         view(Y,K,J) .= a .* view(X,K,J) .+ view(Y,K,J)
     end
     Y