BBBArrowheadMatrices module (#63)

* BBBArrowheadMatrices module

* tests pass

* Update test_arrowhead.jl

* banded interface

* eigvals via banded matrices

* Update test_arrowhead.jl

* v0.4.1

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "PiecewiseOrthogonalPolynomials"
 uuid = "4461d12d-4663-4550-8580-cb764c85e20f"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.4"
+version = "0.4.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/PiecewiseOrthogonalPolynomials.jl

@@ -10,15 +10,15 @@ import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBa
 import LazyArrays: paddeddata, AbstractLazyLayout
 import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
 import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv, show, summary
-import LinearAlgebra: BlasInt
 import InfiniteArrays: OneToInf
 import FillArrays: AbstractFill
-import MatrixFactorizations: reversecholcopy
 import FillArrays: SquareEye
 
-export PiecewisePolynomial, ContinuousPolynomial, DirichletPolynomial, Derivative, Block, weaklaplacian, grammatrix
+export PiecewisePolynomial, ContinuousPolynomial, DirichletPolynomial, Derivative, Block, weaklaplacian, grammatrix, BBBArrowheadMatrix
 
 include("arrowhead.jl")
+using .BBBArrowheadMatrices
+
 include("piecewisepolynomial.jl")
 include("continuouspolynomial.jl")
 include("dirichlet.jl")

src/arrowhead.jl

@@ -1,3 +1,20 @@
+module BBBArrowheadMatrices
+using LinearAlgebra, BlockArrays, BlockBandedMatrices, BandedMatrices, MatrixFactorizations, LazyBandedMatrices, LazyArrays, ArrayLayouts, InfiniteArrays, FillArrays
+import ArrayLayouts: MemoryLayout, sublayout, sub_materialize, symmetriclayout, transposelayout, SymmetricLayout, HermitianLayout, TriangularLayout, layout_getindex, materialize!, MatLdivVec, AbstractStridedLayout, triangulardata, MatMulMatAdd, MatMulVecAdd, _fill_lmul!, layout_replace_in_print_matrix
+import BandedMatrices: isbanded, bandwidths
+import BlockArrays: BlockSlice, block, blockindex, blockvec
+import BlockBandedMatrices: subblockbandwidths, blockbandwidths, AbstractBandedBlockBandedLayout, AbstractBandedBlockBandedMatrix
+import Base: size, axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail, _sum, inv, show, summary
+import LazyArrays: paddeddata, AbstractLazyLayout
+import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
+import LinearAlgebra: BlasInt, eigvals
+import MatrixFactorizations: reversecholcopy
+import FillArrays: AbstractFill
+import FillArrays: SquareEye
+import InfiniteArrays: OneToInf
+
+
+export BBBArrowheadMatrix
 
 """
 BBBArrowheadMatrix
@@ -498,4 +515,22 @@ for Tri in (:LowerTriangular, :UnitLowerTriangular)
             getfield(F, d)
         end
     end
-end
+end
+
+####
+# banded interface
+#####
+bbbbandwidth(k, A) = bandwidth(A, k)
+bbbbandwidth(k, A, C...) = +(size(A,k), size.(Base.front(C), k)...) + bandwidth(last(C),k)
+bandwidths(A::BBBArrowheadMatrix) = bbbbandwidth(1, A.A, A.C...),  bbbbandwidth(2, A.A, A.B...)
+isbanded(A::BBBArrowheadMatrix) = true
+
+
+####
+# eigvals
+#####
+
+eigvals(A::Symmetric{<:Real,<:BBBArrowheadMatrix}) = eigvals(Symmetric(BandedMatrix(parent(A))))
+eigvals(A::Symmetric{<:Real,<:BBBArrowheadMatrix}, B::Symmetric{<:Real,<:BBBArrowheadMatrix}) = eigvals(Symmetric(BandedMatrix(parent(A))), Symmetric(BandedMatrix(parent(B))))
+
+end #module

test/test_arrowhead.jl

@@ -1,4 +1,4 @@
-using PiecewiseOrthogonalPolynomials, FillArrays, BandedMatrices, MatrixFactorizations, BlockBandedMatrices, Base64, ClassicalOrthogonalPolynomials, Test
+using PiecewiseOrthogonalPolynomials, FillArrays, BandedMatrices, MatrixFactorizations, BlockBandedMatrices, BandedMatrices, Base64, ClassicalOrthogonalPolynomials, LinearAlgebra, Test
 using PiecewiseOrthogonalPolynomials: BBBArrowheadMatrix
 using InfiniteArrays, BlockArrays
 using BandedMatrices: _BandedMatrix
@@ -144,4 +144,34 @@ import Base: oneto, OneTo
         @test U ≈ reversecholesky(Matrix(Symmetric(A))).U
         @test U*U' ≈ Symmetric(A)
     end
+
+    @testset "banded" begin
+        n = 5; p = 5;
+        A = BBBArrowheadMatrix(BandedMatrix(0 => 1:n, 1 => 1:n-1, -1 => 1:n-1),
+                                ntuple(_ -> BandedMatrix((0 => randn(n-1), -1 => randn(n-1)), (n,n-1)), 2),
+                                ntuple(_ -> BandedMatrix((0 => randn(n), 1 => randn(n-1)), (n-1,n)), 3),
+                            fill(BandedMatrix((0 => randn(p) .+ 10, 2 => randn(p-2), -1=> randn(p-1)), (p, p)), n-1))
+
+
+        @test blockbandwidths(A) == (3,2)
+        @test subblockbandwidths(A) == (1,1)
+        @test BandedMatrices.isbanded(A)
+        @test bandwidths(A) == (13,9)
+        @test BandedMatrix(A) == A
+    end
+
+    @testset "eigvals" begin
+        r = range(-1,1; length=4)
+        C = ContinuousPolynomial{1}(r)
+        KR = Block.(OneTo(5))
+        Δ = weaklaplacian(C)[KR,KR]
+        M = grammatrix(C)[KR,KR]
+        ΔB = BandedMatrix(Δ)
+        MB = BandedMatrix(M)
+        @test bandwidths(Δ) == bandwidths(ΔB) == (1,1)
+        @test bandwidths(M) == bandwidths(MB) == (7,7)
+        @test eigvals(Symmetric(Δ)) == eigvals(Symmetric(ΔB))
+        @test eigvals(Symmetric(M)) == eigvals(Symmetric(MB))
+        @test eigvals(Symmetric(Δ), Symmetric(M)) == eigvals(Symmetric(ΔB), Symmetric(MB))
+    end
 end