From 2993370b25c57fbddc17a4117d1d56469049ddbd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felipe=20A=2E=20V=2E=20de=20Bragan=C3=A7a=20Alves?= Date: Fri, 5 Jan 2018 00:58:50 -0500 Subject: [PATCH 01/12] Added rfft! and irfft! functionality through PaddedRFFTArray type. --- src/FFTW.jl | 1 + src/rfft!.jl | 216 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 + test/test_rfft!.jl | 66 ++++++++++++++ 4 files changed, 285 insertions(+) create mode 100644 src/rfft!.jl create mode 100644 test/test_rfft!.jl diff --git a/src/FFTW.jl b/src/FFTW.jl index 425e57f..a44b519 100644 --- a/src/FFTW.jl +++ b/src/FFTW.jl @@ -58,5 +58,6 @@ end include("fft.jl") include("dct.jl") +include("rfft!.jl") end # module diff --git a/src/rfft!.jl b/src/rfft!.jl new file mode 100644 index 0000000..dbe5320 --- /dev/null +++ b/src/rfft!.jl @@ -0,0 +1,216 @@ +import Base: IndexStyle, getindex, setindex!, eltype, \, similar, copy, real, read! + +#For compatibility between Julia v0.6 and v0.7 - begin +if VERSION >= v"0.7-" + using Base.@gc_preserve +else + macro gc_preserve(s::Symbol,ex::Expr) + return esc(ex) + end +end +#For compatibility between Julia v0.6 and v0.7 - end + +export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irfft! + + + +# As the time this code was written the new `ReinterpretArray` introduced in +# Julia v0.7 had major performace issues. Those issues were bypassed with the usage of the +# unsafe_wrap for the complex view of the data. As the name sugest, this is +# unsafe: whenever this complex view is called one must be careful to gc_preserve +# the "parent" PaddedRFFTArray. Therefore this view should not be "exported" to the +# user and the PaddedRFFTArray itself behaves like the complex view. +# Since it is not possible to not export a particular field of a exported type, +# a "hack" was used to name the unsafe field "#c", a fieldname that a non-advanced user +# will likely not be able to call. + +const c = Symbol("#c") +@eval struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} + r::SubArray{T,N,Array{T,N},NTuple{N,UnitRange{Int}},L} # Real view skipping padding + ($c)::Array{Complex{T},N} + + function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} + rrsize = size(rr) + fsize = rrsize[1] + iseven(fsize) || throw( + ArgumentError("First dimension of allocated array must have even number of elements")) + (nx == fsize-2 || nx == fsize-1) || throw( + ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) + fsize = fsize÷2 + csize = (fsize, rrsize[2:end]...) + if VERSION >= v"0.7-" + @gc_preserve rr c = unsafe_wrap(Array{Complex{T},N}, + reinterpret(Ptr{Complex{T}},pointer(rr)), + csize) + else + c = reinterpret(Complex{T}, rr, csize) + end + rsize = (nx,rrsize[2:end]...) + r = view(rr,(1:l for l in rsize)...) + return @gc_preserve rr new{T, N, N === 1 ? true : false}(r,c) + end # function +end # struct + +@inline real(S::PaddedRFFTArray) = S.r + +@inline unsafe_complex_view(S::PaddedRFFTArray) = getfield(S,c) + +copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(parent(real(S))),size(real(S))[1]) + +similar(f::PaddedRFFTArray,::Type{T},dims::Tuple{Vararg{Int,N}}) where {T, N} = + PaddedRFFTArray{T}(dims) +similar(f::PaddedRFFTArray{T,N,L},dims::NTuple{N2,Int}) where {T,N,L,N2} = + PaddedRFFTArray{T}(dims) +similar(f::PaddedRFFTArray,::Type{T}) where {T} = + PaddedRFFTArray{T}(size(real(f))) +similar(f::PaddedRFFTArray{T,N}) where {T,N} = + PaddedRFFTArray{T,N}(similar(parent(real(f))), size(real(f))[1]) + +size(S::PaddedRFFTArray) = + @gc_preserve S size(unsafe_complex_view(S)) + +IndexStyle(::Type{T}) where {T<:PaddedRFFTArray} = + IndexLinear() + +Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray, i::Int) = + @gc_preserve S getindex(unsafe_complex_view(S),i) + +Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray{T,N}, I::Vararg{Int, N}) where {T,N} = + @gc_preserve S getindex(unsafe_complex_view(S),I...) + +Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray,v,i::Int) = + @gc_preserve S setindex!(unsafe_complex_view(S),v,i) + +Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray{T,N},v,I::Vararg{Int,N}) where {T,N} = + @gc_preserve S setindex!(unsafe_complex_view(S),v,I...) + + +PaddedRFFTArray(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} = PaddedRFFTArray{T,N}(rr,nx) + +function PaddedRFFTArray{T}(ndims::Vararg{Integer,N}) where {T,N} + fsize = (ndims[1]÷2 + 1)*2 + a = Array{T,N}((fsize, ndims[2:end]...)) + PaddedRFFTArray{T,N}(a, ndims[1]) +end + +PaddedRFFTArray{T}(ndims::NTuple{N,Integer}) where {T,N} = + PaddedRFFTArray{T}(ndims...) + +PaddedRFFTArray(ndims::Vararg{Integer,N}) where N = + PaddedRFFTArray{Float64}(ndims...) + +PaddedRFFTArray(ndims::NTuple{N,Integer}) where N = + PaddedRFFTArray{Float64}(ndims...) + +function PaddedRFFTArray{T}(a::AbstractArray{<:Real,N}) where {T<:fftwReal,N} + t = PaddedRFFTArray{T}(size(a)) + @inbounds copy!(t.r, a) + return t +end + +PaddedRFFTArray(a::AbstractArray{<:Real}) = PaddedRFFTArray{Float64}(a) + +function PaddedRFFTArray(stream, dims) + field = PaddedRFFTArray(dims) + return read!(stream,field) +end + +function PaddedRFFTArray{T}(stream, dims) where T + field = PaddedRFFTArray{T}(dims) + return read!(stream,field) +end + +function read!(file::AbstractString, field::PaddedRFFTArray) + open(file) do io + return read!(io,field) + end +end + +# Read a binary file of an unpaded array directly to a PaddedRFFT array, without the need +# of the creation of a intermediary Array. If the data is already padded then the user +# should just use PaddedRFFTArray{T}(read("file",unpaddeddim),nx) +function read!(stream::IO, field::PaddedRFFTArray{T,N,L}) where {T,N,L} + rr = parent(field.r) + dims = size(real(field)) + nx = dims[1] + nb = sizeof(T)*nx + npencils = prod(dims)÷nx + npad = iseven(nx) ? 2 : 1 + for i=0:(npencils-1) + unsafe_read(stream,Ref(rr,Int((nx+npad)*i+1)),nb) + end + return field +end + + +########################################################################################### +# Foward plans + +function plan_rfft!(X::PaddedRFFTArray{T,N}, region; + flags::Integer=ESTIMATE, + timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} + + (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) + if flags&ESTIMATE != 0 + @gc_preserve X p = + rFFTWPlan{T,FORWARD,true,N}(real(X), unsafe_complex_view(X), region, flags, timelimit) + else + x = similar(X) + @gc_preserve x p = + rFFTWPlan{T,FORWARD,true,N}(real(x), unsafe_complex_view(x), region, flags, timelimit) + end + return p +end + +plan_rfft!(f::PaddedRFFTArray;kws...) = plan_rfft!(f, 1:ndims(f); kws...) + +*(p::rFFTWPlan{T,FORWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = + (@gc_preserve f A_mul_B!(unsafe_complex_view(f), p, real(f)); f) + +rfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_rfft!(f, region) * f + +function \(p::rFFTWPlan{T,FORWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} + isdefined(p,:pinv) || (p.pinv = plan_irfft!(f,p.region)) + return p.pinv * f +end + + +########################################################################################## +# Inverse plans + +function plan_brfft!(X::PaddedRFFTArray{T,N}, region; + flags::Integer=PRESERVE_INPUT, + timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} + (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) + if flags&PRESERVE_INPUT != 0 + a = similar(X) + return @gc_preserve a rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(a), real(a), region, flags,timelimit) + else + return @gc_preserve X rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(X), real(X), region, flags,timelimit) + end +end + +plan_brfft!(f::PaddedRFFTArray;kws...) = plan_brfft!(f,1:ndims(f);kws...) + +*(p::rFFTWPlan{Complex{T},BACKWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = + (@gc_preserve f A_mul_B!(real(f), p, unsafe_complex_view(f)); real(f)) + +brfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_brfft!(f, region) * f + +function plan_irfft!(x::PaddedRFFTArray{T,N}, region; kws...) where {T,N} + ScaledPlan(plan_brfft!(x, region; kws...),normalization(T, size(real(x)), region)) +end + +plan_irfft!(f::PaddedRFFTArray;kws...) = plan_irfft!(f,1:ndims(f);kws...) + +*(p::ScaledPlan,f::PaddedRFFTArray) = begin + p.p * f + scale!(parent(real(f)), p.scale) + real(f) +end + +irfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_irfft!(f,region) * f + + + + diff --git a/test/runtests.jl b/test/runtests.jl index 2874b11..b0623a2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -457,3 +457,5 @@ if fftw_vendor() != :mkl @test psXdct![i] ≈ true_Xdct[i] end end # fftw_vendor() != :mkl + +include("test_rfft!.jl") \ No newline at end of file diff --git a/test/test_rfft!.jl b/test/test_rfft!.jl new file mode 100644 index 0000000..39a26be --- /dev/null +++ b/test/test_rfft!.jl @@ -0,0 +1,66 @@ +a = rand(Float64,(8,4,4)) +b = PaddedRFFTArray(a) +c = copy(b) + +@testset "PaddedRFFTArray creation" begin + @test a == real(b) + @test c == b + @test c.r == b.r + @test typeof(similar(b)) === typeof(b) + @test size(similar(b,Float32)) === size(b) + @test size(similar(b,Float32).r) === size(b.r) + @test size(similar(b,(4,4,4)).r) === (4,4,4) + @test size(similar(b,Float32,(4,4,4)).r) === (4,4,4) +end + +@testset "rfft! and irfft!" begin + @test rfft(a) ≈ rfft!(b) + @test a ≈ irfft!(b) + @test rfft(a,1:2) ≈ rfft!(b,1:2) + @test a ≈ irfft!(b,1:2) + @test rfft(a,(1,3)) ≈ rfft!(b,(1,3)) + @test a ≈ irfft!(b,(1,3)) + + p = plan_rfft!(c) + @test p*c ≈ rfft!(b) + @test p\c ≈ irfft!(b) + + a = rand(9,4,4) + b = PaddedRFFTArray(a) + @test a == real(b) + @test rfft(a) ≈ rfft!(b) + @test a ≈ irfft!(b) + @test rfft(a,1:2) ≈ rfft!(b,1:2) + @test a ≈ irfft!(b,1:2) + @test rfft(a,(1,3)) ≈ rfft!(b,(1,3)) + @test a ≈ irfft!(b,(1,3)) +end + +@testset "Read binary file to PaddedRFFTArray" begin + for s in ((8,4,4),(9,4,4),(8,),(9,)) + aa = rand(s) + f = Base.Filesystem.tempname() + write(f,aa) + @test aa == real(PaddedRFFTArray(f,s)) + aa = rand(Float32,s) + write(f,aa) + @test aa == real(PaddedRFFTArray{Float32}(f,s)) + end +end + +@testset "brfft!" begin + a = rand(4,4) + b = PaddedRFFTArray(a) + rfft!(b) + @test (brfft!(b) ./ 16) ≈ a +end + +@testset "FFTW MEASURE flag" begin + c = similar(b) + p = plan_rfft!(c,flags=FFTW.MEASURE) + p.pinv = plan_irfft!(c,flags=FFTW.MEASURE) + c .= b + @test c == b + @test p*c ≈ rfft!(b) + @test p\c ≈ irfft!(b) +end \ No newline at end of file From 76e5bceb32115cdf2af61a379f6b041116c8f0ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Felipe=20A=2E=20V=2E=20de=20Bragan=C3=A7a=20Alves?= Date: Fri, 5 Jan 2018 00:58:50 -0500 Subject: [PATCH 02/12] Added rfft! and irfft! functionality through PaddedRFFTArray type. --- src/FFTW.jl | 1 + src/rfft!.jl | 216 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 + test/test_rfft!.jl | 66 ++++++++++++++ 4 files changed, 285 insertions(+) create mode 100644 src/rfft!.jl create mode 100644 test/test_rfft!.jl diff --git a/src/FFTW.jl b/src/FFTW.jl index ac0c5cb..5801e04 100644 --- a/src/FFTW.jl +++ b/src/FFTW.jl @@ -55,5 +55,6 @@ end include("fft.jl") include("dct.jl") +include("rfft!.jl") end # module diff --git a/src/rfft!.jl b/src/rfft!.jl new file mode 100644 index 0000000..dbe5320 --- /dev/null +++ b/src/rfft!.jl @@ -0,0 +1,216 @@ +import Base: IndexStyle, getindex, setindex!, eltype, \, similar, copy, real, read! + +#For compatibility between Julia v0.6 and v0.7 - begin +if VERSION >= v"0.7-" + using Base.@gc_preserve +else + macro gc_preserve(s::Symbol,ex::Expr) + return esc(ex) + end +end +#For compatibility between Julia v0.6 and v0.7 - end + +export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irfft! + + + +# As the time this code was written the new `ReinterpretArray` introduced in +# Julia v0.7 had major performace issues. Those issues were bypassed with the usage of the +# unsafe_wrap for the complex view of the data. As the name sugest, this is +# unsafe: whenever this complex view is called one must be careful to gc_preserve +# the "parent" PaddedRFFTArray. Therefore this view should not be "exported" to the +# user and the PaddedRFFTArray itself behaves like the complex view. +# Since it is not possible to not export a particular field of a exported type, +# a "hack" was used to name the unsafe field "#c", a fieldname that a non-advanced user +# will likely not be able to call. + +const c = Symbol("#c") +@eval struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} + r::SubArray{T,N,Array{T,N},NTuple{N,UnitRange{Int}},L} # Real view skipping padding + ($c)::Array{Complex{T},N} + + function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} + rrsize = size(rr) + fsize = rrsize[1] + iseven(fsize) || throw( + ArgumentError("First dimension of allocated array must have even number of elements")) + (nx == fsize-2 || nx == fsize-1) || throw( + ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) + fsize = fsize÷2 + csize = (fsize, rrsize[2:end]...) + if VERSION >= v"0.7-" + @gc_preserve rr c = unsafe_wrap(Array{Complex{T},N}, + reinterpret(Ptr{Complex{T}},pointer(rr)), + csize) + else + c = reinterpret(Complex{T}, rr, csize) + end + rsize = (nx,rrsize[2:end]...) + r = view(rr,(1:l for l in rsize)...) + return @gc_preserve rr new{T, N, N === 1 ? true : false}(r,c) + end # function +end # struct + +@inline real(S::PaddedRFFTArray) = S.r + +@inline unsafe_complex_view(S::PaddedRFFTArray) = getfield(S,c) + +copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(parent(real(S))),size(real(S))[1]) + +similar(f::PaddedRFFTArray,::Type{T},dims::Tuple{Vararg{Int,N}}) where {T, N} = + PaddedRFFTArray{T}(dims) +similar(f::PaddedRFFTArray{T,N,L},dims::NTuple{N2,Int}) where {T,N,L,N2} = + PaddedRFFTArray{T}(dims) +similar(f::PaddedRFFTArray,::Type{T}) where {T} = + PaddedRFFTArray{T}(size(real(f))) +similar(f::PaddedRFFTArray{T,N}) where {T,N} = + PaddedRFFTArray{T,N}(similar(parent(real(f))), size(real(f))[1]) + +size(S::PaddedRFFTArray) = + @gc_preserve S size(unsafe_complex_view(S)) + +IndexStyle(::Type{T}) where {T<:PaddedRFFTArray} = + IndexLinear() + +Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray, i::Int) = + @gc_preserve S getindex(unsafe_complex_view(S),i) + +Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray{T,N}, I::Vararg{Int, N}) where {T,N} = + @gc_preserve S getindex(unsafe_complex_view(S),I...) + +Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray,v,i::Int) = + @gc_preserve S setindex!(unsafe_complex_view(S),v,i) + +Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray{T,N},v,I::Vararg{Int,N}) where {T,N} = + @gc_preserve S setindex!(unsafe_complex_view(S),v,I...) + + +PaddedRFFTArray(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} = PaddedRFFTArray{T,N}(rr,nx) + +function PaddedRFFTArray{T}(ndims::Vararg{Integer,N}) where {T,N} + fsize = (ndims[1]÷2 + 1)*2 + a = Array{T,N}((fsize, ndims[2:end]...)) + PaddedRFFTArray{T,N}(a, ndims[1]) +end + +PaddedRFFTArray{T}(ndims::NTuple{N,Integer}) where {T,N} = + PaddedRFFTArray{T}(ndims...) + +PaddedRFFTArray(ndims::Vararg{Integer,N}) where N = + PaddedRFFTArray{Float64}(ndims...) + +PaddedRFFTArray(ndims::NTuple{N,Integer}) where N = + PaddedRFFTArray{Float64}(ndims...) + +function PaddedRFFTArray{T}(a::AbstractArray{<:Real,N}) where {T<:fftwReal,N} + t = PaddedRFFTArray{T}(size(a)) + @inbounds copy!(t.r, a) + return t +end + +PaddedRFFTArray(a::AbstractArray{<:Real}) = PaddedRFFTArray{Float64}(a) + +function PaddedRFFTArray(stream, dims) + field = PaddedRFFTArray(dims) + return read!(stream,field) +end + +function PaddedRFFTArray{T}(stream, dims) where T + field = PaddedRFFTArray{T}(dims) + return read!(stream,field) +end + +function read!(file::AbstractString, field::PaddedRFFTArray) + open(file) do io + return read!(io,field) + end +end + +# Read a binary file of an unpaded array directly to a PaddedRFFT array, without the need +# of the creation of a intermediary Array. If the data is already padded then the user +# should just use PaddedRFFTArray{T}(read("file",unpaddeddim),nx) +function read!(stream::IO, field::PaddedRFFTArray{T,N,L}) where {T,N,L} + rr = parent(field.r) + dims = size(real(field)) + nx = dims[1] + nb = sizeof(T)*nx + npencils = prod(dims)÷nx + npad = iseven(nx) ? 2 : 1 + for i=0:(npencils-1) + unsafe_read(stream,Ref(rr,Int((nx+npad)*i+1)),nb) + end + return field +end + + +########################################################################################### +# Foward plans + +function plan_rfft!(X::PaddedRFFTArray{T,N}, region; + flags::Integer=ESTIMATE, + timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} + + (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) + if flags&ESTIMATE != 0 + @gc_preserve X p = + rFFTWPlan{T,FORWARD,true,N}(real(X), unsafe_complex_view(X), region, flags, timelimit) + else + x = similar(X) + @gc_preserve x p = + rFFTWPlan{T,FORWARD,true,N}(real(x), unsafe_complex_view(x), region, flags, timelimit) + end + return p +end + +plan_rfft!(f::PaddedRFFTArray;kws...) = plan_rfft!(f, 1:ndims(f); kws...) + +*(p::rFFTWPlan{T,FORWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = + (@gc_preserve f A_mul_B!(unsafe_complex_view(f), p, real(f)); f) + +rfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_rfft!(f, region) * f + +function \(p::rFFTWPlan{T,FORWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} + isdefined(p,:pinv) || (p.pinv = plan_irfft!(f,p.region)) + return p.pinv * f +end + + +########################################################################################## +# Inverse plans + +function plan_brfft!(X::PaddedRFFTArray{T,N}, region; + flags::Integer=PRESERVE_INPUT, + timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} + (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) + if flags&PRESERVE_INPUT != 0 + a = similar(X) + return @gc_preserve a rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(a), real(a), region, flags,timelimit) + else + return @gc_preserve X rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(X), real(X), region, flags,timelimit) + end +end + +plan_brfft!(f::PaddedRFFTArray;kws...) = plan_brfft!(f,1:ndims(f);kws...) + +*(p::rFFTWPlan{Complex{T},BACKWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = + (@gc_preserve f A_mul_B!(real(f), p, unsafe_complex_view(f)); real(f)) + +brfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_brfft!(f, region) * f + +function plan_irfft!(x::PaddedRFFTArray{T,N}, region; kws...) where {T,N} + ScaledPlan(plan_brfft!(x, region; kws...),normalization(T, size(real(x)), region)) +end + +plan_irfft!(f::PaddedRFFTArray;kws...) = plan_irfft!(f,1:ndims(f);kws...) + +*(p::ScaledPlan,f::PaddedRFFTArray) = begin + p.p * f + scale!(parent(real(f)), p.scale) + real(f) +end + +irfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_irfft!(f,region) * f + + + + diff --git a/test/runtests.jl b/test/runtests.jl index 40fa645..884c605 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -485,3 +485,5 @@ if fftw_vendor() != :mkl @test psXdct![i] ≈ true_Xdct[i] end end # fftw_vendor() != :mkl + +include("test_rfft!.jl") \ No newline at end of file diff --git a/test/test_rfft!.jl b/test/test_rfft!.jl new file mode 100644 index 0000000..39a26be --- /dev/null +++ b/test/test_rfft!.jl @@ -0,0 +1,66 @@ +a = rand(Float64,(8,4,4)) +b = PaddedRFFTArray(a) +c = copy(b) + +@testset "PaddedRFFTArray creation" begin + @test a == real(b) + @test c == b + @test c.r == b.r + @test typeof(similar(b)) === typeof(b) + @test size(similar(b,Float32)) === size(b) + @test size(similar(b,Float32).r) === size(b.r) + @test size(similar(b,(4,4,4)).r) === (4,4,4) + @test size(similar(b,Float32,(4,4,4)).r) === (4,4,4) +end + +@testset "rfft! and irfft!" begin + @test rfft(a) ≈ rfft!(b) + @test a ≈ irfft!(b) + @test rfft(a,1:2) ≈ rfft!(b,1:2) + @test a ≈ irfft!(b,1:2) + @test rfft(a,(1,3)) ≈ rfft!(b,(1,3)) + @test a ≈ irfft!(b,(1,3)) + + p = plan_rfft!(c) + @test p*c ≈ rfft!(b) + @test p\c ≈ irfft!(b) + + a = rand(9,4,4) + b = PaddedRFFTArray(a) + @test a == real(b) + @test rfft(a) ≈ rfft!(b) + @test a ≈ irfft!(b) + @test rfft(a,1:2) ≈ rfft!(b,1:2) + @test a ≈ irfft!(b,1:2) + @test rfft(a,(1,3)) ≈ rfft!(b,(1,3)) + @test a ≈ irfft!(b,(1,3)) +end + +@testset "Read binary file to PaddedRFFTArray" begin + for s in ((8,4,4),(9,4,4),(8,),(9,)) + aa = rand(s) + f = Base.Filesystem.tempname() + write(f,aa) + @test aa == real(PaddedRFFTArray(f,s)) + aa = rand(Float32,s) + write(f,aa) + @test aa == real(PaddedRFFTArray{Float32}(f,s)) + end +end + +@testset "brfft!" begin + a = rand(4,4) + b = PaddedRFFTArray(a) + rfft!(b) + @test (brfft!(b) ./ 16) ≈ a +end + +@testset "FFTW MEASURE flag" begin + c = similar(b) + p = plan_rfft!(c,flags=FFTW.MEASURE) + p.pinv = plan_irfft!(c,flags=FFTW.MEASURE) + c .= b + @test c == b + @test p*c ≈ rfft!(b) + @test p\c ≈ irfft!(b) +end \ No newline at end of file From 9c8776138cdedf1761d4f426bcf958c048b63e49 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Thu, 29 Mar 2018 15:32:24 -0400 Subject: [PATCH 03/12] Removed unsafe complex view. Using custom getindex and setindex! --- src/rfft!.jl | 132 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 78 insertions(+), 54 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index dbe5320..9655b6d 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -1,33 +1,18 @@ import Base: IndexStyle, getindex, setindex!, eltype, \, similar, copy, real, read! -#For compatibility between Julia v0.6 and v0.7 - begin -if VERSION >= v"0.7-" - using Base.@gc_preserve -else - macro gc_preserve(s::Symbol,ex::Expr) - return esc(ex) - end -end -#For compatibility between Julia v0.6 and v0.7 - end - export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irfft! # As the time this code was written the new `ReinterpretArray` introduced in # Julia v0.7 had major performace issues. Those issues were bypassed with the usage of the -# unsafe_wrap for the complex view of the data. As the name sugest, this is -# unsafe: whenever this complex view is called one must be careful to gc_preserve -# the "parent" PaddedRFFTArray. Therefore this view should not be "exported" to the -# user and the PaddedRFFTArray itself behaves like the complex view. -# Since it is not possible to not export a particular field of a exported type, -# a "hack" was used to name the unsafe field "#c", a fieldname that a non-advanced user -# will likely not be able to call. - -const c = Symbol("#c") -@eval struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} +# custom getindex and setindex! below. Hopefully, once the performance issues with ReinterpretArray +# are solved we can just index the reinterpret array directly. + +struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} + data::Array{T,N} r::SubArray{T,N,Array{T,N},NTuple{N,UnitRange{Int}},L} # Real view skipping padding - ($c)::Array{Complex{T},N} + c::(@static VERSION >= v"0.7-" ? Base.ReinterpretArray{Complex{T},N,T,Array{T,N}} : Array{Complex{T},N}) function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} rrsize = size(rr) @@ -38,24 +23,22 @@ const c = Symbol("#c") ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) fsize = fsize÷2 csize = (fsize, rrsize[2:end]...) - if VERSION >= v"0.7-" - @gc_preserve rr c = unsafe_wrap(Array{Complex{T},N}, - reinterpret(Ptr{Complex{T}},pointer(rr)), - csize) - else - c = reinterpret(Complex{T}, rr, csize) - end + c = @static VERSION >= v"0.7-" ? + reinterpret(Complex{T}, rr) : + reinterpret(Complex{T}, rr, csize) rsize = (nx,rrsize[2:end]...) r = view(rr,(1:l for l in rsize)...) - return @gc_preserve rr new{T, N, N === 1 ? true : false}(r,c) + return new{T, N, N === 1 ? true : false}(rr,r,c) end # function end # struct @inline real(S::PaddedRFFTArray) = S.r -@inline unsafe_complex_view(S::PaddedRFFTArray) = getfield(S,c) +@inline complex_view(S::PaddedRFFTArray) = S.c -copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(parent(real(S))),size(real(S))[1]) +@inline data(S::PaddedRFFTArray) = S.data + +copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(data(S)),size(real(S),1)) similar(f::PaddedRFFTArray,::Type{T},dims::Tuple{Vararg{Int,N}}) where {T, N} = PaddedRFFTArray{T}(dims) @@ -64,32 +47,75 @@ similar(f::PaddedRFFTArray{T,N,L},dims::NTuple{N2,Int}) where {T,N,L,N2} = similar(f::PaddedRFFTArray,::Type{T}) where {T} = PaddedRFFTArray{T}(size(real(f))) similar(f::PaddedRFFTArray{T,N}) where {T,N} = - PaddedRFFTArray{T,N}(similar(parent(real(f))), size(real(f))[1]) + PaddedRFFTArray{T,N}(similar(data(f)), size(real(f),1)) size(S::PaddedRFFTArray) = - @gc_preserve S size(unsafe_complex_view(S)) + size(complex_view(S)) IndexStyle(::Type{T}) where {T<:PaddedRFFTArray} = IndexLinear() -Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray, i::Int) = - @gc_preserve S getindex(unsafe_complex_view(S),i) - -Base.@propagate_inbounds @inline getindex(S::PaddedRFFTArray{T,N}, I::Vararg{Int, N}) where {T,N} = - @gc_preserve S getindex(unsafe_complex_view(S),I...) +@inline function getindex(A::PaddedRFFTArray{T,N}, i2::Integer) where {T,N} + d = data(A) + i = 2i2 + @boundscheck checkbounds(d,i) + @inbounds begin + return Complex{T}(d[i-1],d[i]) + end +end -Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray,v,i::Int) = - @gc_preserve S setindex!(unsafe_complex_view(S),v,i) +@inline @generated function getindex(A::PaddedRFFTArray{T,N}, I2::Vararg{Integer,N}) where {T,N} + ip = :(2*I2[1]) + t = Expr(:tuple) + for i=2:N + push!(t.args,:(I2[$i])) + end + quote + d = data(A) + i = $ip + I = $t + @boundscheck checkbounds(d,i,I...) + @inbounds begin + return Complex{T}(d[i-1,I...],d[i,I...]) + end + end +end -Base.@propagate_inbounds @inline setindex!(S::PaddedRFFTArray{T,N},v,I::Vararg{Int,N}) where {T,N} = - @gc_preserve S setindex!(unsafe_complex_view(S),v,I...) +@inline function setindex!(A::PaddedRFFTArray{T,N},x, i2::Integer) where {T,N} + d = data(A) + i = 2i2 + @boundscheck checkbounds(d,i) + @inbounds begin + d[i-1] = real(x) + d[i] = imag(x) + end + A +end +@inline @generated function setindex!(A::PaddedRFFTArray{T,N}, x, I2::Vararg{Integer,N}) where {T,N} + ip = :(2*I2[1]) + t = Expr(:tuple) + for i=2:N + push!(t.args,:(I2[$i])) + end + quote + d = data(A) + i = $ip + I = $t + @boundscheck checkbounds(d,i,I...) + @inbounds begin + d[i-1,I...] = real(x) + d[i,I...] = imag(x) + end + A + end +end PaddedRFFTArray(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} = PaddedRFFTArray{T,N}(rr,nx) function PaddedRFFTArray{T}(ndims::Vararg{Integer,N}) where {T,N} fsize = (ndims[1]÷2 + 1)*2 - a = Array{T,N}((fsize, ndims[2:end]...)) + a = zeros(T,(fsize, ndims[2:end]...)) PaddedRFFTArray{T,N}(a, ndims[1]) end @@ -104,7 +130,7 @@ PaddedRFFTArray(ndims::NTuple{N,Integer}) where N = function PaddedRFFTArray{T}(a::AbstractArray{<:Real,N}) where {T<:fftwReal,N} t = PaddedRFFTArray{T}(size(a)) - @inbounds copy!(t.r, a) + @inbounds copyto!(t.r, a) return t end @@ -130,7 +156,7 @@ end # of the creation of a intermediary Array. If the data is already padded then the user # should just use PaddedRFFTArray{T}(read("file",unpaddeddim),nx) function read!(stream::IO, field::PaddedRFFTArray{T,N,L}) where {T,N,L} - rr = parent(field.r) + rr = data(field) dims = size(real(field)) nx = dims[1] nb = sizeof(T)*nx @@ -152,12 +178,10 @@ function plan_rfft!(X::PaddedRFFTArray{T,N}, region; (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) if flags&ESTIMATE != 0 - @gc_preserve X p = - rFFTWPlan{T,FORWARD,true,N}(real(X), unsafe_complex_view(X), region, flags, timelimit) + p = rFFTWPlan{T,FORWARD,true,N}(real(X), complex_view(X), region, flags, timelimit) else x = similar(X) - @gc_preserve x p = - rFFTWPlan{T,FORWARD,true,N}(real(x), unsafe_complex_view(x), region, flags, timelimit) + p = rFFTWPlan{T,FORWARD,true,N}(real(x), complex_view(x), region, flags, timelimit) end return p end @@ -165,7 +189,7 @@ end plan_rfft!(f::PaddedRFFTArray;kws...) = plan_rfft!(f, 1:ndims(f); kws...) *(p::rFFTWPlan{T,FORWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = - (@gc_preserve f A_mul_B!(unsafe_complex_view(f), p, real(f)); f) + (mul!(complex_view(f), p, real(f)); f) rfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_rfft!(f, region) * f @@ -184,16 +208,16 @@ function plan_brfft!(X::PaddedRFFTArray{T,N}, region; (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) if flags&PRESERVE_INPUT != 0 a = similar(X) - return @gc_preserve a rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(a), real(a), region, flags,timelimit) + return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(a), real(a), region, flags,timelimit) else - return @gc_preserve X rFFTWPlan{Complex{T},BACKWARD,true,N}(unsafe_complex_view(X), real(X), region, flags,timelimit) + return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(X), real(X), region, flags,timelimit) end end plan_brfft!(f::PaddedRFFTArray;kws...) = plan_brfft!(f,1:ndims(f);kws...) *(p::rFFTWPlan{Complex{T},BACKWARD,true,N},f::PaddedRFFTArray{T,N}) where {T<:fftwReal,N} = - (@gc_preserve f A_mul_B!(real(f), p, unsafe_complex_view(f)); real(f)) + (mul!(real(f), p, complex_view(f)); real(f)) brfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_brfft!(f, region) * f @@ -205,7 +229,7 @@ plan_irfft!(f::PaddedRFFTArray;kws...) = plan_irfft!(f,1:ndims(f);kws...) *(p::ScaledPlan,f::PaddedRFFTArray) = begin p.p * f - scale!(parent(real(f)), p.scale) + rmul!(data(f), p.scale) real(f) end From 87327dd9bd092cea600914239504d0fa4792176d Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Tue, 17 Jul 2018 16:24:40 -0400 Subject: [PATCH 04/12] Remove compatibility with v0.6. --- src/rfft!.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index 9655b6d..5261837 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -12,7 +12,7 @@ export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irf struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} data::Array{T,N} r::SubArray{T,N,Array{T,N},NTuple{N,UnitRange{Int}},L} # Real view skipping padding - c::(@static VERSION >= v"0.7-" ? Base.ReinterpretArray{Complex{T},N,T,Array{T,N}} : Array{Complex{T},N}) + c::Base.ReinterpretArray{Complex{T},N,T,Array{T,N}} function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} rrsize = size(rr) @@ -23,9 +23,7 @@ struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) fsize = fsize÷2 csize = (fsize, rrsize[2:end]...) - c = @static VERSION >= v"0.7-" ? - reinterpret(Complex{T}, rr) : - reinterpret(Complex{T}, rr, csize) + c = reinterpret(Complex{T}, rr) rsize = (nx,rrsize[2:end]...) r = view(rr,(1:l for l in rsize)...) return new{T, N, N === 1 ? true : false}(rr,r,c) @@ -234,7 +232,3 @@ plan_irfft!(f::PaddedRFFTArray;kws...) = plan_irfft!(f,1:ndims(f);kws...) end irfft!(f::PaddedRFFTArray, region=1:ndims(f)) = plan_irfft!(f,region) * f - - - - From e7e994aa55a37950b5b3f9e3c5adf50e131c4370 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Tue, 17 Jul 2018 16:39:36 -0400 Subject: [PATCH 05/12] Using let block for tests. --- test/test_rfft!.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_rfft!.jl b/test/test_rfft!.jl index 39a26be..60f9479 100644 --- a/test/test_rfft!.jl +++ b/test/test_rfft!.jl @@ -1,6 +1,4 @@ -a = rand(Float64,(8,4,4)) -b = PaddedRFFTArray(a) -c = copy(b) +let a = rand(Float64,(8,4,4)), b = PaddedRFFTArray(a), c = copy(b) @testset "PaddedRFFTArray creation" begin @test a == real(b) @@ -63,4 +61,5 @@ end @test c == b @test p*c ≈ rfft!(b) @test p\c ≈ irfft!(b) -end \ No newline at end of file +end +end #let block From 6fe7cd72e1574b9b8b2432d1865a152535283781 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Tue, 17 Jul 2018 16:49:14 -0400 Subject: [PATCH 06/12] Fix rand deprecation warning on tests. --- test/test_rfft!.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_rfft!.jl b/test/test_rfft!.jl index 60f9479..a4948f0 100644 --- a/test/test_rfft!.jl +++ b/test/test_rfft!.jl @@ -23,7 +23,7 @@ end @test p*c ≈ rfft!(b) @test p\c ≈ irfft!(b) - a = rand(9,4,4) + a = rand(Float64,(9,4,4)) b = PaddedRFFTArray(a) @test a == real(b) @test rfft(a) ≈ rfft!(b) @@ -36,7 +36,7 @@ end @testset "Read binary file to PaddedRFFTArray" begin for s in ((8,4,4),(9,4,4),(8,),(9,)) - aa = rand(s) + aa = rand(Float64,s) f = Base.Filesystem.tempname() write(f,aa) @test aa == real(PaddedRFFTArray(f,s)) @@ -47,7 +47,7 @@ end end @testset "brfft!" begin - a = rand(4,4) + a = rand(Float64,(4,4)) b = PaddedRFFTArray(a) rfft!(b) @test (brfft!(b) ./ 16) ≈ a From 833d973ab91257c654c05ea528c5d70998196528 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Thu, 19 Jul 2018 12:19:32 -0400 Subject: [PATCH 07/12] fix plan_brfft! flags. --- src/rfft!.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index 5261837..cd1c415 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -201,14 +201,14 @@ end # Inverse plans function plan_brfft!(X::PaddedRFFTArray{T,N}, region; - flags::Integer=PRESERVE_INPUT, + flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) - if flags&PRESERVE_INPUT != 0 + if flags&ESTIMATE != 0 + return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(X), real(X), region, flags,timelimit) + else a = similar(X) return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(a), real(a), region, flags,timelimit) - else - return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(X), real(X), region, flags,timelimit) end end From b7ec5bfd752ee9d524aeced69215c62beed09ed7 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Thu, 19 Jul 2018 12:37:09 -0400 Subject: [PATCH 08/12] Use IOBuffer instead of tem_file for tests. --- test/test_rfft!.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/test_rfft!.jl b/test/test_rfft!.jl index a4948f0..31ffa14 100644 --- a/test/test_rfft!.jl +++ b/test/test_rfft!.jl @@ -37,12 +37,13 @@ end @testset "Read binary file to PaddedRFFTArray" begin for s in ((8,4,4),(9,4,4),(8,),(9,)) aa = rand(Float64,s) - f = Base.Filesystem.tempname() + f = IOBuffer() write(f,aa) - @test aa == real(PaddedRFFTArray(f,s)) + @test aa == real(PaddedRFFTArray(seekstart(f),s)) aa = rand(Float32,s) + f = IOBuffer() write(f,aa) - @test aa == real(PaddedRFFTArray{Float32}(f,s)) + @test aa == real(PaddedRFFTArray{Float32}(seekstart(f),s)) end end From 044fc263cb778b2abe8afcb26842ba2f95e0f742 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Thu, 19 Jul 2018 14:48:59 -0400 Subject: [PATCH 09/12] Let array be overridden by plan. --- src/rfft!.jl | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index cd1c415..263d685 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -175,13 +175,7 @@ function plan_rfft!(X::PaddedRFFTArray{T,N}, region; timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) - if flags&ESTIMATE != 0 - p = rFFTWPlan{T,FORWARD,true,N}(real(X), complex_view(X), region, flags, timelimit) - else - x = similar(X) - p = rFFTWPlan{T,FORWARD,true,N}(real(x), complex_view(x), region, flags, timelimit) - end - return p + return rFFTWPlan{T,FORWARD,true,N}(real(X), complex_view(X), region, flags, timelimit) end plan_rfft!(f::PaddedRFFTArray;kws...) = plan_rfft!(f, 1:ndims(f); kws...) @@ -204,12 +198,7 @@ function plan_brfft!(X::PaddedRFFTArray{T,N}, region; flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwReal,N} (1 in region) || throw(ArgumentError("The first dimension must always be transformed")) - if flags&ESTIMATE != 0 - return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(X), real(X), region, flags,timelimit) - else - a = similar(X) - return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(a), real(a), region, flags,timelimit) - end + return rFFTWPlan{Complex{T},BACKWARD,true,N}(complex_view(X), real(X), region, flags,timelimit) end plan_brfft!(f::PaddedRFFTArray;kws...) = plan_brfft!(f,1:ndims(f);kws...) From de2657e0a1ac56bd2fa3941c2572571d2e81c902 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Thu, 19 Jul 2018 15:23:24 -0400 Subject: [PATCH 10/12] Use Colon instead of UnitRange for real view. --- src/rfft!.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index 263d685..771b8ed 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -9,27 +9,27 @@ export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irf # custom getindex and setindex! below. Hopefully, once the performance issues with ReinterpretArray # are solved we can just index the reinterpret array directly. -struct PaddedRFFTArray{T<:fftwReal,N,L} <: DenseArray{Complex{T},N} +struct PaddedRFFTArray{T<:fftwReal,N,Nm1,L} <: DenseArray{Complex{T},N} data::Array{T,N} - r::SubArray{T,N,Array{T,N},NTuple{N,UnitRange{Int}},L} # Real view skipping padding + r::SubArray{T,N,Array{T,N},Tuple{UnitRange{Int},Vararg{Base.Slice{Base.OneTo{Int}},Nm1}},L} # Real view skipping padding c::Base.ReinterpretArray{Complex{T},N,T,Array{T,N}} - function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} - rrsize = size(rr) - fsize = rrsize[1] + function PaddedRFFTArray{T,N,Nm1,L}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N,Nm1,L} + fsize = size(rr)[1] iseven(fsize) || throw( ArgumentError("First dimension of allocated array must have even number of elements")) (nx == fsize-2 || nx == fsize-1) || throw( ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) - fsize = fsize÷2 - csize = (fsize, rrsize[2:end]...) c = reinterpret(Complex{T}, rr) - rsize = (nx,rrsize[2:end]...) - r = view(rr,(1:l for l in rsize)...) - return new{T, N, N === 1 ? true : false}(rr,r,c) + r = view(rr, 1:nx, ntuple(i->Colon(),Nm1)...) + return new{T, N, Nm1, L}(rr,r,c) end # function end # struct +@generated function PaddedRFFTArray{T,N}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N} + :(PaddedRFFTArray{T,N,$(N-1),$(N === 1 ? true : false)}(rr,nx)) +end + @inline real(S::PaddedRFFTArray) = S.r @inline complex_view(S::PaddedRFFTArray) = S.c From a88455b289e8fb83acadf96ba35780c0a88c13b9 Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Fri, 20 Jul 2018 11:28:14 -0400 Subject: [PATCH 11/12] Val(Nm1) instead of Nm1 for tuple creation. --- src/rfft!.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index 771b8ed..614f1ec 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -21,7 +21,7 @@ struct PaddedRFFTArray{T<:fftwReal,N,Nm1,L} <: DenseArray{Complex{T},N} (nx == fsize-2 || nx == fsize-1) || throw( ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) c = reinterpret(Complex{T}, rr) - r = view(rr, 1:nx, ntuple(i->Colon(),Nm1)...) + r = view(rr, 1:nx, ntuple(i->Colon(),Val(Nm1))...) return new{T, N, Nm1, L}(rr,r,c) end # function end # struct From 102af34e95ad42897df67f6cb2446bcc4b4b66cd Mon Sep 17 00:00:00 2001 From: Felipe Alves Date: Fri, 20 Jul 2018 11:35:46 -0400 Subject: [PATCH 12/12] micro-optimization: OneTo for first dimension of real view. --- src/rfft!.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/rfft!.jl b/src/rfft!.jl index 614f1ec..25f3e6e 100644 --- a/src/rfft!.jl +++ b/src/rfft!.jl @@ -11,7 +11,7 @@ export PaddedRFFTArray, plan_rfft!, rfft!, plan_irfft!, plan_brfft!, brfft!, irf struct PaddedRFFTArray{T<:fftwReal,N,Nm1,L} <: DenseArray{Complex{T},N} data::Array{T,N} - r::SubArray{T,N,Array{T,N},Tuple{UnitRange{Int},Vararg{Base.Slice{Base.OneTo{Int}},Nm1}},L} # Real view skipping padding + r::SubArray{T,N,Array{T,N},Tuple{Base.OneTo{Int},Vararg{Base.Slice{Base.OneTo{Int}},Nm1}},L} # Real view skipping padding c::Base.ReinterpretArray{Complex{T},N,T,Array{T,N}} function PaddedRFFTArray{T,N,Nm1,L}(rr::Array{T,N},nx::Int) where {T<:fftwReal,N,Nm1,L} @@ -21,7 +21,7 @@ struct PaddedRFFTArray{T<:fftwReal,N,Nm1,L} <: DenseArray{Complex{T},N} (nx == fsize-2 || nx == fsize-1) || throw( ArgumentError("Number of elements on the first dimension of array must be either 1 or 2 less than the number of elements on the first dimension of the allocated array")) c = reinterpret(Complex{T}, rr) - r = view(rr, 1:nx, ntuple(i->Colon(),Val(Nm1))...) + r = view(rr, Base.OneTo(nx), ntuple(i->Colon(),Val(Nm1))...) return new{T, N, Nm1, L}(rr,r,c) end # function end # struct