Improvements to ModalInterlace

examples/diskhelmholtz.jl

@@ -6,10 +6,11 @@ W = Weighted(Z)
 xy = axes(Z,1); x,y = first.(xy),last.(xy)
 Δ = Z \ (Laplacian(xy) * W)
 S = Z \ W
+
 k = 2
 f = @.(cos(x*exp(y)))
 F = factorize(Δ + k^2 * S)
-c = (Z \ f)
+c = Z \ f
 F \ c
 
 u = W * ((Δ + k^2 * S) \ (Z \ f))

src/ModalInterlace.jl

@@ -0,0 +1,58 @@
+"""
+ModalInterlace
+"""
+struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
+    ops
+    MN::MMNN
+    bandwidths::NTuple{2,Int}
+end
+
+ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)
+ModalInterlace(ops::AbstractVector{<:AbstractMatrix{T}}, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T}(ops, MN, bandwidths)
+
+axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
+
+blockbandwidths(R::ModalInterlace) = R.bandwidths
+subblockbandwidths(::ModalInterlace) = (0,0)
+
+
+function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
+    K,J = KJ.n
+    dat = Matrix{T}(undef,1,J)
+    l,u = blockbandwidths(R)
+    if iseven(J-K) && -l ≤ J - K ≤ u
+        sh = (J-K)÷2
+        if isodd(K)
+            k = K÷2+1
+            dat[1,1] = R.ops[1][k,k+sh]
+        end
+        for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
+            k = K÷2-m÷2+isodd(K)
+            dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
+        end
+    else
+        fill!(dat, zero(T))
+    end
+    _BandedMatrix(dat, K, 0, 0)
+end
+
+getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
+
+struct ModalInterlaceLayout <: AbstractBandedBlockBandedLayout end
+struct LazyModalInterlaceLayout <: AbstractLazyBandedBlockBandedLayout end
+
+MemoryLayout(::Type{<:ModalInterlace}) = ModalInterlaceLayout()
+MemoryLayout(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyModalInterlaceLayout()
+sublayout(::Union{ModalInterlaceLayout,LazyModalInterlaceLayout}, ::Type{<:NTuple{2,BlockSlice{<:BlockOneTo}}}) = ModalInterlaceLayout()
+
+
+function sub_materialize(::ModalInterlaceLayout, V::AbstractMatrix{T}) where T
+    kr,jr = parentindices(V)
+    KR,JR = kr.block,jr.block
+    M,N = Int(last(KR)), Int(last(JR))
+    R = parent(V)
+    ModalInterlace{T}([R.ops[m][1:(M-m+2)÷2,1:(N-m+2)÷2] for m=1:min(N,M)], (M,N), R.bandwidths)
+end
+
+# act like lazy array
+Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()

src/MultivariateOrthogonalPolynomials.jl

@@ -10,12 +10,12 @@ import DomainSets: boundary
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
 import ContinuumArrays: @simplify, Weight, grid, plotgrid, TransformFactorization, Expansion
 
-import ArrayLayouts: MemoryLayout
+import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
 import BlockArrays: block, blockindex, BlockSlice, viewblock
 import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
 import LinearAlgebra: factorize
 import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout
-import LazyBandedMatrices: LazyBandedBlockBandedLayout
+import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout
 import InfiniteArrays: InfiniteCardinal
 
 import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted
@@ -26,7 +26,7 @@ export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle,
       MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
       zerniker, zernikez, Weighted, Block, ZernikeWeight
 
-
+include("ModalInterlace.jl")
 include("disk.jl")
 include("triangle.jl")
 

src/disk.jl

@@ -249,55 +249,6 @@ getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,
 end
 
 
-"""
-    ModalInterlace
-"""
-struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
-    ops
-    MN::MMNN
-    bandwidths::NTuple{2,Int}
-end
-
-ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = 
-    ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)
-
-# act like lazy array
-MemoryLayout(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyBandedBlockBandedLayout()
-Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()
-
-axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
-
-blockbandwidths(R::ModalInterlace) = R.bandwidths
-subblockbandwidths(::ModalInterlace) = (0,0)
-
-
-function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
-    K,J = KJ.n
-    dat = Matrix{T}(undef,1,J)
-    l,u = blockbandwidths(R)
-    if iseven(J-K) && -l ≤ J - K ≤ u
-        sh = (J-K)÷2
-        if isodd(K)
-            k = K÷2+1
-            dat[1,1] = R.ops[1][k,k+sh]
-        end
-        for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
-            k = K÷2-m÷2+isodd(K)
-            dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
-        end
-    else
-        fill!(dat, zero(T))
-    end
-    _BandedMatrix(dat, K, 0, 0)
-end
-
-getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
-
-function getindex(R::ModalInterlace{T}, KR::BlockOneTo, JR::BlockOneTo) where T
-    M,N = Int(last(KR)), Int(last(JR))
-    ModalInterlace{T}([R.ops[m][1:(M-m+2)÷2,1:(N-m+2)÷2] for m=1:min(N,M)], (M,N), R.bandwidths)
-end
-
 function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
     TV = promote_type(T,V)
     A.a == B.a && A.b == B.b && return Eye{TV}(∞)

test/runtests.jl

@@ -1,5 +1,6 @@
 using MultivariateOrthogonalPolynomials, Test
 
+include("test_modalinterlace.jl")
 include("test_disk.jl")
 include("test_triangle.jl")
 # include("test_dirichlettriangle.jl")

test/test_modalinterlace.jl

@@ -0,0 +1,11 @@
+using MultivariateOrthogonalPolynomials, ArrayLayouts, BandedMatrices, Test
+import MultivariateOrthogonalPolynomials: ModalInterlace, ModalInterlaceLayout
+
+@testset "ModalInterlace" begin
+    ops = [brand(2,3,1,2), brand(1,2,1,1), brand(1,2,1,2)]
+    A = ModalInterlace(ops, (3,5), (2,4))
+    @test MemoryLayout(A) isa ModalInterlaceLayout
+    @test A[[1,4],[1,4,11]] == ops[1]
+    @test A[[2],[2,7]] == ops[2]
+    @test A[[5],[5,12]] == A[[6],[6,13]] == ops[3]
+end