Introduce ArrowheadMatrix (#32)

* start arrowhead

* Work on LeftArrowhead

* change order for LeftArrowHead

* start symarrowhead

* Forget SymArrowhead

* more general arrowhead

* Update Project.toml

* LeftArrowhead -> Arrowhrad

* update tests

* Start reversecholesky for ArrowheadMatrix

* B blocks updated

* reversecholesky! works!

* tests bandwidths

* increase cov

* Use Arrowhead for building operators

* start explicit weaklaplacian

* Add mass matrix

* grammatrix

* Fast cholesky works!

* Update test_arrowhead.jl

* Getting tests to pass

* Update test_arrowhead.jl

* Update Project.toml

* Start UpperTriangular

* adjtrans for Arrowhead

* add algebra tests

Author: dlfivefifty

.github/workflows/ci.yml

@@ -10,8 +10,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.7'
-          - '1'
+          - '1.9'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,9 +1,11 @@
 name = "PiecewiseOrthogonalPolynomials"
 uuid = "4461d12d-4663-4550-8580-cb764c85e20f"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.0.3"
+version = "0.1.0"
 
 [deps]
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
@@ -14,25 +16,29 @@ InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 
 [compat]
 BlockArrays = "0.16.25"
-BlockBandedMatrices = "0.11.10, 0.12"
-ClassicalOrthogonalPolynomials = "0.7, 0.8, 0.9"
-ContinuumArrays = "0.12"
-FillArrays = "0.13, 1.0"
+BlockBandedMatrices = "0.12"
+ClassicalOrthogonalPolynomials = "0.11"
+ContinuumArrays = "0.14"
+FillArrays = "1.0"
 InfiniteArrays = "0.12.6"
 InfiniteLinearAlgebra = "0.6.16"
-LazyArrays = "0.22.9, 1"
-LazyBandedMatrices = "0.8.7"
-QuasiArrays = "0.9"
-julia = "1.7"
+LazyArrays = "1.0"
+LazyBandedMatrices = "0.8.13"
+MatrixFactorizations = "2.0"
+QuasiArrays = "0.11"
+julia = "1.9"
 
 [extras]
+Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "StaticArrays", "ForwardDiff"]
+test = ["Base64", "ForwardDiff", "Random", "StaticArrays", "Test"]

src/PiecewiseOrthogonalPolynomials.jl

@@ -1,15 +1,19 @@
 module PiecewiseOrthogonalPolynomials
-using ClassicalOrthogonalPolynomials, LinearAlgebra, BlockArrays, BlockBandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays
+using ClassicalOrthogonalPolynomials, LinearAlgebra, BlockArrays, BlockBandedMatrices, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays, MatrixFactorizations, ArrayLayouts
 
+import ArrayLayouts: sublayout, sub_materialize, symmetriclayout, transposelayout, SymmetricLayout, HermitianLayout, TriangularLayout, layout_getindex, materialize!, MatLdivVec, AbstractStridedLayout, triangulardata
+import BandedMatrices: _BandedMatrix
 import BlockArrays: BlockSlice, block, blockindex, blockvec
-import BlockBandedMatrices: _BandedBlockBandedMatrix
-import ClassicalOrthogonalPolynomials: grid, massmatrix, ldiv, pad, adaptivetransform_ldiv
-import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_grid_transform
+import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, subblockbandwidths, blockbandwidths, AbstractBandedBlockBandedLayout, layout_replace_in_print_matrix
+import ClassicalOrthogonalPolynomials: grid, ldiv, pad, adaptivetransform_ldiv, grammatrix
+import ContinuumArrays: @simplify, factorize, TransformFactorization, AbstractBasisLayout, MemoryLayout, layout_broadcasted, ExpansionLayout, basis, plan_grid_transform, grammatrix
 import LazyArrays: paddeddata
-import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec
-import Base: axes, getindex, ==, \, OneTo
+import LazyBandedMatrices: BlockBroadcastMatrix, BlockVec, BandedLazyLayouts, AbstractLazyBandedBlockBandedLayout, UpperOrLowerTriangular
+import Base: axes, getindex, +, -, *, /, ==, \, OneTo, oneto, replace_in_print_matrix, copy, diff, getproperty, adjoint, transpose, tail
+import LinearAlgebra: BlasInt
+import MatrixFactorizations: reversecholcopy
 
-export PiecewisePolynomial, ContinuousPolynomial, Derivative, Block
+export PiecewisePolynomial, ContinuousPolynomial, Derivative, Block, weaklaplacian, grammatrix
 
 abstract type AbstractPiecewisePolynomial{order,T,P<:AbstractVector} <: Basis{T} end
 
@@ -204,35 +208,83 @@ function \(P::ContinuousPolynomial{0}, C::ContinuousPolynomial{1})
     _BandedBlockBandedMatrix(dat, axes(P, 2), (1, 1), (0, 1))
 end
 
+function \(P::ContinuousPolynomial{0, <:Any, <:AbstractRange}, C::ContinuousPolynomial{1, <:Any, <:AbstractRange})
+    T = promote_type(eltype(P), eltype(C))
+    @assert P.points == C.points
+    v = (convert(T, 2):2:∞) ./ (3:2:∞)
+    N = length(P.points)
+    L = ArrowheadMatrix(_BandedMatrix(Ones{T}(2, N)/2, oneto(N-1), 0, 1),
+        (_BandedMatrix(Fill(v[1], 1, N-1), oneto(N-1), 0, 0),),
+        (_BandedMatrix(Vcat(Ones{T}(1, N)/2, -Ones{T}(1, N)/2), oneto(N-1), 0, 1),),
+        Fill(_BandedMatrix(Hcat(v, Zeros{T}(∞), -v)', axes(v,1), 1, 1), N-1))
+end
+
+
+
 ######
 # Gram matrix
 ######
 
-@simplify function *(Ac::QuasiAdjoint{<:Any,<:ContinuousPolynomial{0}}, B::ContinuousPolynomial{0})
-    A = Ac'
-    T = promote_type(eltype(A), eltype(B))
+function grammatrix(A::ContinuousPolynomial{0,T}) where T
+    r = A.points
+    N = length(r) - 1
+    hs = diff(r)
+    M = grammatrix(Legendre{T}())
+    ArrowheadMatrix{T}(Diagonal(Fill(hs[1], N)), (), (), [Diagonal(M.diag[2:end] * h/2) for h in hs])
+end
+
+function grammatrix(A::ContinuousPolynomial{0,T, <:AbstractRange}) where T
     r = A.points
-    @assert r == B.points
     N = length(r)
-    M = massmatrix(Legendre{T}())
+    M = grammatrix(Legendre{T}())
     Diagonal(mortar(Fill.((step(r) / 2) .* M.diag, N - 1)))
 end
 
+function grammatrix(C::ContinuousPolynomial{1, T, <:AbstractRange}) where T
+    r = C.points
+
+    N = length(r) - 1
+    h = step(r) # 2/N
+    a = ((convert(T,4):4:∞) .* (convert(T,-2):2:∞)) ./ ((1:2:∞) .* (3:2:∞) .* (-1:2:∞))
+    b = (((convert(T,2):2:∞) ./ (3:2:∞)).^2 .* (convert(T,2) ./ (1:2:∞) .+ convert(T,2) ./ (5:2:∞)))
+
+    a11 = LazyBandedMatrices.Bidiagonal(Vcat(h/3, Fill(2h/3, N-1), h/3), Fill(h/6, N), :U)
+    a21 = _BandedMatrix(Fill(h/3, 2, N), N+1, 1, 0)
+    a31 = _BandedMatrix(Vcat(Fill(-2h/15, 1, N), Fill(2h/15, 1, N)), N+1, 1, 0)
+
+    Symmetric(ArrowheadMatrix(a11, (a21, a31), (),
+                Fill(_BandedMatrix(Vcat((-h*a/2)',
+                Zeros(1,∞),
+                (h*b/2)'), ∞, 0, 2), N)))
+end
+
+
+function grammatrix(C::ContinuousPolynomial)
+    P = ContinuousPolynomial{0}(C)
+    L = P \ C
+    L' * grammatrix(P) * L
+end
 
 @simplify function *(Ac::QuasiAdjoint{<:Any,<:ContinuousPolynomial}, B::ContinuousPolynomial)
     A = Ac'
+    A == B && return grammatrix(A)
     P = ContinuousPolynomial{0}(A)
-    Q = ContinuousPolynomial{0}(B)
-    (P \ A)' * (P'Q) * (Q \ B)
+    (P \ A)' * grammatrix(P) * (P \ B)
+end
+
+@simplify function *(Ac::QuasiAdjoint{<:Any,<:ContinuousPolynomial{0}}, B::ContinuousPolynomial{0})
+    A = Ac'
+    @assert A == B
+    grammatrix(A)
 end
 
+
+
 #####
 # Derivative
 #####
 
-@simplify function *(D::Derivative, C::ContinuousPolynomial{1})
-    T = promote_type(eltype(D), eltype(C))
-
+function diff(C::ContinuousPolynomial{1,T}; dims=1) where T
     # Legendre() \ (D*Weighted(Jacobi(1,1)))
     r = C.points
     N = length(r)
@@ -242,7 +294,18 @@ end
     H = BlockBroadcastArray(hcat, z, v)
     M = BlockVcat(Hcat(Ones{T}(N) .* [zero(T); s] , -Ones{T}(N) .* [s; zero(T)] ), H)
     P = ContinuousPolynomial{0}(C)
-    P * _BandedBlockBandedMatrix(M', (axes(P, 2), axes(C, 2)), (0, 0), (0, 1))
+    ApplyQuasiMatrix(*, P, _BandedBlockBandedMatrix(M', (axes(P, 2), axes(C, 2)), (0, 0), (0, 1)))
+end
+
+function weaklaplacian(C::ContinuousPolynomial{1,T,<:AbstractRange}) where T
+    r = C.points
+    N = length(r)
+    s = step(r)
+    si = inv(s)
+    t1 = Vcat(-si, Fill(-2si, N-2), -si)
+    t2 = Fill(si, N-1)
+    Symmetric(ArrowheadMatrix(LazyBandedMatrices.Bidiagonal(t1, t2, :U), (), (),
+        Fill(Diagonal(convert(T, -16) .* (1:∞) .^ 2 ./ (s .* ((2:2:∞) .+ 1))), N-1)))
 end
 
 
@@ -268,5 +331,7 @@ function layout_broadcasted(::Tuple{ExpansionLayout{PiecewisePolynomialLayout{0}
     (a .* P) * (P \ C)
 end
 
+include("arrowhead.jl")
+
 
 end # module