11# This file was formerly a part of Julia. License is MIT: https://julialang.org/license
22
3- using Base. LinAlg: BlasReal
3+ using LinearAlgebra
4+ using LinearAlgebra: BlasReal
45import Base: show, summary, size, ndims, length, eltype,
5- * , A_mul_B!, inv, \ , A_ldiv_B!
6+ * , inv, \
67
78# DFT plan where the inputs are an array of eltype T
89abstract type Plan{T} end
@@ -33,11 +34,11 @@ realfloat(x::AbstractArray{T}) where {T<:Real} = copy1(typeof(fftfloat(zero(T)))
3334
3435# copy to a 1-based array, using circular permutation
3536function copy1 (:: Type{T} , x) where T
36- y = Array {T} (map (length, indices (x)))
37+ y = Array {T} (uninitialized, map (length, axes (x)))
3738 Base. circcopy! (y, x)
3839end
3940
40- to1 (x:: AbstractArray ) = _to1 (indices (x), x)
41+ to1 (x:: AbstractArray ) = _to1 (axes (x), x)
4142_to1 (:: Tuple{Base.OneTo,Vararg{Base.OneTo}} , x) = x
4243_to1 (:: Tuple , x) = copy1 (eltype (x), x)
4344
@@ -93,11 +94,11 @@ contains all of the information needed to compute `fft(A, dims)` quickly.
9394To apply `P` to an array `A`, use `P * A`; in general, the syntax for applying plans is much
9495like that of matrices. (A plan can only be applied to arrays of the same size as the `A`
9596for which the plan was created.) You can also apply a plan with a preallocated output array `Â`
96- by calling `A_mul_B !(Â, plan, A)`. (For `A_mul_B !`, however, the input array `A` must
97+ by calling `mul !(Â, plan, A)`. (For `mul !`, however, the input array `A` must
9798be a complex floating-point array like the output `Â`.) You can compute the inverse-transform plan by `inv(P)`
9899and apply the inverse plan with `P \\ Â` (the inverse plan is cached and reused for
99100subsequent calls to `inv` or `\\ `), and apply the inverse plan to a pre-allocated output
100- array `A` with `A_ldiv_B !(A, P, Â)`.
101+ array `A` with `ldiv !(A, P, Â)`.
101102
102103The `flags` argument is a bitwise-or of FFTW planner flags, defaulting to `FFTW.ESTIMATE`.
103104e.g. passing `FFTW.MEASURE` or `FFTW.PATIENT` will instead spend several seconds (or more)
@@ -206,12 +207,12 @@ plan_rfft(x::AbstractArray, region; kws...) = plan_rfft(realfloat(x), region; kw
206207# only require implementation to provide *(::Plan{T}, ::Array{T})
207208* (p:: Plan{T} , x:: AbstractArray ) where {T} = p * copy1 (T, x)
208209
209- # Implementations should also implement A_mul_B !(Y, plan, X) so as to support
210- # pre-allocated output arrays. We don't define * in terms of A_mul_B !
210+ # Implementations should also implement mul !(Y, plan, X) so as to support
211+ # pre-allocated output arrays. We don't define * in terms of mul !
211212# generically here, however, because of subtleties for in-place and rfft plans.
212213
213214# #############################################################################
214- # To support inv, \, and A_ldiv_B !(y, p, x), we require Plan subtypes
215+ # To support inv, \, and ldiv !(y, p, x), we require Plan subtypes
215216# to have a pinv::Plan field, which caches the inverse plan, and which
216217# should be initially undefined. They should also implement
217218# plan_inv(p) to construct the inverse of a plan p.
@@ -224,7 +225,7 @@ pinv_type(p::Plan) = eltype(_pinv_type(p))
224225inv (p:: Plan ) =
225226 isdefined (p, :pinv ) ? p. pinv:: pinv_type (p) : (p. pinv = plan_inv (p))
226227\ (p:: Plan , x:: AbstractArray ) = inv (p) * x
227- A_ldiv_B ! (y:: AbstractArray , p:: Plan , x:: AbstractArray ) = A_mul_B ! (y, inv (p), x)
228+ LinearAlgebra . ldiv ! (y:: AbstractArray , p:: Plan , x:: AbstractArray ) = LinearAlgebra . mul ! (y, inv (p), x)
228229
229230# #############################################################################
230231# implementations only need to provide the unnormalized backwards FFT,
@@ -245,7 +246,13 @@ size(p::ScaledPlan) = size(p.p)
245246show (io:: IO , p:: ScaledPlan ) = print (io, p. scale, " * " , p. p)
246247summary (p:: ScaledPlan ) = string (p. scale, " * " , summary (p. p))
247248
248- * (p:: ScaledPlan , x:: AbstractArray ) = scale! (p. p * x, p. scale)
249+ if VERSION >= v " 0.7.0-DEV.3665"
250+ * (p:: ScaledPlan , x:: AbstractArray ) = LinearAlgebra. rmul! (p. p * x, p. scale)
251+ elseif VERSION >= v " 0.7.0-DEV.3563"
252+ * (p:: ScaledPlan , x:: AbstractArray ) = LinearAlgebra. mul1! (p. p * x, p. scale)
253+ else
254+ * (p:: ScaledPlan , x:: AbstractArray ) = scale! (p. p * x, p. scale)
255+ end
249256
250257* (α:: Number , p:: Plan ) = ScaledPlan (p, α)
251258* (p:: Plan , α:: Number ) = ScaledPlan (p, α)
@@ -265,8 +272,16 @@ plan_ifft!(x::AbstractArray, region; kws...) =
265272
266273plan_inv (p:: ScaledPlan ) = ScaledPlan (plan_inv (p. p), inv (p. scale))
267274
268- A_mul_B! (y:: AbstractArray , p:: ScaledPlan , x:: AbstractArray ) =
269- scale! (p. scale, A_mul_B! (y, p. p, x))
275+ if VERSION >= v " 0.7.0-DEV.3665"
276+ LinearAlgebra. mul! (y:: AbstractArray , p:: ScaledPlan , x:: AbstractArray ) =
277+ LinearAlgebra. lmul! (p. scale, LinearAlgebra. mul! (y, p. p, x))
278+ elseif VERSION >= v " 0.7.0-DEV.3563"
279+ LinearAlgebra. mul! (y:: AbstractArray , p:: ScaledPlan , x:: AbstractArray ) =
280+ LinearAlgebra. mul2! (p. scale, LinearAlgebra. mul! (y, p. p, x))
281+ else
282+ LinearAlgebra. mul! (y:: AbstractArray , p:: ScaledPlan , x:: AbstractArray ) =
283+ scale! (p. scale, LinearAlgebra. mul! (y, p. p, x))
284+ end
270285
271286# #############################################################################
272287# Real-input DFTs are annoying because the output has a different size
0 commit comments