Support block banded cholesky (#130)

* Support block banded cholesky

* Update test_linalg.jl

* v0.11.7

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.11.6"
+version = "0.11.7"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -13,11 +13,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ArrayLayouts = "0.8.2"
-BandedMatrices = "0.16.8, 0.17"
-BlockArrays = "0.16.16"
+ArrayLayouts = "0.8.7"
+BandedMatrices = "0.17"
+BlockArrays = "0.16.18"
 FillArrays = "0.13"
-MatrixFactorizations = "0.8.5, 0.9"
+MatrixFactorizations = "0.9"
 julia = "1.6"
 
 [extras]

src/BlockBandedMatrices.jl

@@ -33,7 +33,7 @@ import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal,
                         Block, BlockSlice, unblock, block, blockindex,
                         _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                         blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
-                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle, BlockSlices
+                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle, BlockSlices, _blockkron
 
 import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                         inbands_setindex!, inbands_getindex, banded_setindex!,

src/interfaceimpl.jl

@@ -164,4 +164,17 @@ end
 
 # ambiguity
 sub_materialize(::AbstractBandedLayout, V, ::Tuple{BlockedUnitRange,Base.OneTo{Int}}) = BandedMatrix(V)
-sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},BlockedUnitRange}) = BandedMatrix(V)
+sub_materialize(::AbstractBandedLayout, V, ::Tuple{Base.OneTo{Int},BlockedUnitRange}) = BandedMatrix(V)
+
+
+
+#####
+# BlockKron
+#####
+
+isblockbanded(K::BlockKron) = isbanded(first(K.args))
+isbandedblockbanded(K::BlockKron) = all(isbanded, K.args)
+blockbandwidths(K::BlockKron) = bandwidths(first(K.args))
+subblockbandwidths(K::BlockKron) = bandwidths(last(K.args))
+
+_blockkron(::Tuple{Vararg{AbstractBandedLayout}}, A) = BandedBlockBandedMatrix(BlockKron(A...))

src/linalg.jl

@@ -145,3 +145,11 @@ strides(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) whe
     (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block.block)])
 
 MemoryLayout(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU} = ColumnMajor()
+
+
+
+###
+# cholesky
+##
+
+ArrayLayouts._cholesky(::AbstractBandedBlockBandedLayout, axes, A::Symmetric, piv::ArrayLayouts.CNoPivot=ArrayLayouts.CNoPivot(); check::Bool = true) = cholesky!(Symmetric(BlockBandedMatrix(parent(A)), Symbol(A.uplo)), piv; check = check)

test/test_linalg.jl

@@ -175,3 +175,13 @@ end
     B.data .= randn.()
     @test A*B ≈ Matrix(A)*Matrix(B)
 end
+
+
+@testset "Cholesky" begin
+    n = 5
+    D² = SymTridiagonal(fill(2.0,n), -ones(n-1))
+    Δ = blockkron(D²,I(n)) + blockkron(I(n),D²)
+    @test Δ isa BandedBlockBandedMatrix
+    @test cholesky(Symmetric(Δ)).U ≈ cholesky(Matrix(Δ)).U
+    @test cholesky(Symmetric(Δ,:L)).U ≈ cholesky(Matrix(Δ)).U
+end