Start itransform

Author: dlfivefifty

examples/diskhelmholtz.jl

@@ -1,14 +1,30 @@
-using MultivariateOrthogonalPolynomials
+using MultivariateOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Plots
 
 WZ = Weighted(Zernike(1))
 xy = axes(WZ,1)
 x,y = first.(xy),last.(xy)
-u = WZ \ @. (exp(x*cos(y))*(1-x^2-y^2))
+u = Zernike(1) \ @. exp(x*cos(y))
 
-Δ = (Zernike(1) \ (Laplacian(xy) * WZ))
-C = Zernike(1) \ WZ
-N = 40; @time C[Block.(1:N), Block.(1:N)];
-k = 100
+N = 50
+KR = Block.(Base.OneTo(N))
+Δ = BandedBlockBandedMatrix((Zernike(1) \ (Laplacian(xy) * WZ))[KR,KR],(0,0),(0,0))
+C = (Zernike(1) \ WZ)[KR,KR]
+k = 5
 L = Δ - k^2 * C
 
-MemoryLayout(L)
+v = (L \ u[KR])
+
+F = factorize(Zernike(1)[:,KR])
+F.plan \ v
+F |> typeof |> fieldnames
+
+grid(WZ[:,KR])
+
+F |>typeof |> fieldnames
+F.F * v
+
+m = DiskTrav(v).matrix
+
+plan_disk2cxf(m, 0, 0) * m
+
+

src/disk.jl

@@ -8,15 +8,39 @@
 struct DiskTrav{T, AA<:AbstractMatrix{T}} <: AbstractBlockVector{T}
     matrix::AA
     function DiskTrav{T, AA}(matrix::AA) where {T,AA<:AbstractMatrix{T}}
-        n,m = size(matrix)
-        isodd(m) && n == m ÷ 4 + 1 || throw(ArgumentError("size must match"))
+        m,n = size(matrix)
+        isodd(n) && m == n ÷ 4 + 1 || throw(ArgumentError("size must match"))
         new{T,AA}(matrix)
     end
 end
 
 DiskTrav{T}(matrix::AbstractMatrix{T}) where T = DiskTrav{T,typeof(matrix)}(matrix)
 DiskTrav(matrix::AbstractMatrix{T}) where T = DiskTrav{T}(matrix)
 
+function DiskTrav(v::AbstractVector{T}) where T
+    N = blocksize(v,1)
+    m = N ÷ 2 + 1
+    n = 4(m-1) + 1
+    mat = zeros(T, m, n)
+    for K in blockaxes(v,1)
+        K̃ = Int(K)
+        w = v[K]
+        if isodd(K̃)
+            mat[K̃÷2 + 1,1] = w[1]
+            for j = 2:2:K̃-1
+                mat[K̃÷2-j÷2+1,2(j-1)+2] = w[j]
+                mat[K̃÷2-j÷2+1,2(j-1)+3] = w[j+1]
+            end
+        else
+            for j = 1:2:K̃
+                mat[K̃÷2-j÷2,2(j-1)+2] = w[j]
+                mat[K̃÷2-j÷2,2(j-1)+3] = w[j+1]
+            end
+        end
+    end
+    DiskTrav(mat)
+end
+
 axes(A::DiskTrav) = (blockedrange(oneto(div(size(A.matrix,2),2,RoundUp))),)
 
 function getindex(A::DiskTrav, K::Block{1})
@@ -165,6 +189,7 @@ function ZernikeTransform{T}(N::Int, a::Number, b::Number) where T<:Real
     ZernikeTransform{T}(N, plan_disk2cxf(T, Ñ, a, b), plan_disk_analysis(T, Ñ, 4Ñ-3))
 end
 *(P::ZernikeTransform{T}, f::Matrix{T}) where T = DiskTrav(P.disk2cxf \ (P.analysis * f))[Block.(1:P.N)]
+\(P::ZernikeTransform, f::AbstractVector) = P.analysis \ (P.disk2cxf * DiskTrav(f).matrix)
 
 factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))
 

test/test_disk.jl

@@ -46,6 +46,15 @@ import ClassicalOrthogonalPolynomials: HalfWeighted
         @test_throws ArgumentError DiskTrav([1 2 3 4])
         @test_throws ArgumentError DiskTrav([1 2 3; 4 5 6])
         @test_throws ArgumentError DiskTrav([1 2 3 4; 5 6 7 8])
+
+        for N = 1:10
+            v = PseudoBlockArray(1:sum(1:N),1:N)
+            if iseven(N)
+                @test DiskTrav(v) == [v; zeros(N+1)]
+            else
+                @test DiskTrav(v) == v
+            end
+        end
     end
 
     @testset "Transform" begin