Skip to content

Commit c790df0

Browse files
dlfivefiftydlfivefifty
authored
Improvements to ModalInterlace (#77)
* Improvements to ModalInterlace * Update ModalInterlace.jl * fix indenting
1 parent 097dc70 commit c790df0

File tree

6 files changed

+75
-53
lines changed

6 files changed

+75
-53
lines changed

‎examples/diskhelmholtz.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,11 @@ W = Weighted(Z)
66
xy = axes(Z,1); x,y = first.(xy),last.(xy)
77
Δ = Z \ (Laplacian(xy) * W)
88
S = Z \ W
9+
910
k = 2
1011
f = @.(cos(x*exp(y)))
1112
F = factorize(Δ + k^2 * S)
12-
c = (Z \ f)
13+
c = Z \ f
1314
F \ c
1415

1516
u = W * ((Δ + k^2 * S) \ (Z \ f))

‎src/ModalInterlace.jl

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
"""
2+
ModalInterlace
3+
"""
4+
struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
5+
ops
6+
MN::MMNN
7+
bandwidths::NTuple{2,Int}
8+
end
9+
10+
ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)
11+
ModalInterlace(ops::AbstractVector{<:AbstractMatrix{T}}, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ModalInterlace{T}(ops, MN, bandwidths)
12+
13+
axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
14+
15+
blockbandwidths(R::ModalInterlace) = R.bandwidths
16+
subblockbandwidths(::ModalInterlace) = (0,0)
17+
18+
19+
function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
20+
K,J = KJ.n
21+
dat = Matrix{T}(undef,1,J)
22+
l,u = blockbandwidths(R)
23+
if iseven(J-K) && -l ≤ J - K ≤ u
24+
sh = (J-K)÷2
25+
if isodd(K)
26+
k = K÷2+1
27+
dat[1,1] = R.ops[1][k,k+sh]
28+
end
29+
for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
30+
k = K÷2-m÷2+isodd(K)
31+
dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
32+
end
33+
else
34+
fill!(dat, zero(T))
35+
end
36+
_BandedMatrix(dat, K, 0, 0)
37+
end
38+
39+
getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
40+
41+
struct ModalInterlaceLayout <: AbstractBandedBlockBandedLayout end
42+
struct LazyModalInterlaceLayout <: AbstractLazyBandedBlockBandedLayout end
43+
44+
MemoryLayout(::Type{<:ModalInterlace}) = ModalInterlaceLayout()
45+
MemoryLayout(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyModalInterlaceLayout()
46+
sublayout(::Union{ModalInterlaceLayout,LazyModalInterlaceLayout}, ::Type{<:NTuple{2,BlockSlice{<:BlockOneTo}}}) = ModalInterlaceLayout()
47+
48+
49+
function sub_materialize(::ModalInterlaceLayout, V::AbstractMatrix{T}) where T
50+
kr,jr = parentindices(V)
51+
KR,JR = kr.block,jr.block
52+
M,N = Int(last(KR)), Int(last(JR))
53+
R = parent(V)
54+
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)
55+
end
56+
57+
# act like lazy array
58+
Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -10,12 +10,12 @@ import DomainSets: boundary
1010
import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle
1111
import ContinuumArrays: @simplify, Weight, grid, plotgrid, TransformFactorization, Expansion
1212

13-
import ArrayLayouts: MemoryLayout
13+
import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
1414
import BlockArrays: block, blockindex, BlockSlice, viewblock
1515
import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedMatrix, _BandedMatrix, blockbandwidths, subblockbandwidths
1616
import LinearAlgebra: factorize
1717
import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout
18-
import LazyBandedMatrices: LazyBandedBlockBandedLayout
18+
import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout
1919
import InfiniteArrays: InfiniteCardinal
2020

2121
import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted
@@ -26,7 +26,7 @@ export UnitTriangle, UnitDisk, JacobiTriangle, TriangleWeight, WeightedTriangle,
2626
MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial, Zernike, RadialCoordinate,
2727
zerniker, zernikez, Weighted, Block, ZernikeWeight
2828

29-
29+
include("ModalInterlace.jl")
3030
include("disk.jl")
3131
include("triangle.jl")
3232

‎src/disk.jl

Lines changed: 0 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -249,55 +249,6 @@ getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,
249249
end
250250

251251

252-
"""
253-
ModalInterlace
254-
"""
255-
struct ModalInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
256-
ops
257-
MN::MMNN
258-
bandwidths::NTuple{2,Int}
259-
end
260-
261-
ModalInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T =
262-
ModalInterlace{T,typeof(MN)}(ops, MN, bandwidths)
263-
264-
# act like lazy array
265-
MemoryLayout(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyBandedBlockBandedLayout()
266-
Base.BroadcastStyle(::Type{<:ModalInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()
267-
268-
axes(Z::ModalInterlace) = blockedrange.(oneto.(Z.MN))
269-
270-
blockbandwidths(R::ModalInterlace) = R.bandwidths
271-
subblockbandwidths(::ModalInterlace) = (0,0)
272-
273-
274-
function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
275-
K,J = KJ.n
276-
dat = Matrix{T}(undef,1,J)
277-
l,u = blockbandwidths(R)
278-
if iseven(J-K) && -l ≤ J - K ≤ u
279-
sh = (J-K)÷2
280-
if isodd(K)
281-
k = K÷2+1
282-
dat[1,1] = R.ops[1][k,k+sh]
283-
end
284-
for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
285-
k = K÷2-m÷2+isodd(K)
286-
dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
287-
end
288-
else
289-
fill!(dat, zero(T))
290-
end
291-
_BandedMatrix(dat, K, 0, 0)
292-
end
293-
294-
getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
295-
296-
function getindex(R::ModalInterlace{T}, KR::BlockOneTo, JR::BlockOneTo) where T
297-
M,N = Int(last(KR)), Int(last(JR))
298-
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)
299-
end
300-
301252
function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
302253
TV = promote_type(T,V)
303254
A.a == B.a && A.b == B.b && return Eye{TV}(∞)

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
using MultivariateOrthogonalPolynomials, Test
22

3+
include("test_modalinterlace.jl")
34
include("test_disk.jl")
45
include("test_triangle.jl")
56
# include("test_dirichlettriangle.jl")

‎test/test_modalinterlace.jl

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
using MultivariateOrthogonalPolynomials, ArrayLayouts, BandedMatrices, Test
2+
import MultivariateOrthogonalPolynomials: ModalInterlace, ModalInterlaceLayout
3+
4+
@testset "ModalInterlace" begin
5+
ops = [brand(2,3,1,2), brand(1,2,1,1), brand(1,2,1,2)]
6+
A = ModalInterlace(ops, (3,5), (2,4))
7+
@test MemoryLayout(A) isa ModalInterlaceLayout
8+
@test A[[1,4],[1,4,11]] == ops[1]
9+
@test A[[2],[2,7]] == ops[2]
10+
@test A[[5],[5,12]] == A[[6],[6,13]] == ops[3]
11+
end

0 commit comments

Comments
 (0)