Add support for BlockBidiagonal (#77)

* Add support for BlockBidiagonal

* v0.8.8

* Fix I - A

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.8.7"
+version = "0.8.8"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,8 +17,8 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 ArrayLayouts = "0.3.3"
-BandedMatrices = "0.15.10"
-BlockArrays = "0.12.7"
+BandedMatrices = "0.15.14"
+BlockArrays = "0.12.9"
 FillArrays = "0.8.10"
 MatrixFactorizations = "0.4.1"
 julia = "1.2"

src/BlockBandedMatrices.jl

@@ -45,7 +45,7 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
                         _banded_colval, _banded_rowval, _banded_nzval # for sparse
 
 export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockbandwidth, blockbandwidths,
-        subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, isblockbanded
+        subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, BlockBidiagonal, isblockbanded
 
 
 const Block1 = Block{1,Int}

src/interfaceimpl.jl

@@ -36,10 +36,12 @@ function sizes_from_blocks(A::Diagonal, _)
 end    
 
 
-# Block Tridiagonal
+# Block Bi/Tridiagonal
 const BlockTridiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Tridiagonal{VT}}
+const BlockBidiagonal{T,VT<:Matrix{T}} = BlockMatrix{T,<:Bidiagonal{VT}}
 
 BlockTridiagonal(A,B,C) = mortar(Tridiagonal(A,B,C))
+BlockBidiagonal(A, B, uplo) = mortar(Bidiagonal(A,B,uplo))
 
 function sizes_from_blocks(A::Tridiagonal, _) 
     # for k = 1:length(A.du)
@@ -51,9 +53,26 @@ function sizes_from_blocks(A::Tridiagonal, _)
     (size.(A.d, 1), size.(A.d,2))
 end
 
+function sizes_from_blocks(A::Bidiagonal, _) 
+    # for k = 1:length(A.du)
+    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
+    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
+    # end
+    (size.(A.dv, 1), size.(A.dv,2))
+end
+
 blockbandwidths(A::BlockArray) = bandwidths(A.blocks)
 isblockbanded(A::BlockArray) = isbanded(A.blocks)
 
+@inline function getblock(block_arr::BlockBidiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
+    @boundscheck blockcheckbounds(block_arr, K, J)
+    l,u = blockbandwidths(block_arr)
+    -l ≤ (J-K) ≤ u || return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
+    block_arr.blocks[K,J]
+end
+
 @inline function getblock(block_arr::BlockTridiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
     @boundscheck blockcheckbounds(block_arr, K, J)
     abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
@@ -70,7 +89,15 @@ for op in (:-, :+)
         end
         function $op(λ::UniformScaling, A::BlockTridiagonal) 
             checksquareblocks(A)
-            mortar(Tridiagonal(A.blocks.dl, broadcast($op, Ref(λ), A.blocks.d), A.blocks.du))
+            mortar(Tridiagonal(broadcast($op,A.blocks.dl), broadcast($op, Ref(λ), A.blocks.d), broadcast($op,A.blocks.du)))
+        end
+        function $op(A::BlockBidiagonal, λ::UniformScaling) 
+            checksquareblocks(A)
+            mortar(Bidiagonal(broadcast($op, A.blocks.dv, Ref(λ)), A.blocks.ev, A.blocks.uplo))
+        end
+        function $op(λ::UniformScaling, A::BlockBidiagonal) 
+            checksquareblocks(A)
+            mortar(Bidiagonal(broadcast($op, Ref(λ), A.blocks.dv), broadcast($op,A.blocks.ev), A.blocks.uplo))
         end
     end
 end

test/test_misc.jl

@@ -21,6 +21,25 @@ import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedStyle, bandeddata
     @test BandedMatrix(V) == V
 end
 
+@testset "Block Bidiagonal" begin
+    Bu = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :U)
+    Bl = BlockBidiagonal(fill([1 2],4), fill([3 4],3), :L)
+    @test blockbandwidths(Bu) == (0,1)
+    @test blockbandwidths(Bl) == (1,0)
+    @test isblockbanded(Bu)
+    @test isblockbanded(Bl)
+    @test Bu[Block(1,1)] == Bl[Block(1,1)] == [1 2]
+    @test @inferred(Bu[Block(1,2)]) == @inferred(Bl[Block(2,1)]) == [3 4]
+    @test @inferred(getblock(Bu,1,3)) == @inferred(Bu[Block(1,3)]) == [0 0]
+    @test_throws DimensionMismatch Bu+I
+    Bu = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :U)
+    Bl = BlockBidiagonal(fill([1 2; 1 2],4), fill([3 4; 3 4],3), :L)
+    @test Bu+I == I+Bu == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) + I
+    @test Bl+I == I+Bl == mortar(Bidiagonal(fill([2 2; 1 3],4), fill([3 4; 3 4],3), :L)) == Matrix(Bl) + I
+    @test Bu-I == mortar(Bidiagonal(fill([0 2; 1 1],4), fill([3 4; 3 4],3), :U)) == Matrix(Bu) - I
+    @test I-Bu == mortar(Bidiagonal(fill([0 -2; -1 -1],4), fill(-[3 4; 3 4],3), :U)) == I - Matrix(Bu)
+end
+
 @testset "Block Tridiagonal" begin
     A = BlockTridiagonal(fill([1 2],3), fill([3 4],4), fill([4 5],3))
     @test blockbandwidths(A) == (1,1)
@@ -31,4 +50,6 @@ end
     @test_throws DimensionMismatch A+I
     A = BlockTridiagonal(fill([1 2; 1 2],3), fill([3 4; 3 4],4), fill([4 5; 4 5],3))
     @test A+I == I+A == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([4 4; 3 5],4), fill([4 5; 4 5],3))) == Matrix(A) + I
+    @test A-I == mortar(Tridiagonal(fill([1 2; 1 2],3), fill([2 4; 3 3],4), fill([4 5; 4 5],3))) == Matrix(A) - I
+    @test I-A == mortar(Tridiagonal(fill(-[1 2; 1 2],3), fill([-2 -4; -3 -3],4), fill(-[4 5; 4 5],3))) == I - Matrix(A)
 end