Add more blockbanded Kron tests

Author: dlfivefifty

Project.toml

@@ -14,7 +14,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 ArrayLayouts = "0.8.7"
-BandedMatrices = "0.17"
+BandedMatrices = "0.17.2"
 BlockArrays = "0.16.18"
 FillArrays = "0.13"
 MatrixFactorizations = "0.9"

test/test_misc.jl

@@ -29,6 +29,23 @@ function Base.getindex(A::MyBlockBandedMatrix, k::Int, j::Int)
     A.A[k,j]
 end
 
+struct FiniteDifference{T} <: AbstractBandedMatrix{T}
+    n::Int
+end
+
+FiniteDifference(n) = FiniteDifference{Float64}(n)
+
+Base.getindex(F::FiniteDifference{T}, k::Int, j::Int) where T =
+    if k == j
+        -2*one(T)*F.n^2
+    elseif abs(k-j) == 1
+        one(T)*F.n^2
+    else
+        zero(T)
+    end
+
+BandedMatrices.bandwidths(F::FiniteDifference) = (1,1)
+Base.size(F::FiniteDifference) = (F.n,F.n)
 
 @testset "Interfaces" begin
     @testset "MyBandedBlockBandedMatrix" begin
@@ -182,4 +199,40 @@ end
             @test blockbandwidths(B) == (4,0)
         end
     end
+
+    @testset "Block banded Kron" begin
+        n = 10
+        h = 1/n
+        D² = BandedMatrix(0 => Fill(-2,n), 1 => Fill(1,n-1), -1 => Fill(1,n-1))/h^2
+
+        @time D_xx = BandedBlockBandedMatrix(BlockKron(D², Eye(n)))
+        @time D_yy = BandedBlockBandedMatrix(BlockKron(Eye(n),D²))
+        @test D_xx == blockkron(D², Eye(n))
+        @time Δ = D_xx + D_yy
+
+        @test Δ isa BandedBlockBandedMatrix
+        @test blockbandwidths(Δ) == subblockbandwidths(Δ) == (1,1)
+        @test Δ == blockkron(Matrix(D²), Matrix(I,n,n)) + blockkron(Matrix(I,n,n), Matrix(D²))
+
+        n = 10
+        D² = FiniteDifference(n)
+        D̃_xx = BlockKron(D², Eye(n))
+        @test blockbandwidths(D̃_xx) == (1,1)
+        @test subblockbandwidths(D̃_xx) == (0,0)
+
+        V = view(D̃_xx, Block(1,1))
+        @test bandwidths(V) == (0,0)
+
+        @test BandedBlockBandedMatrix(D̃_xx) ≈ D_xx
+
+        D̃_yy = BlockKron(Eye(n), D²)
+        @test blockbandwidths(D̃_yy) == (0,0)
+        @test subblockbandwidths(D̃_yy) == (1,1)
+
+        V = view(D̃_yy, Block(1,1))
+        @test bandwidths(V) == (1,1)
+
+        @test BandedBlockBandedMatrix(D̃_yy) ≈ D_yy
+    end
 end
+