diff --git a/Project.toml b/Project.toml index aa473c0..640b8df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FourierSeriesEvaluators" uuid = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5" authors = ["lxvm "] -version = "1.1.1-DEV" +version = "1.1.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/ext/FourierSeriesEvaluatorsStaticArraysExt.jl b/ext/FourierSeriesEvaluatorsStaticArraysExt.jl index 9c70729..d9ea948 100644 --- a/ext/FourierSeriesEvaluatorsStaticArraysExt.jl +++ b/ext/FourierSeriesEvaluatorsStaticArraysExt.jl @@ -1,7 +1,24 @@ module FourierSeriesEvaluatorsStaticArraysExt import FourierSeriesEvaluators: _hermitianpart using StaticArrays: SHermitianCompact, StaticMatrix +using LinearAlgebra: Hermitian, checksquare -_hermitianpart(A::StaticMatrix) = SHermitianCompact(A) - +_hermitianpart(A::StaticMatrix) = SHermitianCompact{checksquare(A)}(Hermitian(A, 'L')) +#= +# the code above runs as efficiently as this manually written version +function _hermitianpart(A::StaticMatrix) + n = checksquare(A) + M = ntuple(Val(n)) do j + o = (j-1)*(n+1) + ntuple(Val(n-j+1)) do i + a = A[i+o] + a isa Complex || return a # avoid type instability for real coefficients + i == 1 ? complex(real(a)) : a + end + end + SHermitianCompact(SVector(_tcat(M...))) +end +_tcat(a, b) = (a..., b...) +_tcat(a,b,c...) = _tcat(_tcat(a, b), c...) +=# end diff --git a/src/hermitian.jl b/src/hermitian.jl index 5cc18e2..c5b36ae 100644 --- a/src/hermitian.jl +++ b/src/hermitian.jl @@ -158,7 +158,7 @@ Inplace series are not currently supported function HermitianFourierSeries(f::FourierSeries) isinplace(f) && throw(ArgumentError("inplace series are not supported")) _check_coeffs_hermsym(f.c, f.o) || throw(ArgumentError("Fourier coefficients do not have Hermitian symmetry")) - _check_coeffs_centered(f.c, f.o) + _check_coeffs_centered(f.c, f.o) || throw(ArgumentError("Fourier coefficients are not centered")) ax = axes(f.c) _ax = ax[ndims(f.c)-length(f.o)+2:ndims(f.c)] ax_ = ax[1:ndims(f.c)-length(f.o)] @@ -209,4 +209,4 @@ frequency(s::HermitianFourierSeries) = s.f function nextderivative(s::HermitianFourierSeries, dim) a = raise_multiplier(s.a, dim) return HermitianFourierSeries{ndims(s),isinplace(s),typeof(s.c),typeof(a),eltype(s.t),eltype(s.f)}(s.c, a, s.t, s.f) -end \ No newline at end of file +end diff --git a/test/hermitian.jl b/test/hermitian.jl index 21dd121..4e4dc16 100644 --- a/test/hermitian.jl +++ b/test/hermitian.jl @@ -9,12 +9,22 @@ using FourierSeriesEvaluators: workspace_allocate, workspace_evaluate @testset "hermitian checks" begin n = 3 - c = rand(SMatrix{3,3,ComplexF64,9}, 2n+1,2n+1) + T = SMatrix{3,3,ComplexF64,9} + c = rand(T, 2n+1,2n+1) @test !_check_coeffs_hermsym(c, (nothing, nothing)) c = c + reverse(adjoint.(c)) @test _check_coeffs_hermsym(c, (nothing, nothing)) @test !_check_coeffs_centered(c, (-n, -n)) @test _check_coeffs_centered(c, (-n-1, -n-1)) + f = FourierSeries(rand(T, 2n+1, 2n+1); period=1.0, offset=-n-1) + @test_throws ArgumentError HermitianFourierSeries(f) + f = FourierSeries(c; period=1.0, offset=-n) + @test_throws ArgumentError HermitianFourierSeries(f) + # even-numbered + c = rand(SMatrix{3,3,ComplexF64,9}, 2n, 2n) + c = c + reverse(adjoint.(c)) + f = FourierSeries(c; period=1.0, offset=-n-1) + @test_throws ArgumentError HermitianFourierSeries(f) #= ninner=4 cc = Array{ComplexF64}(undef, ninner, 2n+1,2n+1) @@ -33,13 +43,12 @@ end n = 5 ntest = 10 T = Float64 - idx = Iterators.product(repeat([1:(2n+1)], d)...) + # idx = Iterators.product(repeat([1:(2n+1)], d)...) for a in [ - [rand(Complex{T}) for _ in idx], - # [rand(Complex{T}, 2, 2) for _ in idx], + rand(Complex{T}, ntuple(_->2n+1,d)...), # [rand(Complex{T}, 3, 3) for _ in idx], - [rand(SMatrix{2,2,Complex{T},4}) for _ in idx], - [rand(SMatrix{3,3,Complex{T},9}) for _ in idx], + rand(SMatrix{2,2,Complex{T},4}, ntuple(_->2n+1,d)...), + rand(SMatrix{3,3,Complex{T},9}, ntuple(_->2n+1,d)...), ] c = a + reverse(adjoint.(a)) f = FourierSeries(c; period=1.0, offset=-n-1) @@ -50,7 +59,7 @@ end for _ in 1:ntest x = rand(T, d) hfx = @inferred hf(x) - @test ishermitian(hfx) + @test ishermitian(hfx isa AbstractMatrix ? Matrix(hfx) : hfx) @test f(x) ≈ hfx if first(c) isa AbstractMatrix && d == 1 @test hfx ≈ Hermitian(f2(x)) @@ -64,12 +73,12 @@ end ntest = 10 ninner = 4 T = Float64 - idx = Iterators.product(ninner, repeat([1:(2n+1)], d)...) + # idx = Iterators.product(ninner, repeat([1:(2n+1)], d)...) for a in [ - [rand(Complex{T}) for _ in idx], + rand(Complex{T}, ntuple(_->2n+1,d)...), # [rand(Complex{T}, 3, 3) for _ in idx], - [rand(SMatrix{2,2,Complex{T},4}) for _ in idx], - [rand(SMatrix{3,3,Complex{T},9}) for _ in idx], + rand(SMatrix{2,2,Complex{T},4}, ntuple(_->2n+1,d)...), + rand(SMatrix{3,3,Complex{T},9}, ntuple(_->2n+1,d)...), ] c = Array{eltype(a)}(undef, ninner, 2n+1,2n+1) for i in axes(c,1) @@ -93,9 +102,8 @@ end d = 2 n = 5 ntest = 3 - idx = Iterators.product(repeat([1:(2n+1)], d)...) T = Float64 - c = [rand(SMatrix{2,2,Complex{T},4}) for _ in idx] + c = rand(SMatrix{2,2,Complex{T},4}, ntuple(_->2n+1,d)...) c = c + reverse(adjoint.(c)) f = FourierSeries(c; period=1.0, offset=-n-1) hf = HermitianFourierSeries(f) @@ -141,7 +149,7 @@ end d = 3 n = 3 T = SMatrix{2,2,ComplexF64,4} - C = rand(T, ntuple(_ -> n, d)...) + C = rand(T, ntuple(_ -> 2n+1, d)...) #= # inplace for nvar in 1:d-1 @@ -162,11 +170,11 @@ end =# C = C + reverse(adjoint.(C)) periods = rand(d) - f = FourierSeries(C, period=periods) + f = FourierSeries(C; period=periods, offset=-n-1) hf = HermitianFourierSeries(f) x = tuple(rand(d)...) for s in (hf, ManyFourierSeries(hf,hf,period=periods), JacobianSeries(hf)) ws = workspace_allocate(s, x) @test workspace_evaluate(ws, x) == s(x) end -end \ No newline at end of file +end