From 3f4a89bb91cf28f9cc8d3f1af82d7463e17ff64b Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Sun, 2 May 2021 13:36:29 +0800 Subject: [PATCH 01/11] Fix for revertation to "fftw" With out `force = true`, the revertation from "mkl" to "fftw" won't work. ```julia ERROR: ArgumentError: Cannot set preference 'provider' to 'fftw' for FFTW in C:\Users\XXX\.julia\environments\v1.6\LocalPreferences.toml: preference already set to 'mkl'! ``` --- src/providers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/providers.jl b/src/providers.jl index 7f8c4c8..01933bf 100644 --- a/src/providers.jl +++ b/src/providers.jl @@ -31,7 +31,7 @@ function set_provider!(provider; export_prefs::Bool = false) if provider !== nothing && provider !== missing && provider ∉ ("fftw", "mkl") throw(ArgumentError("Invalid provider value '$(provider)'")) end - set_preferences!(@__MODULE__, "provider" => provider; export_prefs) + set_preferences!(@__MODULE__, "provider" => provider; export_prefs, force = true) if provider != fftw_provider # Re-fetch to get default values in the event that `nothing` or `missing` was passed in. provider = get_provider() From 85bef3afe3fb7e774c052f8dc54636139362fb57 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Sun, 2 May 2021 14:57:21 +0800 Subject: [PATCH 02/11] make `UNALIGNED` on for MKL --- src/fft.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/fft.jl b/src/fft.jl index e59342f..dcc2940 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -76,6 +76,11 @@ const ESTIMATE = UInt32(1 << 6) const WISDOM_ONLY = UInt32(1 << 21) const NO_SIMD = UInt32(1 << 17) # disable SIMD, useful for benchmarking +@static if fftw_provider == "mkl" + addflags(flags) = UNALIGNED | flags +else + addflags(flags) = flags +end ## R2R transform kinds const R2HC = 0 @@ -682,13 +687,13 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} cFFTWPlan{T,$direction,false,N}(X, fakesimilar(flags, X, T), - region, flags, timelimit) + region, addflags(flags), timelimit) end function $plan_f!(X::StridedArray{T,N}, region; flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} - cFFTWPlan{T,$direction,true,N}(X, X, region, flags, timelimit) + cFFTWPlan{T,$direction,true,N}(X, X, region, addflags(flags), timelimit) end $plan_f(X::StridedArray{<:fftwComplex}; kws...) = $plan_f(X, 1:ndims(X); kws...) @@ -699,7 +704,7 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) X = Array{T}(undef, p.sz) Y = inplace ? X : fakesimilar(p.flags, X, T) ScaledPlan(cFFTWPlan{T,$idirection,inplace,N}(X, Y, p.region, - p.flags, NO_TIMELIMIT), + addflags(p.flags), NO_TIMELIMIT), normalization(X, p.region)) end end From f7f2e6f0b8e8db9e54406b4b5ed16b5dd44a36b1 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 14:13:59 +0800 Subject: [PATCH 03/11] Update fft.jl omit alignment check for mkl --- src/fft.jl | 80 ++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 50 insertions(+), 30 deletions(-) diff --git a/src/fft.jl b/src/fft.jl index dcc2940..cea6c31 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -76,11 +76,6 @@ const ESTIMATE = UInt32(1 << 6) const WISDOM_ONLY = UInt32(1 << 21) const NO_SIMD = UInt32(1 << 17) # disable SIMD, useful for benchmarking -@static if fftw_provider == "mkl" - addflags(flags) = UNALIGNED | flags -else - addflags(flags) = flags -end ## R2R transform kinds const R2HC = 0 @@ -422,30 +417,55 @@ end # Check whether a FFTWPlan is applicable to a given input array, and # throw an informative error if not: -function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T - if size(X) != p.sz - throw(ArgumentError("FFTW plan applied to wrong-size array")) - elseif strides(X) != p.istride - throw(ArgumentError("FFTW plan applied to wrong-strides array")) - elseif alignment_of(X) != p.ialign && p.flags & UNALIGNED == 0 - throw(ArgumentError("FFTW plan applied to array with wrong memory alignment")) +@static if fftw_provider == "fftw" + function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T + if size(X) != p.sz + throw(ArgumentError("FFTW plan applied to wrong-size array")) + elseif strides(X) != p.istride + throw(ArgumentError("FFTW plan applied to wrong-strides array")) + elseif alignment_of(X) != p.ialign && p.flags & UNALIGNED == 0 + throw(ArgumentError("FFTW plan applied to array with wrong memory alignment")) + end + end + + function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} + assert_applicable(p, X) + if size(Y) != p.osz + throw(ArgumentError("FFTW plan applied to wrong-size output")) + elseif strides(Y) != p.ostride + throw(ArgumentError("FFTW plan applied to wrong-strides output")) + elseif alignment_of(Y) != p.oalign && p.flags & UNALIGNED == 0 + throw(ArgumentError("FFTW plan applied to output with wrong memory alignment")) + elseif inplace != (pointer(X) == pointer(Y)) + throw(ArgumentError(string("FFTW ", + inplace ? "in-place" : "out-of-place", + " plan applied to ", + inplace ? "out-of-place" : "in-place", + " data"))) + end + end +else + function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T + if size(X) != p.sz + throw(ArgumentError("FFTW plan applied to wrong-size array")) + elseif strides(X) != p.istride + throw(ArgumentError("FFTW plan applied to wrong-strides array")) + end end -end -function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} - assert_applicable(p, X) - if size(Y) != p.osz - throw(ArgumentError("FFTW plan applied to wrong-size output")) - elseif strides(Y) != p.ostride - throw(ArgumentError("FFTW plan applied to wrong-strides output")) - elseif alignment_of(Y) != p.oalign && p.flags & UNALIGNED == 0 - throw(ArgumentError("FFTW plan applied to output with wrong memory alignment")) - elseif inplace != (pointer(X) == pointer(Y)) - throw(ArgumentError(string("FFTW ", - inplace ? "in-place" : "out-of-place", - " plan applied to ", - inplace ? "out-of-place" : "in-place", - " data"))) + function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} + assert_applicable(p, X) + if size(Y) != p.osz + throw(ArgumentError("FFTW plan applied to wrong-size output")) + elseif strides(Y) != p.ostride + throw(ArgumentError("FFTW plan applied to wrong-strides output")) + elseif inplace != (pointer(X) == pointer(Y)) + throw(ArgumentError(string("FFTW ", + inplace ? "in-place" : "out-of-place", + " plan applied to ", + inplace ? "out-of-place" : "in-place", + " data"))) + end end end @@ -687,13 +707,13 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} cFFTWPlan{T,$direction,false,N}(X, fakesimilar(flags, X, T), - region, addflags(flags), timelimit) + region, flags, timelimit) end function $plan_f!(X::StridedArray{T,N}, region; flags::Integer=ESTIMATE, timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N} - cFFTWPlan{T,$direction,true,N}(X, X, region, addflags(flags), timelimit) + cFFTWPlan{T,$direction,true,N}(X, X, region, flags, timelimit) end $plan_f(X::StridedArray{<:fftwComplex}; kws...) = $plan_f(X, 1:ndims(X); kws...) @@ -704,7 +724,7 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) X = Array{T}(undef, p.sz) Y = inplace ? X : fakesimilar(p.flags, X, T) ScaledPlan(cFFTWPlan{T,$idirection,inplace,N}(X, Y, p.region, - addflags(p.flags), NO_TIMELIMIT), + p.flags, NO_TIMELIMIT), normalization(X, p.region)) end end From 18162c1adcb66d89227738fa1306b0887cbe139d Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 19:41:10 +0800 Subject: [PATCH 04/11] Create mkl_walkplan.jl MKL's FFTW3 wrapper don't support multidimension vector fft (ndims(x) must <= length(region) + 1). This file implent a walk plan to support more multidimension fft. --- src/mkl_walkplan.jl | 116 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 src/mkl_walkplan.jl diff --git a/src/mkl_walkplan.jl b/src/mkl_walkplan.jl new file mode 100644 index 0000000..77d7e7f --- /dev/null +++ b/src/mkl_walkplan.jl @@ -0,0 +1,116 @@ +struct PickDim{T,N,P<:StridedArray{T}} <: DenseArray{T,N} + data::P + dims::NTuple{N,Int} + PickDim(data::StridedArray{T}, dims::NTuple{N,Int}) where{T,N} = begin + new{T,N,typeof(data)}(data, dims) + end +end +Base.size(A::PickDim) = map(i -> size(A.data)[i], A.dims) +Base.strides(A::PickDim) = map(i -> strides(A.data)[i], A.dims) + +mutable struct cWalkPlan{T,K,inplace,N,NW} <: FFTWPlan{T,K,inplace} + plan::PlanPtr + sz::NTuple{N,Int} + osz::NTuple{N,Int} + istride::NTuple{N,Int} + ostride::NTuple{N,Int} + walkdims::NTuple{NW,Int} + function cWalkPlan{T,K,inplace}(plan::PlanPtr, X, Y, walkdims) where {T,K,inplace} + p = new{T,K,inplace,ndims(X),length(walkdims)}(plan, + size(X), size(Y), strides(X), strides(Y), walkdims) + finalizer(maybe_destroy_plan, p) + p + end +end +for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3), + (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f)) + @eval @exclusive function cWalkPlan{$Tc,direction,inplace}(X::StridedArray{$Tc,N}, + Y::StridedArray{$Tc,N}, + region, timelimit::Real) where {direction,inplace,N} + length(region) == length(unique(region)) || + throw(ArgumentError("each dimension can be transformed at most once")) + plandims, walkdims = split_dim(X, region) + X′, Y′ = PickDim(X, plandims), PickDim(Y, plandims) + dims, howmany = dims_howmany(X′, Y′, [size(X′)...], 1:length(region)) + plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib), + PlanPtr, + (Int32, Ptr{Int}, Int32, Ptr{Int}, + Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32), + size(dims,2), dims, size(howmany,2), howmany, + X, Y, direction, 0) + if plan == C_NULL + error("FFTW could not create plan") # shouldn't normally happen + end + return cWalkPlan{$Tc,K,inplace}(plan, X, Y, walkdims) + end +end + +function split_dim(X, region) + dims = filter(i -> !in(i, region), 1:ndims(X)) + # Determine which dimension to be planned by MKL's Dfti + # I believe the first dimension shoule be selected anyhow for columnmajor layout. + # Other dimensions might be selected based on size for better thread performance. + id = !in(1,region) || size(X,dims[1]) > 10Threads.nthreads() ? + 1 : argmax(map(i -> size(X,i), dims)) + walkdims = dims[[1:id-1;id+1:end]] |> NTuple{ndims(X) - length(region) - 1, Int} + plandims = (region...,dims[id]) + plandims, walkdims +end + +show(io::IO, p::cWalkPlan{T,K,inplace}) where {T,K,inplace} = begin + print(io, inplace ? "FFTW in-place " : "FFTW ", + K < 0 ? "forward" : "backward", " plan for ") + showfftdims(io, p.sz, p.istride, T) + has_sprint_plan && print(io, "\n", sprint_plan(p)) +end + +unsafe_walk_execute!(plan::cWalkPlan{T}, X::Ptr{T}, Y::Ptr{T}) where {T<:fftwComplex} = begin + if T <:fftwSingle + @ccall libfftw3f.fftwf_execute_dft(plan::PlanPtr, X::Ptr{T}, Y::Ptr{T})::Cvoid + else + @ccall libfftw3.fftw_execute_dft(plan::PlanPtr, X::Ptr{T}, Y::Ptr{T})::Cvoid + end +end + +function unsafe_execute!(p::cWalkPlan{T}, x::StridedArray{T}, + y::StridedArray{T}) where {T <: fftwComplex} + walkistride = map(i -> p.istride[i], p.walkdims) + walkostride = map(i -> p.ostride[i], p.walkdims) + ax = map(i -> axes(x,i) .- first(axes(x,i)), p.walkdims) + Elsz, pointerˣ, pointerʸ = sizeof(T), pointer(x), pointer(y) + for ind in Iterators.product(ax...) + pointerˣ′ = pointerˣ + Elsz * sum(walkistride .* ind) + pointerʸ′ = pointerʸ + Elsz * sum(walkostride .* ind) + unsafe_walk_execute!(p, pointerˣ′, pointerʸ′) + end +end + +function mul!(y::StridedArray{T,N}, p::cWalkPlan{T}, x::StridedArray{T,N}) where {T,N} + assert_applicable(p, x, y) + unsafe_execute!(p, x, y) + y +end + +function *(p::cWalkPlan{T,K,inplace}, x::StridedArray{T,N}) where {T,K,N,inplace} + assert_applicable(p, x) + y = inplace ? x : Array{T,N}(undef, p.osz) + unsafe_execute!(p, x, y) + y +end + +const MaybeHighDimension{T} = Union{StridedArray{T,3}, StridedArray{T,4}, StridedArray{T,5}, StridedArray{T,6}, StridedArray{T,7}} +for inplace in (false,true) + for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) + plan_f = inplace ? Symbol("plan_",f,"!") : Symbol("plan_",f) + @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}, region; kws...) = $plan_f(X, Tuple(region); kws...) + @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}; kws...) = $plan_f(X, ntuple(identity, ndims(X)); kws...) + @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}, region::Union{Integer,Tuple}; + flags::Integer=ESTIMATE, + timelimit::Real=NO_TIMELIMIT) = begin + T, N = eltype(X), ndims(X) + X′ = $inplace ? X : FakeArray{T}(size(X)) + N <= length(region) + 1 && return cFFTWPlan{T,$direction,$inplace,N}(X, X′, region, flags, timelimit) + cWalkPlan{T,$direction,$inplace}(X, X′, region) + end + end +end From d7534c3043404d07c1af15c253d3cb71a4f2af63 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 19:42:24 +0800 Subject: [PATCH 05/11] Update FFTW.jl --- src/FFTW.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/FFTW.jl b/src/FFTW.jl index ec92095..e6466c5 100644 --- a/src/FFTW.jl +++ b/src/FFTW.jl @@ -60,6 +60,9 @@ end include("fft.jl") include("dct.jl") +@static if fftw_provider == "mkl" + include("mkl_walkplan.jl") +end include("precompile.jl") _precompile_() From 82769d663cf77de6d633767b2556f93764facc0a Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 19:45:28 +0800 Subject: [PATCH 06/11] Update mkl_walkplan.jl --- src/mkl_walkplan.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mkl_walkplan.jl b/src/mkl_walkplan.jl index 77d7e7f..e70ed0f 100644 --- a/src/mkl_walkplan.jl +++ b/src/mkl_walkplan.jl @@ -25,8 +25,7 @@ end for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3), (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f)) @eval @exclusive function cWalkPlan{$Tc,direction,inplace}(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tc,N}, - region, timelimit::Real) where {direction,inplace,N} + Y::StridedArray{$Tc,N},region) where {direction,inplace,N} length(region) == length(unique(region)) || throw(ArgumentError("each dimension can be transformed at most once")) plandims, walkdims = split_dim(X, region) From 23d05708d0d8e8ee1f9e5202cc8434d0ed86ddbc Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 19:50:25 +0800 Subject: [PATCH 07/11] Update mkl_walkplan.jl --- src/mkl_walkplan.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mkl_walkplan.jl b/src/mkl_walkplan.jl index e70ed0f..fdb427f 100644 --- a/src/mkl_walkplan.jl +++ b/src/mkl_walkplan.jl @@ -24,8 +24,8 @@ mutable struct cWalkPlan{T,K,inplace,N,NW} <: FFTWPlan{T,K,inplace} end for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3), (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f)) - @eval @exclusive function cWalkPlan{$Tc,direction,inplace}(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tc,N},region) where {direction,inplace,N} + @eval @exclusive function cWalkPlan{$Tc,K,inplace}(X::StridedArray{$Tc,N}, + Y::StridedArray{$Tc,N},region) where {K,inplace,N} length(region) == length(unique(region)) || throw(ArgumentError("each dimension can be transformed at most once")) plandims, walkdims = split_dim(X, region) @@ -36,7 +36,7 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3), (Int32, Ptr{Int}, Int32, Ptr{Int}, Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32), size(dims,2), dims, size(howmany,2), howmany, - X, Y, direction, 0) + X, Y, K, ESTIMATE) ## flags is useless if plan == C_NULL error("FFTW could not create plan") # shouldn't normally happen end From e508b45a6759fc0a061e7316133c696e2812328d Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 20:14:49 +0800 Subject: [PATCH 08/11] Update FFTW.jl --- src/FFTW.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/FFTW.jl b/src/FFTW.jl index e6466c5..ec92095 100644 --- a/src/FFTW.jl +++ b/src/FFTW.jl @@ -60,9 +60,6 @@ end include("fft.jl") include("dct.jl") -@static if fftw_provider == "mkl" - include("mkl_walkplan.jl") -end include("precompile.jl") _precompile_() From bb86fe882943834d1ab233d76d0ea1e1bddbc852 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 20:16:23 +0800 Subject: [PATCH 09/11] Update fft.jl --- src/fft.jl | 68 ++++++++++++++++++------------------------------------ 1 file changed, 22 insertions(+), 46 deletions(-) diff --git a/src/fft.jl b/src/fft.jl index cea6c31..cdcf61e 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -417,55 +417,31 @@ end # Check whether a FFTWPlan is applicable to a given input array, and # throw an informative error if not: -@static if fftw_provider == "fftw" - function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T - if size(X) != p.sz - throw(ArgumentError("FFTW plan applied to wrong-size array")) - elseif strides(X) != p.istride - throw(ArgumentError("FFTW plan applied to wrong-strides array")) - elseif alignment_of(X) != p.ialign && p.flags & UNALIGNED == 0 - throw(ArgumentError("FFTW plan applied to array with wrong memory alignment")) - end - end - function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} - assert_applicable(p, X) - if size(Y) != p.osz - throw(ArgumentError("FFTW plan applied to wrong-size output")) - elseif strides(Y) != p.ostride - throw(ArgumentError("FFTW plan applied to wrong-strides output")) - elseif alignment_of(Y) != p.oalign && p.flags & UNALIGNED == 0 - throw(ArgumentError("FFTW plan applied to output with wrong memory alignment")) - elseif inplace != (pointer(X) == pointer(Y)) - throw(ArgumentError(string("FFTW ", - inplace ? "in-place" : "out-of-place", - " plan applied to ", - inplace ? "out-of-place" : "in-place", - " data"))) - end - end -else - function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T - if size(X) != p.sz - throw(ArgumentError("FFTW plan applied to wrong-size array")) - elseif strides(X) != p.istride - throw(ArgumentError("FFTW plan applied to wrong-strides array")) - end +function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T + if size(X) != p.sz + throw(ArgumentError("FFTW plan applied to wrong-size array")) + elseif strides(X) != p.istride + throw(ArgumentError("FFTW plan applied to wrong-strides array")) + elseif alignment_of(X) != p.ialign && p.flags & UNALIGNED == 0 + throw(ArgumentError("FFTW plan applied to array with wrong memory alignment")) end +end - function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} - assert_applicable(p, X) - if size(Y) != p.osz - throw(ArgumentError("FFTW plan applied to wrong-size output")) - elseif strides(Y) != p.ostride - throw(ArgumentError("FFTW plan applied to wrong-strides output")) - elseif inplace != (pointer(X) == pointer(Y)) - throw(ArgumentError(string("FFTW ", - inplace ? "in-place" : "out-of-place", - " plan applied to ", - inplace ? "out-of-place" : "in-place", - " data"))) - end +function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace} + assert_applicable(p, X) + if size(Y) != p.osz + throw(ArgumentError("FFTW plan applied to wrong-size output")) + elseif strides(Y) != p.ostride + throw(ArgumentError("FFTW plan applied to wrong-strides output")) + elseif alignment_of(Y) != p.oalign && p.flags & UNALIGNED == 0 + throw(ArgumentError("FFTW plan applied to output with wrong memory alignment")) + elseif inplace != (pointer(X) == pointer(Y)) + throw(ArgumentError(string("FFTW ", + inplace ? "in-place" : "out-of-place", + " plan applied to ", + inplace ? "out-of-place" : "in-place", + " data"))) end end From 977cfa5257deac3e7832fcf980326bfd7d130bbd Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 20:18:05 +0800 Subject: [PATCH 10/11] Delete mkl_walkplan.jl --- src/mkl_walkplan.jl | 115 -------------------------------------------- 1 file changed, 115 deletions(-) delete mode 100644 src/mkl_walkplan.jl diff --git a/src/mkl_walkplan.jl b/src/mkl_walkplan.jl deleted file mode 100644 index fdb427f..0000000 --- a/src/mkl_walkplan.jl +++ /dev/null @@ -1,115 +0,0 @@ -struct PickDim{T,N,P<:StridedArray{T}} <: DenseArray{T,N} - data::P - dims::NTuple{N,Int} - PickDim(data::StridedArray{T}, dims::NTuple{N,Int}) where{T,N} = begin - new{T,N,typeof(data)}(data, dims) - end -end -Base.size(A::PickDim) = map(i -> size(A.data)[i], A.dims) -Base.strides(A::PickDim) = map(i -> strides(A.data)[i], A.dims) - -mutable struct cWalkPlan{T,K,inplace,N,NW} <: FFTWPlan{T,K,inplace} - plan::PlanPtr - sz::NTuple{N,Int} - osz::NTuple{N,Int} - istride::NTuple{N,Int} - ostride::NTuple{N,Int} - walkdims::NTuple{NW,Int} - function cWalkPlan{T,K,inplace}(plan::PlanPtr, X, Y, walkdims) where {T,K,inplace} - p = new{T,K,inplace,ndims(X),length(walkdims)}(plan, - size(X), size(Y), strides(X), strides(Y), walkdims) - finalizer(maybe_destroy_plan, p) - p - end -end -for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3), - (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f)) - @eval @exclusive function cWalkPlan{$Tc,K,inplace}(X::StridedArray{$Tc,N}, - Y::StridedArray{$Tc,N},region) where {K,inplace,N} - length(region) == length(unique(region)) || - throw(ArgumentError("each dimension can be transformed at most once")) - plandims, walkdims = split_dim(X, region) - X′, Y′ = PickDim(X, plandims), PickDim(Y, plandims) - dims, howmany = dims_howmany(X′, Y′, [size(X′)...], 1:length(region)) - plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib), - PlanPtr, - (Int32, Ptr{Int}, Int32, Ptr{Int}, - Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32), - size(dims,2), dims, size(howmany,2), howmany, - X, Y, K, ESTIMATE) ## flags is useless - if plan == C_NULL - error("FFTW could not create plan") # shouldn't normally happen - end - return cWalkPlan{$Tc,K,inplace}(plan, X, Y, walkdims) - end -end - -function split_dim(X, region) - dims = filter(i -> !in(i, region), 1:ndims(X)) - # Determine which dimension to be planned by MKL's Dfti - # I believe the first dimension shoule be selected anyhow for columnmajor layout. - # Other dimensions might be selected based on size for better thread performance. - id = !in(1,region) || size(X,dims[1]) > 10Threads.nthreads() ? - 1 : argmax(map(i -> size(X,i), dims)) - walkdims = dims[[1:id-1;id+1:end]] |> NTuple{ndims(X) - length(region) - 1, Int} - plandims = (region...,dims[id]) - plandims, walkdims -end - -show(io::IO, p::cWalkPlan{T,K,inplace}) where {T,K,inplace} = begin - print(io, inplace ? "FFTW in-place " : "FFTW ", - K < 0 ? "forward" : "backward", " plan for ") - showfftdims(io, p.sz, p.istride, T) - has_sprint_plan && print(io, "\n", sprint_plan(p)) -end - -unsafe_walk_execute!(plan::cWalkPlan{T}, X::Ptr{T}, Y::Ptr{T}) where {T<:fftwComplex} = begin - if T <:fftwSingle - @ccall libfftw3f.fftwf_execute_dft(plan::PlanPtr, X::Ptr{T}, Y::Ptr{T})::Cvoid - else - @ccall libfftw3.fftw_execute_dft(plan::PlanPtr, X::Ptr{T}, Y::Ptr{T})::Cvoid - end -end - -function unsafe_execute!(p::cWalkPlan{T}, x::StridedArray{T}, - y::StridedArray{T}) where {T <: fftwComplex} - walkistride = map(i -> p.istride[i], p.walkdims) - walkostride = map(i -> p.ostride[i], p.walkdims) - ax = map(i -> axes(x,i) .- first(axes(x,i)), p.walkdims) - Elsz, pointerˣ, pointerʸ = sizeof(T), pointer(x), pointer(y) - for ind in Iterators.product(ax...) - pointerˣ′ = pointerˣ + Elsz * sum(walkistride .* ind) - pointerʸ′ = pointerʸ + Elsz * sum(walkostride .* ind) - unsafe_walk_execute!(p, pointerˣ′, pointerʸ′) - end -end - -function mul!(y::StridedArray{T,N}, p::cWalkPlan{T}, x::StridedArray{T,N}) where {T,N} - assert_applicable(p, x, y) - unsafe_execute!(p, x, y) - y -end - -function *(p::cWalkPlan{T,K,inplace}, x::StridedArray{T,N}) where {T,K,N,inplace} - assert_applicable(p, x) - y = inplace ? x : Array{T,N}(undef, p.osz) - unsafe_execute!(p, x, y) - y -end - -const MaybeHighDimension{T} = Union{StridedArray{T,3}, StridedArray{T,4}, StridedArray{T,5}, StridedArray{T,6}, StridedArray{T,7}} -for inplace in (false,true) - for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD)) - plan_f = inplace ? Symbol("plan_",f,"!") : Symbol("plan_",f) - @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}, region; kws...) = $plan_f(X, Tuple(region); kws...) - @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}; kws...) = $plan_f(X, ntuple(identity, ndims(X)); kws...) - @eval $plan_f(X::MaybeHighDimension{<:fftwComplex}, region::Union{Integer,Tuple}; - flags::Integer=ESTIMATE, - timelimit::Real=NO_TIMELIMIT) = begin - T, N = eltype(X), ndims(X) - X′ = $inplace ? X : FakeArray{T}(size(X)) - N <= length(region) + 1 && return cFFTWPlan{T,$direction,$inplace,N}(X, X′, region, flags, timelimit) - cWalkPlan{T,$direction,$inplace}(X, X′, region) - end - end -end From b7f239928142737cf1d1cdc7dc20f350f373e652 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Mon, 3 May 2021 20:18:40 +0800 Subject: [PATCH 11/11] Update fft.jl --- src/fft.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/fft.jl b/src/fft.jl index cdcf61e..e59342f 100644 --- a/src/fft.jl +++ b/src/fft.jl @@ -417,7 +417,6 @@ end # Check whether a FFTWPlan is applicable to a given input array, and # throw an informative error if not: - function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T if size(X) != p.sz throw(ArgumentError("FFTW plan applied to wrong-size array"))