From 7ca43ceb57c1c061f9c2fea635a04cd7ced7254b Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Fri, 15 Mar 2024 10:43:23 -0500 Subject: [PATCH] add Symmetric(Banded)ToeplitzPlusHankel --- Project.toml | 4 +- src/FastTransforms.jl | 8 +- src/SymmetricToeplitzPlusHankel.jl | 135 +++++++++++++++++++++++ test/runtests.jl | 3 +- test/symmetrictoeplitzplushankeltests.jl | 39 +++++++ 5 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 src/SymmetricToeplitzPlusHankel.jl create mode 100644 test/symmetrictoeplitzplushankeltests.jl diff --git a/Project.toml b/Project.toml index 2f5b3613..8bba0be1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.16" +version = "0.16.0" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FastTransforms_jll = "34b6f7d7-08f9-5794-9e10-3819e4c7e49a" @@ -17,6 +18,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" [compat] AbstractFFTs = "1.0" +BandedMatrices = "1.5" FFTW = "1.7" FastGaussQuadrature = "0.4, 0.5, 1" FastTransforms_jll = "0.6.2" diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl index e9bf6fac..82c34990 100644 --- a/src/FastTransforms.jl +++ b/src/FastTransforms.jl @@ -1,6 +1,6 @@ module FastTransforms -using FastGaussQuadrature, FillArrays, LinearAlgebra, +using BandedMatrices, FastGaussQuadrature, FillArrays, LinearAlgebra, Reexport, SpecialFunctions, ToeplitzMatrices @reexport using AbstractFFTs @@ -19,6 +19,8 @@ import AbstractFFTs: Plan, ScaledPlan, fftshift, ifftshift, rfft_output_size, brfft_output_size, normalization +import BandedMatrices: bandwidths + import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!, plan_dct, plan_idct, fftwNumber @@ -26,7 +28,7 @@ import FastGaussQuadrature: unweightedgausshermite import FillArrays: AbstractFill, getindex_value -import LinearAlgebra: mul!, lmul!, ldiv! +import LinearAlgebra: mul!, lmul!, ldiv!, cholesky import GenericFFT: interlace # imported in downstream packages @@ -112,6 +114,8 @@ include("specialfunctions.jl") include("toeplitzplans.jl") include("toeplitzhankel.jl") +include("SymmetricToeplitzPlusHankel.jl") + # following use libfasttransforms by default for f in (:jac2jac, :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb, diff --git a/src/SymmetricToeplitzPlusHankel.jl b/src/SymmetricToeplitzPlusHankel.jl new file mode 100644 index 00000000..a9bc4013 --- /dev/null +++ b/src/SymmetricToeplitzPlusHankel.jl @@ -0,0 +1,135 @@ +struct SymmetricToeplitzPlusHankel{T} <: AbstractMatrix{T} + v::Vector{T} + n::Int +end + +function SymmetricToeplitzPlusHankel(v::Vector{T}) where T + n = (length(v)+1)÷2 + SymmetricToeplitzPlusHankel{T}(v, n) +end + +size(A::SymmetricToeplitzPlusHankel{T}) where T = (A.n, A.n) +getindex(A::SymmetricToeplitzPlusHankel{T}, i::Integer, j::Integer) where T = A.v[abs(i-j)+1] + A.v[i+j-1] + +struct SymmetricBandedToeplitzPlusHankel{T} <: BandedMatrices.AbstractBandedMatrix{T} + v::Vector{T} + n::Int + b::Int +end + +function SymmetricBandedToeplitzPlusHankel(v::Vector{T}, n::Integer) where T + SymmetricBandedToeplitzPlusHankel{T}(v, n, length(v)-1) +end + +size(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.n, A.n) +function getindex(A::SymmetricBandedToeplitzPlusHankel{T}, i::Integer, j::Integer) where T + v = A.v + if abs(i-j) < length(v) + if i+j-1 ≤ length(v) + v[abs(i-j)+1] + v[i+j-1] + else + v[abs(i-j)+1] + end + else + zero(T) + end +end +bandwidths(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.b, A.b) + +# +# Jac*W-W*Jac' = G*J*G' +# This returns G and J, where J = [0 I; -I 0], respecting the skew-symmetry of the right-hand side. +# +function compute_skew_generators(A::SymmetricToeplitzPlusHankel{T}) where T + v = A.v + n = size(A, 1) + J = [zero(T) one(T); -one(T) zero(T)] + G = zeros(T, n, 2) + G[n, 1] = one(T) + u2 = reverse(v[2:n+1]) + u2[1:n-1] .+= v[n+1:2n-1] + G[:, 2] .= -u2 + G, J +end + +function cholesky(A::SymmetricToeplitzPlusHankel{T}) where T + n = size(A, 1) + G, J = compute_skew_generators(A) + L = zeros(T, n, n) + r = A[:, 1] + r2 = zeros(T, n) + l = zeros(T, n) + v = zeros(T, n) + col1 = zeros(T, n) + STPHcholesky!(L, G, r, r2, l, v, col1, n) + return Cholesky(L, 'L', 0) +end + +function STPHcholesky!(L::Matrix{T}, G, r, r2, l, v, col1, n) where T + @inbounds @simd for k in 1:n-1 + x = sqrt(r[1]) + for j in 1:n-k+1 + L[j+k-1, k] = l[j] = r[j]/x + end + for j in 1:n-k+1 + v[j] = G[j, 1]*G[1,2]-G[j,2]*G[1,1] + end + F12 = k == 1 ? T(2) : T(1) + r2[1] = (r[2] - v[1])/F12 + for j in 2:n-k + r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j] - v[j])/F12 + end + r2[n-k+1] = (r[n-k] + r[1]*col1[n-k+1] - col1[1]*r[n-k+1] - v[n-k+1])/F12 + cst = r[2]/x + for j in 1:n-k + r[j] = r2[j+1] - cst*l[j+1] + end + for j in 1:n-k + col1[j] = -F12/x*l[j+1] + end + c1 = G[1, 1] + c2 = G[1, 2] + for j in 1:n-k + G[j, 1] = G[j+1, 1] - l[j+1]*c1/x + G[j, 2] = G[j+1, 2] - l[j+1]*c2/x + end + end + L[n, n] = sqrt(r[1]) +end + +function cholesky(A::SymmetricBandedToeplitzPlusHankel{T}) where T + n = size(A, 1) + b = A.b + R = BandedMatrix{T}(undef, (n, n), (0, bandwidth(A, 2))) + r = A[1:b+2, 1] + r2 = zeros(T, b+3) + l = zeros(T, b+3) + col1 = zeros(T, b+2) + SBTPHcholesky!(R, r, r2, l, col1, n, b) + return Cholesky(R, 'U', 0) +end + +function SBTPHcholesky!(R::BandedMatrix{T}, r, r2, l, col1, n, b) where T + @inbounds @simd for k in 1:n + x = sqrt(r[1]) + for j in 1:b+1 + l[j] = r[j]/x + end + for j in 1:min(n-k+1, b+1) + R[k, j+k-1] = l[j] + end + F12 = k == 1 ? T(2) : T(1) + r2[1] = r[2]/F12 + for j in 2:b+1 + r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j])/F12 + end + r2[b+2] = (r[b+1] + r[1]*col1[b+2] - col1[1]*r[b+2])/F12 + cst = r[2]/x + for j in 1:b+2 + r[j] = r2[j+1] - cst*l[j+1] + end + for j in 1:b+2 + col1[j] = -F12/x*l[j+1] + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 1efb21e9..de16f36b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,4 +10,5 @@ include("gaunttests.jl") include("hermitetests.jl") include("clenshawtests.jl") include("toeplitzplanstests.jl") -include("toeplitzhankeltests.jl") \ No newline at end of file +include("toeplitzhankeltests.jl") +include("symmetrictoeplitzplushankeltests.jl") diff --git a/test/symmetrictoeplitzplushankeltests.jl b/test/symmetrictoeplitzplushankeltests.jl new file mode 100644 index 00000000..c26c38d3 --- /dev/null +++ b/test/symmetrictoeplitzplushankeltests.jl @@ -0,0 +1,39 @@ +using BandedMatrices, FastTransforms, LinearAlgebra, ToeplitzMatrices, Test + +import FastTransforms: SymmetricToeplitzPlusHankel, SymmetricBandedToeplitzPlusHankel + +@testset "SymmetricToeplitzPlusHankel" begin + n = 128 + for T in (Float32, Float64, BigFloat) + μ = -FastTransforms.chebyshevlogmoments1(T, 2n-1) + μ[1] += 1 + W = SymmetricToeplitzPlusHankel(μ/2) + SMW = Symmetric(Matrix(W)) + @test W ≈ SymmetricToeplitz(μ[1:(length(μ)+1)÷2]/2) + Hankel(μ/2) + L = cholesky(W).L + R = cholesky(SMW).U + @test L*L' ≈ W + @test L' ≈ R + end +end + +@testset "SymmetricBandedToeplitzPlusHankel" begin + n = 1024 + for T in (Float32, Float64) + μ = T[1.875; 0.00390625; 0.5; 0.0009765625; 0.0625] + W = SymmetricBandedToeplitzPlusHankel(μ/2, n) + SBW = Symmetric(BandedMatrix(W)) + W1 = SymmetricToeplitzPlusHankel([μ/2; zeros(2n-1-length(μ))]) + SMW = Symmetric(Matrix(W)) + U = cholesky(SMW).U + L = cholesky(W1).L + UB = cholesky(SBW).U + R = cholesky(W).U + @test L*L' ≈ W + @test UB'UB ≈ W + @test R'R ≈ W + @test UB ≈ U + @test L' ≈ U + @test R ≈ U + end +end