From aba914dad20682b0330d91cd9eef3a1935ada3f0 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Tue, 28 Sep 2021 17:13:38 -0400
Subject: [PATCH 1/4] Partly works? needs optimizing

---
 src/FFTW.jl |   2 +-
 src/fft.jl  | 280 ++++++++++++++++++++++++++++++++++++++++++----------
 2 files changed, 227 insertions(+), 55 deletions(-)

diff --git a/src/FFTW.jl b/src/FFTW.jl
index 4366ee7..0493cbc 100644
--- a/src/FFTW.jl
+++ b/src/FFTW.jl
@@ -12,7 +12,7 @@ import AbstractFFTs: Plan, ScaledPlan,
                      rfft_output_size, brfft_output_size,
                      plan_inv, normalization
 
-export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!
+export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!, brfft!, irfft!, plan_brfft!, plan_irfft!
 
 include("providers.jl")
 
diff --git a/src/fft.jl b/src/fft.jl
index 4063ea7..94929b6 100644
--- a/src/fft.jl
+++ b/src/fft.jl
@@ -1,8 +1,10 @@
 # This file was formerly a part of Julia. License is MIT: https://julialang.org/license
 
-import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer
+import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, copy, similar, eltype, \
 import LinearAlgebra: mul!
 
+
+
 """
     r2r(A, kind [, dims])
 
@@ -108,6 +110,74 @@ const fftwSingle = Union{Float32,Complex{Float32}}
 const fftwTypeDouble = Union{Type{Float64},Type{Complex{Float64}}}
 const fftwTypeSingle = Union{Type{Float32},Type{Complex{Float32}}}
 
+struct PaddedRFFTArray{T, N, A, Ar, Ac} <: AbstractArray{T, N} where {A <: Union{AbstractArray{T, N}, AbstractArray{Complex{T}, N}}, Ar <: AbstractArray{T, N}, Ac <: AbstractArray{Complex{T}, N}}
+    parent::A
+    r::Ar # real view
+    c::Ac # complex view
+
+    function PaddedRFFTArray(data::StridedArray{Tc, N}) where {Tc <: fftwComplex, N}
+        fsize = 2 * (size(data, 1) - 1)
+        r = view(reinterpret(real(Tc), data), Base.OneTo(fsize), ntuple(i -> Colon(), Val(ndims(data) - 1))...)
+        new{real(Tc), N, typeof(data), typeof(r), typeof(data)}(data, r, data)
+    end
+
+    function PaddedRFFTArray(data::StridedArray{T, N}, nx) where {T <: fftwReal, N}
+        fsize = size(data, 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"))
+        c = reinterpret(Complex{T}, data)
+        r = view(data, Base.OneTo(nx), ntuple(i->Colon(), Val(ndims(data) - 1))...)
+        new{T, N, typeof(data), typeof(r), typeof(c)}(data, r, c)
+    end
+end
+
+FFTStridedArray{T, N} = Union{StridedArray{T, N}, PaddedRFFTArray{T, N}}
+
+@inline real_view(S::PaddedRFFTArray) = S.r
+
+@inline complex_view(S::PaddedRFFTArray) = S.c
+
+copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(parent(S)), size(real_view(S), 1))
+
+# similar(f::PaddedRFFTArray, ::Type{T}, dims::Tuple{Vararg{Int, N}}) where {T <: fftwReal, N} =
+#     PaddedRFFTArray{T, N}(dims)
+# similar(f::PaddedRFFTArray{T}, dims::Vararg{Tuple{Int, N}}) where {T, N} =
+#     PaddedRFFTArray{T, N}(dims)
+# similar(f::PaddedRFFTArray, ::Type{T}) where {T <: fftwReal} =
+#     PaddedRFFTArray{T}(size(real_view(f)))
+# similar(f::PaddedRFFTArray{T, N}) where {T,N} = 
+#     PaddedRFFTArray{T, N}(similar(parent(f)), size(real_view(f), 1)) 
+
+# size(S::PaddedRFFTArray) =
+#     size(complex_view(S))
+
+# IndexStyle(::Type{T}) where {T<:PaddedRFFTArray} = 
+#     IndexLinear()
+
+function PaddedRFFTArray{T}(ndims::Vararg{Integer, N}) where {T <: fftwReal, N}
+    fsize = (ndims[1] ÷ 2 + 1) * 2
+    a = zeros(T, (fsize, ndims[2:end]...))
+    PaddedRFFTArray{T, N}(a, ndims[1])
+end
+
+PaddedRFFTArray{T}(ndims::NTuple{N, Integer}) where {T <: fftwReal, 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 copyto!(t.r, a)
+    return t
+end
+
+
 # For ESTIMATE plans, FFTW allows one to pass NULL for the array pointer,
 # since it is not written to.  Hence, it is convenient to create an
 # array-like type that carries a size and a stride like a "real" array
@@ -138,6 +208,67 @@ alignment_of(A::FakeArray) = Int32(0)
 # around FFTW's internal file i/o buffering [see the BUFSZ constant in
 # FFTW's api/import-wisdom-from-file.c file].
 
+#==
+a = rand(10002, 80000)
+const cccc = rfft(a[1:10000, :])
+function e()
+       d = copy(cccc)
+       irfft!(d, 10000)
+       end
+
+Overhead ╎ [+additional indent] Count File:Line; Function
+=========================================================
+     ╎40027 @Base\task.jl:411; (::VSCodeServer.var"#53#54")()
+     ╎ 40027 @VSCodeServer\src\eval.jl:34; macro expansion
+     ╎  40027 @Base\essentials.jl:706; invokelatest(::Any)
+     ╎   40027 @Base\essentials.jl:708; #invokelatest#2
+     ╎    40027 @VSCodeServer\src\repl.jl:124; (::VSCodeServer.var"#68#70"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt})()
+     ╎     40027 @Base\logging.jl:603; with_logger
+     ╎    ╎ 40027 @Base\logging.jl:491; with_logstate(f::Function, logstate::Any)
+     ╎    ╎  40027 @VSCodeServer\src\repl.jl:123; (::VSCodeServer.var"#69#71"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt})()
+     ╎    ╎   40027 @VSCodeServer\src\repl.jl:157; repleval(m::Module, code::Expr, #unused#::String)
+     ╎    ╎    40027 @Base\Base.jl:39; eval
+     ╎    ╎     40027 @Base\boot.jl:360; eval
+     ╎    ╎    ╎ 931   REPL[40]:2; e()
+  931╎    ╎    ╎  931   @Base\array.jl:349; copy
+     ╎    ╎    ╎ 39087 REPL[40]:3; e()
+     ╎    ╎    ╎  39087 @FFTW\src\fft.jl:812; irfft!(x::Matrix{ComplexF64}, d::Int64)
+     ╎    ╎    ╎   35687 @FFTW\src\fft.jl:923; *
+     ╎    ╎    ╎    35687 @FFTW\src\fft.jl:917; *
+     ╎    ╎    ╎     35687 @FFTW\src\fft.jl:892; mul!
+35687╎    ╎    ╎    ╎ 35687 @FFTW\src\fft.jl:556; unsafe_execute!(plan::FFTW.rFFTWPlan{ComplexF64, 1, true, 2, UnitRange{Int64}}, X::Matrix{ComplexF64}, Y::SubArray{Float64, 2...
+     ╎    ╎    ╎   3400  @FFTW\src\fft.jl:924; *
+     ╎    ╎    ╎    3400  ...worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\generic.jl:182; rmul!(X::SubArray{Float64, 2, Base.ReinterpretArray{Float64, 2, ComplexF64, Matrix{ComplexF64}, false}, Tuple{Base.OneTo{Int64...
+     ╎    ╎    ╎     257   @Base\simdloop.jl:75; macro expansion
+  257╎    ╎    ╎    ╎ 257   @Base\int.jl:83; <
+     ╎    ╎    ╎     3143  @Base\simdloop.jl:77; macro expansion
+     ╎    ╎    ╎    ╎ 3143  ...orker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\generic.jl:183; macro expansion
+     ╎    ╎    ╎    ╎  731   @Base\abstractarray.jl:1170; getindex
+     ╎    ╎    ╎    ╎   731   @Base\abstractarray.jl:1214; _getindex
+     ╎    ╎    ╎    ╎    731   @Base\subarray.jl:276; getindex
+     ╎    ╎    ╎    ╎     731   @Base\reinterpretarray.jl:320; getindex
+     ╎    ╎    ╎    ╎    ╎ 368   @Base\reinterpretarray.jl:355; _getindex_ra
+     ╎    ╎    ╎    ╎    ╎  368   @Base\div.jl:120; divrem
+     ╎    ╎    ╎    ╎    ╎   368   @Base\div.jl:124; divrem
+  368╎    ╎    ╎    ╎    ╎    368   @Base\int.jl:262; rem
+     ╎    ╎    ╎    ╎    ╎ 363   @Base\reinterpretarray.jl:371; _getindex_ra
+  259╎    ╎    ╎    ╎    ╎  259   @Base\array.jl:806; getindex
+     ╎    ╎    ╎    ╎    ╎  104   @Base\refvalue.jl:57; setindex!
+  104╎    ╎    ╎    ╎    ╎   104   @Base\Base.jl:34; setproperty!
+     ╎    ╎    ╎    ╎  2404  @Base\abstractarray.jl:1267; setindex!
+     ╎    ╎    ╎    ╎   2404  @Base\abstractarray.jl:1297; _setindex!
+     ╎    ╎    ╎    ╎    2404  @Base\subarray.jl:317; setindex!
+     ╎    ╎    ╎    ╎     2404  @Base\reinterpretarray.jl:446; setindex!
+     ╎    ╎    ╎    ╎    ╎ 548   @Base\reinterpretarray.jl:495; _setindex_ra!
+     ╎    ╎    ╎    ╎    ╎  548   @Base\refvalue.jl:57; setindex!
+  548╎    ╎    ╎    ╎    ╎   548   @Base\Base.jl:34; setproperty!
+     ╎    ╎    ╎    ╎    ╎ 1856  @Base\reinterpretarray.jl:497; _setindex_ra!
+ 1800╎    ╎    ╎    ╎    ╎  1800  @Base\array.jl:845; setindex!
+     ╎    ╎    ╎    ╎    ╎  56    @Base\refvalue.jl:56; getindex
+   56╎    ╎    ╎    ╎    ╎   56    @Base\Base.jl:33; getproperty
+
+==#
+
 @exclusive function export_wisdom(fname::AbstractString)
     f = ccall(:fopen, Ptr{Cvoid}, (Cstring,Cstring), fname, :w)
     systemerror("could not open wisdom file $fname for writing", f == C_NULL)
@@ -204,14 +335,14 @@ unsafe_set_timelimit(precision::fftwTypeSingle,seconds) =
 
 
 @static if fftw_provider == "mkl"
-    alignment_of(A::StridedArray{<:fftwDouble}) =
+    alignment_of(A::FFTStridedArray{<:fftwDouble}) =
         convert(Int32, convert(Int64, pointer(A)) % 16)
-    alignment_of(A::StridedArray{<:fftwSingle}) =
+    alignment_of(A::FFTStridedArray{<:fftwSingle}) =
         convert(Int32, convert(Int64, pointer(A)) % 16)
 else
-    alignment_of(A::StridedArray{T}) where {T<:fftwDouble} =
+    alignment_of(A::FFTStridedArray{T}) where {T<:fftwDouble} =
         ccall((:fftw_alignment_of, libfftw3[]), Int32, (Ptr{T},), A)
-    alignment_of(A::StridedArray{T}) where {T<:fftwSingle} =
+    alignment_of(A::FFTStridedArray{T}) where {T<:fftwSingle} =
         ccall((:fftwf_alignment_of, libfftw3f[]), Int32, (Ptr{T},), A)
 end
 
@@ -237,7 +368,7 @@ for P in (:cFFTWPlan, :rFFTWPlan, :r2rFFTWPlan) # complex, r2c/c2r, and r2r
             region::G # region (iterable) of dims that are transormed
             pinv::ScaledPlan
             function $P{T,K,inplace,N,G}(plan::PlanPtr, flags::Integer, R::G,
-                                         X::StridedArray{T,N}, Y::StridedArray) where {T<:fftwNumber,K,inplace,N,G}
+                                         X::FFTStridedArray{T,N}, Y::FFTStridedArray) where {T<:fftwNumber,K,inplace,N,G}
                 p = new(plan, size(X), size(Y), strides(X), strides(Y),
                         alignment_of(X), alignment_of(Y), flags, R)
                 finalizer(maybe_destroy_plan, p)
@@ -246,8 +377,8 @@ for P in (:cFFTWPlan, :rFFTWPlan, :r2rFFTWPlan) # complex, r2c/c2r, and r2r
         end
 
         function $P{T,K,inplace,N}(plan::PlanPtr, flags::Integer, R::G,
-                                   X::StridedArray{T,N},
-                                   Y::StridedArray) where {T<:fftwNumber,K,inplace,N,G}
+                                   X::FFTStridedArray{T,N},
+                                   Y::FFTStridedArray) where {T<:fftwNumber,K,inplace,N,G}
             $P{T,K,inplace,N,G}(plan, flags, R, X, Y)
         end
     end
@@ -421,7 +552,7 @@ 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
+function assert_applicable(p::FFTWPlan{T}, X::FFTStridedArray{T}) where T
     if size(X) != p.sz
         throw(ArgumentError("FFTW plan applied to wrong-size array"))
     elseif strides(X) != p.istride
@@ -431,7 +562,7 @@ function assert_applicable(p::FFTWPlan{T}, X::StridedArray{T}) where T
     end
 end
 
-function assert_applicable(p::FFTWPlan{T,K,inplace}, X::StridedArray{T}, Y::StridedArray) where {T,K,inplace}
+function assert_applicable(p::FFTWPlan{T,K,inplace}, X::FFTStridedArray{T}, Y::FFTStridedArray) where {T,K,inplace}
     assert_applicable(p, X)
     if size(Y) != p.osz
         throw(ArgumentError("FFTW plan applied to wrong-size output"))
@@ -463,42 +594,42 @@ unsafe_execute!(plan::FFTWPlan{<:fftwSingle}) =
     ccall((:fftwf_execute,libfftw3f[]), Cvoid, (PlanPtr,), plan)
 
 unsafe_execute!(plan::cFFTWPlan{T},
-                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} =
+                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwDouble} =
     ccall((:fftw_execute_dft,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::cFFTWPlan{T},
-                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} =
+                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwSingle} =
     ccall((:fftwf_execute_dft,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Float64,FORWARD},
-                X::StridedArray{Float64}, Y::StridedArray{Complex{Float64}}) =
+                X::FFTStridedArray{Float64}, Y::FFTStridedArray{Complex{Float64}}) =
     ccall((:fftw_execute_dft_r2c,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{Float64},Ptr{Complex{Float64}}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Float32,FORWARD},
-                X::StridedArray{Float32}, Y::StridedArray{Complex{Float32}}) =
+                X::FFTStridedArray{Float32}, Y::FFTStridedArray{Complex{Float32}}) =
     ccall((:fftwf_execute_dft_r2c,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{Float32},Ptr{Complex{Float32}}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Complex{Float64},BACKWARD},
-                X::StridedArray{Complex{Float64}}, Y::StridedArray{Float64}) =
+                X::FFTStridedArray{Complex{Float64}}, Y::FFTStridedArray{Float64}) =
     ccall((:fftw_execute_dft_c2r,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{Complex{Float64}},Ptr{Float64}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Complex{Float32},BACKWARD},
-                X::StridedArray{Complex{Float32}}, Y::StridedArray{Float32}) =
+                X::FFTStridedArray{Complex{Float32}}, Y::FFTStridedArray{Float32}) =
     ccall((:fftwf_execute_dft_c2r,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{Complex{Float32}},Ptr{Float32}), plan, X, Y)
 
 unsafe_execute!(plan::r2rFFTWPlan{T},
-                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} =
+                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwDouble} =
     ccall((:fftw_execute_r2r,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::r2rFFTWPlan{T},
-                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} =
+                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwSingle} =
     ccall((:fftwf_execute_r2r,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
@@ -513,9 +644,15 @@ unsafe_execute!(plan::r2rFFTWPlan{T},
 #    re-use the table of trigonometric constants from the first plan.
 
 # Compute dims and howmany for FFTW guru planner
-function dims_howmany(X::StridedArray, Y::StridedArray,
+function dims_howmany(X::FFTStridedArray, Y::FFTStridedArray,
                       sz::Array{Int,1}, region)
-    reg = Int[region...]
+    reg = if isa(region, Int)
+        region:region
+    elseif isa(region, Tuple)
+        Int[region...]
+    else
+        Int.(region)
+    end
     if length(unique(reg)) < length(reg)
         throw(ArgumentError("each dimension can be transformed at most once"))
     end
@@ -558,8 +695,8 @@ end
 
 for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
                          (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f))
-    @eval @exclusive function cFFTWPlan{$Tc,K,inplace,N}(X::StridedArray{$Tc,N},
-                                              Y::StridedArray{$Tc,N},
+    @eval @exclusive function cFFTWPlan{$Tc,K,inplace,N}(X::FFTStridedArray{$Tc,N},
+                                              Y::FFTStridedArray{$Tc,N},
                                               region, flags::Integer, timelimit::Real) where {K,inplace,N}
         direction = K
         unsafe_set_timelimit($Tr, timelimit)
@@ -578,11 +715,11 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return cFFTWPlan{$Tc,K,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function rFFTWPlan{$Tr,$FORWARD,inplace,N}(X::StridedArray{$Tr,N},
-                                                     Y::StridedArray{$Tc,N},
+    @eval @exclusive function rFFTWPlan{$Tr,$FORWARD,inplace,N}(X::FFTStridedArray{$Tr,N},
+                                                     Y::FFTStridedArray{$Tc,N},
                                                      region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
-        region = circshift(Int[region...],-1) # FFTW halves last dim
+        region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
         unsafe_set_timelimit($Tr, timelimit)
         dims, howmany = dims_howmany(X, Y, [size(X)...], region)
         plan = ccall(($(string(fftw,"_plan_guru64_dft_r2c")),$lib[]),
@@ -598,11 +735,11 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return rFFTWPlan{$Tr,$FORWARD,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function rFFTWPlan{$Tc,$BACKWARD,inplace,N}(X::StridedArray{$Tc,N},
-                                                      Y::StridedArray{$Tr,N},
+    @eval @exclusive function rFFTWPlan{$Tc,$BACKWARD,inplace,N}(X::FFTStridedArray{$Tc,N},
+                                                      Y::FFTStridedArray{$Tr,N},
                                                       region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
-        region = circshift(Int[region...],-1) # FFTW halves last dim
+        region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
         unsafe_set_timelimit($Tr, timelimit)
         dims, howmany = dims_howmany(X, Y, [size(Y)...], region)
         plan = ccall(($(string(fftw,"_plan_guru64_dft_c2r")),$lib[]),
@@ -618,8 +755,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return rFFTWPlan{$Tc,$BACKWARD,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function r2rFFTWPlan{$Tr,Any,inplace,N}(X::StridedArray{$Tr,N},
-                                                  Y::StridedArray{$Tr,N},
+    @eval @exclusive function r2rFFTWPlan{$Tr,Any,inplace,N}(X::FFTStridedArray{$Tr,N},
+                                                  Y::FFTStridedArray{$Tr,N},
                                                   region, kinds, flags::Integer,
                                                   timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
@@ -640,8 +777,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
     end
 
     # support r2r transforms of complex = transforms of real & imag parts
-    @eval @exclusive function r2rFFTWPlan{$Tc,Any,inplace,N}(X::StridedArray{$Tc,N},
-                                                  Y::StridedArray{$Tc,N},
+    @eval @exclusive function r2rFFTWPlan{$Tc,Any,inplace,N}(X::FFTStridedArray{$Tc,N},
+                                                  Y::FFTStridedArray{$Tc,N},
                                                   region, kinds, flags::Integer,
                                                   timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
@@ -668,12 +805,12 @@ end
 
 # Convert arrays of numeric types to FFTW-supported packed complex-float types
 # (FIXME: is there a way to use the Julia promotion rules more cleverly here?)
-fftwcomplex(X::StridedArray{<:fftwComplex}) = X
+fftwcomplex(X::FFTStridedArray{<:fftwComplex}) = X
 fftwcomplex(X::AbstractArray{T}) where {T<:fftwReal} =
     copyto!(Array{typeof(complex(zero(T)))}(undef, size(X)), X)
 fftwcomplex(X::AbstractArray{<:Real}) = copyto!(Array{Complex{Float64}}(undef, size(X)),X)
 fftwcomplex(X::AbstractArray{<:Complex}) = copyto!(Array{Complex{Float64}}(undef, size(X)), X)
-fftwfloat(X::StridedArray{<:fftwReal}) = X
+fftwfloat(X::FFTStridedArray{<:fftwReal}) = X
 fftwfloat(X::AbstractArray{<:Real}) = copyto!(Array{Float64}(undef, size(X)), X)
 fftwfloat(X::AbstractArray{<:Complex}) = fftwcomplex(X)
 
@@ -682,21 +819,21 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD))
     plan_f! = Symbol("plan_",f,"!")
     idirection = -direction
     @eval begin
-        function $plan_f(X::StridedArray{T,N}, region;
+        function $plan_f(X::FFTStridedArray{T,N}, region;
                          flags::Integer=ESTIMATE,
                          timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N}
             cFFTWPlan{T,$direction,false,N}(X, fakesimilar(flags, X, T),
                                             region, flags, timelimit)
         end
 
-        function $plan_f!(X::StridedArray{T,N}, region;
+        function $plan_f!(X::FFTStridedArray{T,N}, region;
                          flags::Integer=ESTIMATE,
                          timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N}
             cFFTWPlan{T,$direction,true,N}(X, X, region, flags, timelimit)
         end
-        $plan_f(X::StridedArray{<:fftwComplex}; kws...) =
+        $plan_f(X::FFTStridedArray{<:fftwComplex}; kws...) =
             $plan_f(X, 1:ndims(X); kws...)
-        $plan_f!(X::StridedArray{<:fftwComplex}; kws...) =
+        $plan_f!(X::FFTStridedArray{<:fftwComplex}; kws...) =
             $plan_f!(X, 1:ndims(X); kws...)
 
         function plan_inv(p::cFFTWPlan{T,$direction,inplace,N}) where {T<:fftwComplex,N,inplace}
@@ -709,20 +846,20 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD))
     end
 end
 
-function mul!(y::StridedArray{T}, p::cFFTWPlan{T}, x::StridedArray{T}) where T
+function mul!(y::FFTStridedArray{T}, p::cFFTWPlan{T}, x::FFTStridedArray{T}) where T
     assert_applicable(p, x, y)
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::cFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N}
+function *(p::cFFTWPlan{T,K,false}, x::FFTStridedArray{T,N}) where {T,K,N}
     assert_applicable(p, x)
     y = Array{T}(undef, p.osz)::Array{T,N}
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::cFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
+function *(p::cFFTWPlan{T,K,true}, x::FFTStridedArray{T}) where {T,K}
     assert_applicable(p, x)
     unsafe_execute!(p, x, x)
     return x
@@ -730,10 +867,21 @@ end
 
 # rfft/brfft and planned variants.  No in-place version for now.
 
+for f in (:brfft!, :irfft!)
+    pf = Symbol("plan_", f)
+    @eval begin
+        $f(x::AbstractArray, d::Integer) = (x = PaddedRFFTArray(x); $pf(x, d) * x)
+        $f(x::AbstractArray, d::Integer, region) = (x = PaddedRFFTArray(x); $pf(x, d, region) * x)
+        $pf(x::PaddedRFFTArray, d::Integer; kws...) = $pf(x, d, 1:ndims(x);kws...)
+        $f(x::AbstractArray{<:Real}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
+        $f(x::AbstractArray{<:Complex{<:Union{Integer,Rational}}}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
+    end
+end
+
 for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
     # Note: use $FORWARD and $BACKWARD below because of issue #9775
     @eval begin
-        function plan_rfft(X::StridedArray{$Tr,N}, region;
+        function plan_rfft(X::FFTStridedArray{$Tr,N}, region;
                            flags::Integer=ESTIMATE,
                            timelimit::Real=NO_TIMELIMIT) where N
             osize = rfft_output_size(X, region)
@@ -741,7 +889,7 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             rFFTWPlan{$Tr,$FORWARD,false,N}(X, Y, region, flags, timelimit)
         end
 
-        function plan_brfft(X::StridedArray{$Tc,N}, d::Integer, region;
+        function plan_brfft(X::FFTStridedArray{$Tc,N}, d::Integer, region;
                             flags::Integer=ESTIMATE,
                             timelimit::Real=NO_TIMELIMIT) where N
             osize = brfft_output_size(X, d, region)
@@ -760,8 +908,22 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             end
         end
 
-        plan_rfft(X::StridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...)
-        plan_brfft(X::StridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...)
+        plan_rfft(X::FFTStridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...)
+        plan_brfft(X::FFTStridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...)
+
+        function plan_brfft!(X::PaddedRFFTArray{$Tr,N}, ::Integer, region;
+            flags::Integer=ESTIMATE,
+            timelimit::Real=NO_TIMELIMIT) where N
+            return rFFTWPlan{$Tc,$BACKWARD,true,N}(complex_view(X), real_view(X), region, flags, timelimit)
+        end
+
+        plan_brfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_brfft!(X,1:ndims(X);kws...)
+
+        function plan_irfft!(x::PaddedRFFTArray{$Tr,N}, d::Integer, region; kws...) where N
+            ScaledPlan(plan_brfft!(x, d, region; kws...), normalization($Tr, size(real_view(x)), region))
+        end
+
+        plan_irfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_irfft!(X,1:ndims(X);kws...)
 
         function plan_inv(p::rFFTWPlan{$Tr,$FORWARD,false,N}) where N
             X = Array{$Tr}(undef, p.sz)
@@ -781,25 +943,25 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
                        normalization(Y, p.region))
         end
 
-        function mul!(y::StridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::StridedArray{$Tr})
+        function mul!(y::FFTStridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::FFTStridedArray{$Tr})
             assert_applicable(p, x, y)
             unsafe_execute!(p, x, y)
             return y
         end
-        function mul!(y::StridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::StridedArray{$Tc})
+        function mul!(y::FFTStridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::FFTStridedArray{$Tc})
             assert_applicable(p, x, y)
             unsafe_execute!(p, x, y) # note: may overwrite x as well as y!
             return y
         end
 
-        function *(p::rFFTWPlan{$Tr,$FORWARD,false}, x::StridedArray{$Tr,N}) where N
+        function *(p::rFFTWPlan{$Tr,$FORWARD,false}, x::FFTStridedArray{$Tr,N}) where N
             assert_applicable(p, x)
             y = Array{$Tc}(undef, p.osz)
             unsafe_execute!(p, x, y)
             return y
         end
 
-        function *(p::rFFTWPlan{$Tc,$BACKWARD,false}, x::StridedArray{$Tc,N}) where N
+        function *(p::rFFTWPlan{$Tc,$BACKWARD,false}, x::FFTStridedArray{$Tc,N}) where N
             if p.flags & PRESERVE_INPUT != 0
                 assert_applicable(p, x)
                 y = Array{$Tr}(undef, p.osz)
@@ -812,9 +974,18 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             end
             return y
         end
+
+        *(p::rFFTWPlan{$Tc,$BACKWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N = 
+            (mul!(real_view(f), p, complex_view(f)); real_view(f))
     end
 end
 
+*(p::ScaledPlan, f::PaddedRFFTArray) = begin
+    p.p * f
+    rmul!(real_view(f), p.scale)
+    real_view(f)
+end
+
 # FFTW r2r transforms (low-level interface)
 
 for f in (:r2r, :r2r!)
@@ -830,14 +1001,14 @@ for f in (:r2r, :r2r!)
     end
 end
 
-function plan_r2r(X::StridedArray{T,N}, kinds, region;
+function plan_r2r(X::FFTStridedArray{T,N}, kinds, region;
                   flags::Integer=ESTIMATE,
                   timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N}
     r2rFFTWPlan{T,Any,false,N}(X, fakesimilar(flags, X, T), region, kinds,
                                flags, timelimit)
 end
 
-function plan_r2r!(X::StridedArray{T,N}, kinds, region;
+function plan_r2r!(X::FFTStridedArray{T,N}, kinds, region;
                    flags::Integer=ESTIMATE,
                    timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N}
     r2rFFTWPlan{T,Any,true,N}(X, X, region, kinds, flags, timelimit)
@@ -872,21 +1043,22 @@ function plan_inv(p::r2rFFTWPlan{T,K,inplace,N}) where {T<:fftwNumber,K,inplace,
                              1:length(iK)))
 end
 
-function mul!(y::StridedArray{T}, p::r2rFFTWPlan{T}, x::StridedArray{T}) where T
+function mul!(y::FFTStridedArray{T}, p::r2rFFTWPlan{T}, x::FFTStridedArray{T}) where T
     assert_applicable(p, x, y)
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::r2rFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N}
+function *(p::r2rFFTWPlan{T,K,false}, x::FFTStridedArray{T,N}) where {T,K,N}
     assert_applicable(p, x)
     y = Array{T}(undef, p.osz)::Array{T,N}
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::r2rFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
+function *(p::r2rFFTWPlan{T,K,true}, x::FFTStridedArray{T}) where {T,K}
     assert_applicable(p, x)
     unsafe_execute!(p, x, x)
     return x
 end
+

From 837a333e4b043cf7ae7023f628712e3d85a83db5 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Tue, 28 Sep 2021 22:43:45 -0400
Subject: [PATCH 2/4] Updates

---
 src/FFTW.jl      |   2 +-
 src/fft.jl       | 264 +++++++++++++++++++----------------------------
 test/runtests.jl |  48 +++++++++
 3 files changed, 156 insertions(+), 158 deletions(-)

diff --git a/src/FFTW.jl b/src/FFTW.jl
index 0493cbc..2274354 100644
--- a/src/FFTW.jl
+++ b/src/FFTW.jl
@@ -12,7 +12,7 @@ import AbstractFFTs: Plan, ScaledPlan,
                      rfft_output_size, brfft_output_size,
                      plan_inv, normalization
 
-export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!, brfft!, irfft!, plan_brfft!, plan_irfft!
+export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!, brfft!, irfft!, plan_brfft!, plan_irfft!, rfft!, plan_rfft!, PaddedRFFTArray
 
 include("providers.jl")
 
diff --git a/src/fft.jl b/src/fft.jl
index 94929b6..ad2dde0 100644
--- a/src/fft.jl
+++ b/src/fft.jl
@@ -1,6 +1,6 @@
 # This file was formerly a part of Julia. License is MIT: https://julialang.org/license
 
-import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, copy, similar, eltype, \
+import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, copy, similar, eltype, \, getindex
 import LinearAlgebra: mul!
 
 
@@ -115,12 +115,14 @@ struct PaddedRFFTArray{T, N, A, Ar, Ac} <: AbstractArray{T, N} where {A <: Union
     r::Ar # real view
     c::Ac # complex view
 
+    # wrap existing complex array
     function PaddedRFFTArray(data::StridedArray{Tc, N}) where {Tc <: fftwComplex, N}
         fsize = 2 * (size(data, 1) - 1)
         r = view(reinterpret(real(Tc), data), Base.OneTo(fsize), ntuple(i -> Colon(), Val(ndims(data) - 1))...)
         new{real(Tc), N, typeof(data), typeof(r), typeof(data)}(data, r, data)
     end
 
+    # wrap existing padded reals array
     function PaddedRFFTArray(data::StridedArray{T, N}, nx) where {T <: fftwReal, N}
         fsize = size(data, 1)
         iseven(fsize) || throw(
@@ -131,52 +133,34 @@ struct PaddedRFFTArray{T, N, A, Ar, Ac} <: AbstractArray{T, N} where {A <: Union
         r = view(data, Base.OneTo(nx), ntuple(i->Colon(), Val(ndims(data) - 1))...)
         new{T, N, typeof(data), typeof(r), typeof(c)}(data, r, c)
     end
-end
 
-FFTStridedArray{T, N} = Union{StridedArray{T, N}, PaddedRFFTArray{T, N}}
+    # copy reals array to padded array
+    function PaddedRFFTArray(reals::StridedArray{T, N}) where {T <: fftwReal, N}
+        fsize = (size(reals, 1) ÷ 2 + 1) * 2
+        data = Array{T, N}(undef, (fsize, size(reals)[2:end]...))
+        c = reinterpret(Complex{T}, data)
+        r = view(data, Base.OneTo(size(reals, 1)), ntuple(i->Colon(), Val(ndims(reals) - 1))...)
+        @inbounds copyto!(r, reals)
+        new{T, N, typeof(data), typeof(r), typeof(c)}(data, r, c)
+    end
+
+    # copy existing padded array
+    function PaddedRFFTArray(a::PaddedRFFTArray{T, N, A, Ar, Ac}) where {T, N, A, Ar, Ac}
+        data = copy(parent(a))
+        c = reinterpret(Complex{T}, data)
+        r = view(data, Base.OneTo(size(real_view(a), 1)), ntuple(i->Colon(), Val(ndims(data) - 1))...)
+        new{T, N, A, Ar, Ac}(data, r, c)
+    end
+end
 
 @inline real_view(S::PaddedRFFTArray) = S.r
 
 @inline complex_view(S::PaddedRFFTArray) = S.c
 
-copy(S::PaddedRFFTArray) = PaddedRFFTArray(copy(parent(S)), size(real_view(S), 1))
-
-# similar(f::PaddedRFFTArray, ::Type{T}, dims::Tuple{Vararg{Int, N}}) where {T <: fftwReal, N} =
-#     PaddedRFFTArray{T, N}(dims)
-# similar(f::PaddedRFFTArray{T}, dims::Vararg{Tuple{Int, N}}) where {T, N} =
-#     PaddedRFFTArray{T, N}(dims)
-# similar(f::PaddedRFFTArray, ::Type{T}) where {T <: fftwReal} =
-#     PaddedRFFTArray{T}(size(real_view(f)))
-# similar(f::PaddedRFFTArray{T, N}) where {T,N} = 
-#     PaddedRFFTArray{T, N}(similar(parent(f)), size(real_view(f), 1)) 
-
-# size(S::PaddedRFFTArray) =
-#     size(complex_view(S))
-
-# IndexStyle(::Type{T}) where {T<:PaddedRFFTArray} = 
-#     IndexLinear()
-
-function PaddedRFFTArray{T}(ndims::Vararg{Integer, N}) where {T <: fftwReal, N}
-    fsize = (ndims[1] ÷ 2 + 1) * 2
-    a = zeros(T, (fsize, ndims[2:end]...))
-    PaddedRFFTArray{T, N}(a, ndims[1])
-end
-
-PaddedRFFTArray{T}(ndims::NTuple{N, Integer}) where {T <: fftwReal, 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 copyto!(t.r, a)
-    return t
-end
-
+size(a::PaddedRFFTArray) = size(complex_view(a))
+getindex(a::PaddedRFFTArray, args...) = getindex(complex_view(a), args...)
+copy(a::PaddedRFFTArray) = PaddedRFFTArray(a)
+parent(a::PaddedRFFTArray) = a.parent
 
 # For ESTIMATE plans, FFTW allows one to pass NULL for the array pointer,
 # since it is not written to.  Hence, it is convenient to create an
@@ -208,67 +192,6 @@ alignment_of(A::FakeArray) = Int32(0)
 # around FFTW's internal file i/o buffering [see the BUFSZ constant in
 # FFTW's api/import-wisdom-from-file.c file].
 
-#==
-a = rand(10002, 80000)
-const cccc = rfft(a[1:10000, :])
-function e()
-       d = copy(cccc)
-       irfft!(d, 10000)
-       end
-
-Overhead ╎ [+additional indent] Count File:Line; Function
-=========================================================
-     ╎40027 @Base\task.jl:411; (::VSCodeServer.var"#53#54")()
-     ╎ 40027 @VSCodeServer\src\eval.jl:34; macro expansion
-     ╎  40027 @Base\essentials.jl:706; invokelatest(::Any)
-     ╎   40027 @Base\essentials.jl:708; #invokelatest#2
-     ╎    40027 @VSCodeServer\src\repl.jl:124; (::VSCodeServer.var"#68#70"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt})()
-     ╎     40027 @Base\logging.jl:603; with_logger
-     ╎    ╎ 40027 @Base\logging.jl:491; with_logstate(f::Function, logstate::Any)
-     ╎    ╎  40027 @VSCodeServer\src\repl.jl:123; (::VSCodeServer.var"#69#71"{Module, Expr, REPL.LineEditREPL, REPL.LineEdit.Prompt})()
-     ╎    ╎   40027 @VSCodeServer\src\repl.jl:157; repleval(m::Module, code::Expr, #unused#::String)
-     ╎    ╎    40027 @Base\Base.jl:39; eval
-     ╎    ╎     40027 @Base\boot.jl:360; eval
-     ╎    ╎    ╎ 931   REPL[40]:2; e()
-  931╎    ╎    ╎  931   @Base\array.jl:349; copy
-     ╎    ╎    ╎ 39087 REPL[40]:3; e()
-     ╎    ╎    ╎  39087 @FFTW\src\fft.jl:812; irfft!(x::Matrix{ComplexF64}, d::Int64)
-     ╎    ╎    ╎   35687 @FFTW\src\fft.jl:923; *
-     ╎    ╎    ╎    35687 @FFTW\src\fft.jl:917; *
-     ╎    ╎    ╎     35687 @FFTW\src\fft.jl:892; mul!
-35687╎    ╎    ╎    ╎ 35687 @FFTW\src\fft.jl:556; unsafe_execute!(plan::FFTW.rFFTWPlan{ComplexF64, 1, true, 2, UnitRange{Int64}}, X::Matrix{ComplexF64}, Y::SubArray{Float64, 2...
-     ╎    ╎    ╎   3400  @FFTW\src\fft.jl:924; *
-     ╎    ╎    ╎    3400  ...worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\generic.jl:182; rmul!(X::SubArray{Float64, 2, Base.ReinterpretArray{Float64, 2, ComplexF64, Matrix{ComplexF64}, false}, Tuple{Base.OneTo{Int64...
-     ╎    ╎    ╎     257   @Base\simdloop.jl:75; macro expansion
-  257╎    ╎    ╎    ╎ 257   @Base\int.jl:83; <
-     ╎    ╎    ╎     3143  @Base\simdloop.jl:77; macro expansion
-     ╎    ╎    ╎    ╎ 3143  ...orker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\generic.jl:183; macro expansion
-     ╎    ╎    ╎    ╎  731   @Base\abstractarray.jl:1170; getindex
-     ╎    ╎    ╎    ╎   731   @Base\abstractarray.jl:1214; _getindex
-     ╎    ╎    ╎    ╎    731   @Base\subarray.jl:276; getindex
-     ╎    ╎    ╎    ╎     731   @Base\reinterpretarray.jl:320; getindex
-     ╎    ╎    ╎    ╎    ╎ 368   @Base\reinterpretarray.jl:355; _getindex_ra
-     ╎    ╎    ╎    ╎    ╎  368   @Base\div.jl:120; divrem
-     ╎    ╎    ╎    ╎    ╎   368   @Base\div.jl:124; divrem
-  368╎    ╎    ╎    ╎    ╎    368   @Base\int.jl:262; rem
-     ╎    ╎    ╎    ╎    ╎ 363   @Base\reinterpretarray.jl:371; _getindex_ra
-  259╎    ╎    ╎    ╎    ╎  259   @Base\array.jl:806; getindex
-     ╎    ╎    ╎    ╎    ╎  104   @Base\refvalue.jl:57; setindex!
-  104╎    ╎    ╎    ╎    ╎   104   @Base\Base.jl:34; setproperty!
-     ╎    ╎    ╎    ╎  2404  @Base\abstractarray.jl:1267; setindex!
-     ╎    ╎    ╎    ╎   2404  @Base\abstractarray.jl:1297; _setindex!
-     ╎    ╎    ╎    ╎    2404  @Base\subarray.jl:317; setindex!
-     ╎    ╎    ╎    ╎     2404  @Base\reinterpretarray.jl:446; setindex!
-     ╎    ╎    ╎    ╎    ╎ 548   @Base\reinterpretarray.jl:495; _setindex_ra!
-     ╎    ╎    ╎    ╎    ╎  548   @Base\refvalue.jl:57; setindex!
-  548╎    ╎    ╎    ╎    ╎   548   @Base\Base.jl:34; setproperty!
-     ╎    ╎    ╎    ╎    ╎ 1856  @Base\reinterpretarray.jl:497; _setindex_ra!
- 1800╎    ╎    ╎    ╎    ╎  1800  @Base\array.jl:845; setindex!
-     ╎    ╎    ╎    ╎    ╎  56    @Base\refvalue.jl:56; getindex
-   56╎    ╎    ╎    ╎    ╎   56    @Base\Base.jl:33; getproperty
-
-==#
-
 @exclusive function export_wisdom(fname::AbstractString)
     f = ccall(:fopen, Ptr{Cvoid}, (Cstring,Cstring), fname, :w)
     systemerror("could not open wisdom file $fname for writing", f == C_NULL)
@@ -335,14 +258,14 @@ unsafe_set_timelimit(precision::fftwTypeSingle,seconds) =
 
 
 @static if fftw_provider == "mkl"
-    alignment_of(A::FFTStridedArray{<:fftwDouble}) =
+    alignment_of(A::StridedArray{<:fftwDouble}) =
         convert(Int32, convert(Int64, pointer(A)) % 16)
-    alignment_of(A::FFTStridedArray{<:fftwSingle}) =
+    alignment_of(A::StridedArray{<:fftwSingle}) =
         convert(Int32, convert(Int64, pointer(A)) % 16)
 else
-    alignment_of(A::FFTStridedArray{T}) where {T<:fftwDouble} =
+    alignment_of(A::StridedArray{T}) where {T<:fftwDouble} =
         ccall((:fftw_alignment_of, libfftw3[]), Int32, (Ptr{T},), A)
-    alignment_of(A::FFTStridedArray{T}) where {T<:fftwSingle} =
+    alignment_of(A::StridedArray{T}) where {T<:fftwSingle} =
         ccall((:fftwf_alignment_of, libfftw3f[]), Int32, (Ptr{T},), A)
 end
 
@@ -368,7 +291,7 @@ for P in (:cFFTWPlan, :rFFTWPlan, :r2rFFTWPlan) # complex, r2c/c2r, and r2r
             region::G # region (iterable) of dims that are transormed
             pinv::ScaledPlan
             function $P{T,K,inplace,N,G}(plan::PlanPtr, flags::Integer, R::G,
-                                         X::FFTStridedArray{T,N}, Y::FFTStridedArray) where {T<:fftwNumber,K,inplace,N,G}
+                                         X::StridedArray{T,N}, Y::StridedArray) where {T<:fftwNumber,K,inplace,N,G}
                 p = new(plan, size(X), size(Y), strides(X), strides(Y),
                         alignment_of(X), alignment_of(Y), flags, R)
                 finalizer(maybe_destroy_plan, p)
@@ -377,8 +300,8 @@ for P in (:cFFTWPlan, :rFFTWPlan, :r2rFFTWPlan) # complex, r2c/c2r, and r2r
         end
 
         function $P{T,K,inplace,N}(plan::PlanPtr, flags::Integer, R::G,
-                                   X::FFTStridedArray{T,N},
-                                   Y::FFTStridedArray) where {T<:fftwNumber,K,inplace,N,G}
+                                   X::StridedArray{T,N},
+                                   Y::StridedArray) where {T<:fftwNumber,K,inplace,N,G}
             $P{T,K,inplace,N,G}(plan, flags, R, X, Y)
         end
     end
@@ -552,7 +475,7 @@ 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::FFTStridedArray{T}) where T
+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
@@ -562,7 +485,7 @@ function assert_applicable(p::FFTWPlan{T}, X::FFTStridedArray{T}) where T
     end
 end
 
-function assert_applicable(p::FFTWPlan{T,K,inplace}, X::FFTStridedArray{T}, Y::FFTStridedArray) where {T,K,inplace}
+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"))
@@ -594,42 +517,42 @@ unsafe_execute!(plan::FFTWPlan{<:fftwSingle}) =
     ccall((:fftwf_execute,libfftw3f[]), Cvoid, (PlanPtr,), plan)
 
 unsafe_execute!(plan::cFFTWPlan{T},
-                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwDouble} =
+                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} =
     ccall((:fftw_execute_dft,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::cFFTWPlan{T},
-                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwSingle} =
+                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} =
     ccall((:fftwf_execute_dft,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Float64,FORWARD},
-                X::FFTStridedArray{Float64}, Y::FFTStridedArray{Complex{Float64}}) =
+                X::StridedArray{Float64}, Y::StridedArray{Complex{Float64}}) =
     ccall((:fftw_execute_dft_r2c,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{Float64},Ptr{Complex{Float64}}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Float32,FORWARD},
-                X::FFTStridedArray{Float32}, Y::FFTStridedArray{Complex{Float32}}) =
+                X::StridedArray{Float32}, Y::StridedArray{Complex{Float32}}) =
     ccall((:fftwf_execute_dft_r2c,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{Float32},Ptr{Complex{Float32}}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Complex{Float64},BACKWARD},
-                X::FFTStridedArray{Complex{Float64}}, Y::FFTStridedArray{Float64}) =
+                X::StridedArray{Complex{Float64}}, Y::StridedArray{Float64}) =
     ccall((:fftw_execute_dft_c2r,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{Complex{Float64}},Ptr{Float64}), plan, X, Y)
 
 unsafe_execute!(plan::rFFTWPlan{Complex{Float32},BACKWARD},
-                X::FFTStridedArray{Complex{Float32}}, Y::FFTStridedArray{Float32}) =
+                X::StridedArray{Complex{Float32}}, Y::StridedArray{Float32}) =
     ccall((:fftwf_execute_dft_c2r,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{Complex{Float32}},Ptr{Float32}), plan, X, Y)
 
 unsafe_execute!(plan::r2rFFTWPlan{T},
-                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwDouble} =
+                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwDouble} =
     ccall((:fftw_execute_r2r,libfftw3[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
 unsafe_execute!(plan::r2rFFTWPlan{T},
-                X::FFTStridedArray{T}, Y::FFTStridedArray{T}) where {T<:fftwSingle} =
+                X::StridedArray{T}, Y::StridedArray{T}) where {T<:fftwSingle} =
     ccall((:fftwf_execute_r2r,libfftw3f[]), Cvoid,
           (PlanPtr,Ptr{T},Ptr{T}), plan, X, Y)
 
@@ -644,7 +567,7 @@ unsafe_execute!(plan::r2rFFTWPlan{T},
 #    re-use the table of trigonometric constants from the first plan.
 
 # Compute dims and howmany for FFTW guru planner
-function dims_howmany(X::FFTStridedArray, Y::FFTStridedArray,
+function dims_howmany(X::StridedArray, Y::StridedArray,
                       sz::Array{Int,1}, region)
     reg = if isa(region, Int)
         region:region
@@ -695,8 +618,8 @@ end
 
 for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
                          (:Float32,:(Complex{Float32}),"fftwf",:libfftw3f))
-    @eval @exclusive function cFFTWPlan{$Tc,K,inplace,N}(X::FFTStridedArray{$Tc,N},
-                                              Y::FFTStridedArray{$Tc,N},
+    @eval @exclusive function cFFTWPlan{$Tc,K,inplace,N}(X::StridedArray{$Tc,N},
+                                              Y::StridedArray{$Tc,N},
                                               region, flags::Integer, timelimit::Real) where {K,inplace,N}
         direction = K
         unsafe_set_timelimit($Tr, timelimit)
@@ -715,8 +638,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return cFFTWPlan{$Tc,K,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function rFFTWPlan{$Tr,$FORWARD,inplace,N}(X::FFTStridedArray{$Tr,N},
-                                                     Y::FFTStridedArray{$Tc,N},
+    @eval @exclusive function rFFTWPlan{$Tr,$FORWARD,inplace,N}(X::StridedArray{$Tr,N},
+                                                     Y::StridedArray{$Tc,N},
                                                      region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
         region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
@@ -735,8 +658,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return rFFTWPlan{$Tr,$FORWARD,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function rFFTWPlan{$Tc,$BACKWARD,inplace,N}(X::FFTStridedArray{$Tc,N},
-                                                      Y::FFTStridedArray{$Tr,N},
+    @eval @exclusive function rFFTWPlan{$Tc,$BACKWARD,inplace,N}(X::StridedArray{$Tc,N},
+                                                      Y::StridedArray{$Tr,N},
                                                       region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
         region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
@@ -755,8 +678,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
         return rFFTWPlan{$Tc,$BACKWARD,inplace,N}(plan, flags, R, X, Y)
     end
 
-    @eval @exclusive function r2rFFTWPlan{$Tr,Any,inplace,N}(X::FFTStridedArray{$Tr,N},
-                                                  Y::FFTStridedArray{$Tr,N},
+    @eval @exclusive function r2rFFTWPlan{$Tr,Any,inplace,N}(X::StridedArray{$Tr,N},
+                                                  Y::StridedArray{$Tr,N},
                                                   region, kinds, flags::Integer,
                                                   timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
@@ -777,8 +700,8 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
     end
 
     # support r2r transforms of complex = transforms of real & imag parts
-    @eval @exclusive function r2rFFTWPlan{$Tc,Any,inplace,N}(X::FFTStridedArray{$Tc,N},
-                                                  Y::FFTStridedArray{$Tc,N},
+    @eval @exclusive function r2rFFTWPlan{$Tc,Any,inplace,N}(X::StridedArray{$Tc,N},
+                                                  Y::StridedArray{$Tc,N},
                                                   region, kinds, flags::Integer,
                                                   timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
@@ -805,12 +728,12 @@ end
 
 # Convert arrays of numeric types to FFTW-supported packed complex-float types
 # (FIXME: is there a way to use the Julia promotion rules more cleverly here?)
-fftwcomplex(X::FFTStridedArray{<:fftwComplex}) = X
+fftwcomplex(X::StridedArray{<:fftwComplex}) = X
 fftwcomplex(X::AbstractArray{T}) where {T<:fftwReal} =
     copyto!(Array{typeof(complex(zero(T)))}(undef, size(X)), X)
 fftwcomplex(X::AbstractArray{<:Real}) = copyto!(Array{Complex{Float64}}(undef, size(X)),X)
 fftwcomplex(X::AbstractArray{<:Complex}) = copyto!(Array{Complex{Float64}}(undef, size(X)), X)
-fftwfloat(X::FFTStridedArray{<:fftwReal}) = X
+fftwfloat(X::StridedArray{<:fftwReal}) = X
 fftwfloat(X::AbstractArray{<:Real}) = copyto!(Array{Float64}(undef, size(X)), X)
 fftwfloat(X::AbstractArray{<:Complex}) = fftwcomplex(X)
 
@@ -819,21 +742,21 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD))
     plan_f! = Symbol("plan_",f,"!")
     idirection = -direction
     @eval begin
-        function $plan_f(X::FFTStridedArray{T,N}, region;
+        function $plan_f(X::StridedArray{T,N}, region;
                          flags::Integer=ESTIMATE,
                          timelimit::Real=NO_TIMELIMIT) where {T<:fftwComplex,N}
             cFFTWPlan{T,$direction,false,N}(X, fakesimilar(flags, X, T),
                                             region, flags, timelimit)
         end
 
-        function $plan_f!(X::FFTStridedArray{T,N}, region;
+        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)
         end
-        $plan_f(X::FFTStridedArray{<:fftwComplex}; kws...) =
+        $plan_f(X::StridedArray{<:fftwComplex}; kws...) =
             $plan_f(X, 1:ndims(X); kws...)
-        $plan_f!(X::FFTStridedArray{<:fftwComplex}; kws...) =
+        $plan_f!(X::StridedArray{<:fftwComplex}; kws...) =
             $plan_f!(X, 1:ndims(X); kws...)
 
         function plan_inv(p::cFFTWPlan{T,$direction,inplace,N}) where {T<:fftwComplex,N,inplace}
@@ -846,20 +769,20 @@ for (f,direction) in ((:fft,FORWARD), (:bfft,BACKWARD))
     end
 end
 
-function mul!(y::FFTStridedArray{T}, p::cFFTWPlan{T}, x::FFTStridedArray{T}) where T
+function mul!(y::StridedArray{T}, p::cFFTWPlan{T}, x::StridedArray{T}) where T
     assert_applicable(p, x, y)
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::cFFTWPlan{T,K,false}, x::FFTStridedArray{T,N}) where {T,K,N}
+function *(p::cFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N}
     assert_applicable(p, x)
     y = Array{T}(undef, p.osz)::Array{T,N}
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::cFFTWPlan{T,K,true}, x::FFTStridedArray{T}) where {T,K}
+function *(p::cFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
     assert_applicable(p, x)
     unsafe_execute!(p, x, x)
     return x
@@ -867,21 +790,32 @@ end
 
 # rfft/brfft and planned variants.  No in-place version for now.
 
+for f in (:rfft!,)
+    pf = Symbol("plan_", f)
+    @eval begin
+        $f(x::AbstractArray) = (x = PaddedRFFTArray(x); $pf(x) * x)
+        $f(x::AbstractArray, region) = (x = PaddedRFFTArray(x); $pf(x, region) * x)
+        $f(x::PaddedRFFTArray) = $pf(x) * x
+        $f(x::PaddedRFFTArray, region) = $pf(x, region) * x
+        $pf(x::PaddedRFFTArray; kws...) = $pf(x, 1:ndims(x); kws...)
+    end
+end
+
 for f in (:brfft!, :irfft!)
     pf = Symbol("plan_", f)
     @eval begin
-        $f(x::AbstractArray, d::Integer) = (x = PaddedRFFTArray(x); $pf(x, d) * x)
-        $f(x::AbstractArray, d::Integer, region) = (x = PaddedRFFTArray(x); $pf(x, d, region) * x)
+        $f(x::AbstractArray, d::Integer) = $f(PaddedRFFTArray(x),d)
+        $f(x::AbstractArray, d::Integer, region) = $f(PaddedRFFTArray(x), d, region)
+        $f(x::PaddedRFFTArray, d::Integer) = $pf(x, d) * x
+        $f(x::PaddedRFFTArray, d::Integer, region) = $pf(x, d, region) * x
         $pf(x::PaddedRFFTArray, d::Integer; kws...) = $pf(x, d, 1:ndims(x);kws...)
-        $f(x::AbstractArray{<:Real}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
-        $f(x::AbstractArray{<:Complex{<:Union{Integer,Rational}}}, d::Integer, region=1:ndims(x)) = $f(complexfloat(x), d, region)
     end
 end
 
 for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
     # Note: use $FORWARD and $BACKWARD below because of issue #9775
     @eval begin
-        function plan_rfft(X::FFTStridedArray{$Tr,N}, region;
+        function plan_rfft(X::StridedArray{$Tr,N}, region;
                            flags::Integer=ESTIMATE,
                            timelimit::Real=NO_TIMELIMIT) where N
             osize = rfft_output_size(X, region)
@@ -889,7 +823,7 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             rFFTWPlan{$Tr,$FORWARD,false,N}(X, Y, region, flags, timelimit)
         end
 
-        function plan_brfft(X::FFTStridedArray{$Tc,N}, d::Integer, region;
+        function plan_brfft(X::StridedArray{$Tc,N}, d::Integer, region;
                             flags::Integer=ESTIMATE,
                             timelimit::Real=NO_TIMELIMIT) where N
             osize = brfft_output_size(X, d, region)
@@ -908,15 +842,23 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             end
         end
 
-        plan_rfft(X::FFTStridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...)
-        plan_brfft(X::FFTStridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...)
+        plan_rfft(X::StridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...)
+        plan_brfft(X::StridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...)
 
-        function plan_brfft!(X::PaddedRFFTArray{$Tr,N}, ::Integer, region;
+        function plan_rfft!(X::PaddedRFFTArray{$Tr,N}, region;
+                            flags::Integer=ESTIMATE,
+                            timelimit::Real=NO_TIMELIMIT) where N
+            rFFTWPlan{$Tr,$FORWARD,true,N}(real_view(X), complex_view(X), region, flags, timelimit)
+        end
+
+        function plan_brfft!(X::PaddedRFFTArray{$Tr,N}, d::Integer, region;
             flags::Integer=ESTIMATE,
             timelimit::Real=NO_TIMELIMIT) where N
+            @assert size(complex_view(X))[first(region)] == d>>1 + 1
             return rFFTWPlan{$Tc,$BACKWARD,true,N}(complex_view(X), real_view(X), region, flags, timelimit)
         end
 
+        plan_rfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_rfft!(X,1:ndims(X);kws...)
         plan_brfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_brfft!(X,1:ndims(X);kws...)
 
         function plan_irfft!(x::PaddedRFFTArray{$Tr,N}, d::Integer, region; kws...) where N
@@ -943,25 +885,25 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
                        normalization(Y, p.region))
         end
 
-        function mul!(y::FFTStridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::FFTStridedArray{$Tr})
+        function mul!(y::StridedArray{$Tc}, p::rFFTWPlan{$Tr,$FORWARD}, x::StridedArray{$Tr})
             assert_applicable(p, x, y)
             unsafe_execute!(p, x, y)
             return y
         end
-        function mul!(y::FFTStridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::FFTStridedArray{$Tc})
+        function mul!(y::StridedArray{$Tr}, p::rFFTWPlan{$Tc,$BACKWARD}, x::StridedArray{$Tc})
             assert_applicable(p, x, y)
             unsafe_execute!(p, x, y) # note: may overwrite x as well as y!
             return y
         end
 
-        function *(p::rFFTWPlan{$Tr,$FORWARD,false}, x::FFTStridedArray{$Tr,N}) where N
+        function *(p::rFFTWPlan{$Tr,$FORWARD,false}, x::StridedArray{$Tr,N}) where N
             assert_applicable(p, x)
             y = Array{$Tc}(undef, p.osz)
             unsafe_execute!(p, x, y)
             return y
         end
 
-        function *(p::rFFTWPlan{$Tc,$BACKWARD,false}, x::FFTStridedArray{$Tc,N}) where N
+        function *(p::rFFTWPlan{$Tc,$BACKWARD,false}, x::StridedArray{$Tc,N}) where N
             if p.flags & PRESERVE_INPUT != 0
                 assert_applicable(p, x)
                 y = Array{$Tr}(undef, p.osz)
@@ -975,6 +917,14 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             return y
         end
 
+        *(p::rFFTWPlan{$Tr,$FORWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N = 
+            (mul!(complex_view(f), p, real_view(f)); complex_view(f))
+
+        function \(p::rFFTWPlan{$Tr,$FORWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N
+            isdefined(p, :pinv) || (p.pinv = plan_irfft!(f, p.region))
+            return p.pinv * f
+        end
+
         *(p::rFFTWPlan{$Tc,$BACKWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N = 
             (mul!(real_view(f), p, complex_view(f)); real_view(f))
     end
@@ -1001,14 +951,14 @@ for f in (:r2r, :r2r!)
     end
 end
 
-function plan_r2r(X::FFTStridedArray{T,N}, kinds, region;
+function plan_r2r(X::StridedArray{T,N}, kinds, region;
                   flags::Integer=ESTIMATE,
                   timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N}
     r2rFFTWPlan{T,Any,false,N}(X, fakesimilar(flags, X, T), region, kinds,
                                flags, timelimit)
 end
 
-function plan_r2r!(X::FFTStridedArray{T,N}, kinds, region;
+function plan_r2r!(X::StridedArray{T,N}, kinds, region;
                    flags::Integer=ESTIMATE,
                    timelimit::Real=NO_TIMELIMIT) where {T<:fftwNumber,N}
     r2rFFTWPlan{T,Any,true,N}(X, X, region, kinds, flags, timelimit)
@@ -1043,20 +993,20 @@ function plan_inv(p::r2rFFTWPlan{T,K,inplace,N}) where {T<:fftwNumber,K,inplace,
                              1:length(iK)))
 end
 
-function mul!(y::FFTStridedArray{T}, p::r2rFFTWPlan{T}, x::FFTStridedArray{T}) where T
+function mul!(y::StridedArray{T}, p::r2rFFTWPlan{T}, x::StridedArray{T}) where T
     assert_applicable(p, x, y)
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::r2rFFTWPlan{T,K,false}, x::FFTStridedArray{T,N}) where {T,K,N}
+function *(p::r2rFFTWPlan{T,K,false}, x::StridedArray{T,N}) where {T,K,N}
     assert_applicable(p, x)
     y = Array{T}(undef, p.osz)::Array{T,N}
     unsafe_execute!(p, x, y)
     return y
 end
 
-function *(p::r2rFFTWPlan{T,K,true}, x::FFTStridedArray{T}) where {T,K}
+function *(p::r2rFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
     assert_applicable(p, x)
     unsafe_execute!(p, x, x)
     return x
diff --git a/test/runtests.jl b/test/runtests.jl
index 2291003..7785bcb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -528,3 +528,51 @@ end
         @test occursin("dft-thr", string(p2))
     end
 end
+
+let a = rand(Float64,(8,4,4)), b = PaddedRFFTArray(a), c = copy(b)
+    @testset "PaddedRFFTArray creation" begin
+      @test a == FFTW.real_view(b)
+      @test c == b
+      @test c.r == b.r
+    end
+    
+    @testset "rfft! and irfft!" begin
+      @test rfft(a) ≈ rfft!(b) 
+      @test a ≈ irfft!(b, size(a, 1))
+      @test rfft(a, 1:2) ≈ rfft!(b, 1:2) 
+      @test a ≈ irfft!(b, size(a, 1), 1:2)
+      # @test rfft(a, (1,3)) ≈ rfft!(b, (1,3)) fails 
+      # @test a ≈ irfft!(b, size(a, 1), (1,3))
+      
+      p = plan_rfft!(c)
+      @test p*c ≈ rfft!(b)
+      @test p\c ≈ irfft!(b, size(a, 1))
+    
+      a = rand(Float64, (9,4,4))
+      b = PaddedRFFTArray(a)
+      @test a == FFTW.real_view(b)
+      @test rfft(a) ≈ rfft!(b) 
+      @test a ≈ irfft!(b, size(a, 1))
+      @test rfft(a, 1:2) ≈ rfft!(b, 1:2) 
+      @test a ≈ irfft!(b, size(a, 1), 1:2)
+      # @test rfft(a, (1,3)) ≈ rfft!(b, (1,3)) # fails?
+      # @test a ≈ irfft!(b, size(a, 1), (1,3))
+    end
+    
+    @testset "brfft!" begin
+      a = rand(Float64,(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
+end #let block

From a9b7863069c1acfc82dec373c57563c1fb2683a9 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Wed, 29 Sep 2021 11:00:45 -0400
Subject: [PATCH 3/4] Subset to support irfft!

---
 src/FFTW.jl      |  2 +-
 src/fft.jl       | 90 +++++++++---------------------------------------
 test/runtests.jl | 60 +++++++++++---------------------
 3 files changed, 38 insertions(+), 114 deletions(-)

diff --git a/src/FFTW.jl b/src/FFTW.jl
index 2274354..0493cbc 100644
--- a/src/FFTW.jl
+++ b/src/FFTW.jl
@@ -12,7 +12,7 @@ import AbstractFFTs: Plan, ScaledPlan,
                      rfft_output_size, brfft_output_size,
                      plan_inv, normalization
 
-export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!, brfft!, irfft!, plan_brfft!, plan_irfft!, rfft!, plan_rfft!, PaddedRFFTArray
+export dct, idct, dct!, idct!, plan_dct, plan_idct, plan_dct!, plan_idct!, brfft!, irfft!, plan_brfft!, plan_irfft!
 
 include("providers.jl")
 
diff --git a/src/fft.jl b/src/fft.jl
index ad2dde0..8a6ac13 100644
--- a/src/fft.jl
+++ b/src/fft.jl
@@ -1,6 +1,6 @@
 # This file was formerly a part of Julia. License is MIT: https://julialang.org/license
 
-import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, copy, similar, eltype, \, getindex
+import Base: show, *, convert, unsafe_convert, size, strides, ndims, pointer, copy, getindex
 import LinearAlgebra: mul!
 
 
@@ -110,57 +110,31 @@ const fftwSingle = Union{Float32,Complex{Float32}}
 const fftwTypeDouble = Union{Type{Float64},Type{Complex{Float64}}}
 const fftwTypeSingle = Union{Type{Float32},Type{Complex{Float32}}}
 
-struct PaddedRFFTArray{T, N, A, Ar, Ac} <: AbstractArray{T, N} where {A <: Union{AbstractArray{T, N}, AbstractArray{Complex{T}, N}}, Ar <: AbstractArray{T, N}, Ac <: AbstractArray{Complex{T}, N}}
-    parent::A
+# padded array type to support RFFT in-place operations
+struct PaddedRFFTArray{T, N, Ac, Ar} <: AbstractArray{T, N} where {Ac <: AbstractArray{Complex{T}, N}, Ar <: AbstractArray{T, N}}
+    c::Ac # complex view / underlying array
     r::Ar # real view
-    c::Ac # complex view
 
     # wrap existing complex array
-    function PaddedRFFTArray(data::StridedArray{Tc, N}) where {Tc <: fftwComplex, N}
-        fsize = 2 * (size(data, 1) - 1)
-        r = view(reinterpret(real(Tc), data), Base.OneTo(fsize), ntuple(i -> Colon(), Val(ndims(data) - 1))...)
-        new{real(Tc), N, typeof(data), typeof(r), typeof(data)}(data, r, data)
-    end
-
-    # wrap existing padded reals array
-    function PaddedRFFTArray(data::StridedArray{T, N}, nx) where {T <: fftwReal, N}
-        fsize = size(data, 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"))
-        c = reinterpret(Complex{T}, data)
-        r = view(data, Base.OneTo(nx), ntuple(i->Colon(), Val(ndims(data) - 1))...)
-        new{T, N, typeof(data), typeof(r), typeof(c)}(data, r, c)
-    end
-
-    # copy reals array to padded array
-    function PaddedRFFTArray(reals::StridedArray{T, N}) where {T <: fftwReal, N}
-        fsize = (size(reals, 1) ÷ 2 + 1) * 2
-        data = Array{T, N}(undef, (fsize, size(reals)[2:end]...))
-        c = reinterpret(Complex{T}, data)
-        r = view(data, Base.OneTo(size(reals, 1)), ntuple(i->Colon(), Val(ndims(reals) - 1))...)
-        @inbounds copyto!(r, reals)
-        new{T, N, typeof(data), typeof(r), typeof(c)}(data, r, c)
+    function PaddedRFFTArray(a::StridedArray{Complex{T}, N}) where {T, N}
+        fsize = 2 * (size(a, 1) - 1)
+        r = view(reinterpret(T, a), Base.OneTo(fsize), ntuple(i -> Colon(), Val(ndims(a) - 1))...)
+        new{T, N, typeof(a), typeof(r)}(a, r)
     end
 
     # copy existing padded array
-    function PaddedRFFTArray(a::PaddedRFFTArray{T, N, A, Ar, Ac}) where {T, N, A, Ar, Ac}
-        data = copy(parent(a))
-        c = reinterpret(Complex{T}, data)
-        r = view(data, Base.OneTo(size(real_view(a), 1)), ntuple(i->Colon(), Val(ndims(data) - 1))...)
-        new{T, N, A, Ar, Ac}(data, r, c)
+    function PaddedRFFTArray(a::PaddedRFFTArray{T, N, Ac, Ar}) where {T, N, Ac, Ar}
+        c = copy(complex_view(a))
+        r = view(reinterpret(T, c), Base.OneTo(size(real_view(a), 1)), ntuple(i -> Colon(), Val(ndims(a) - 1))...)
+        new{T, N, Ac, Ar}(c, r)
     end
 end
 
 @inline real_view(S::PaddedRFFTArray) = S.r
-
 @inline complex_view(S::PaddedRFFTArray) = S.c
-
 size(a::PaddedRFFTArray) = size(complex_view(a))
 getindex(a::PaddedRFFTArray, args...) = getindex(complex_view(a), args...)
 copy(a::PaddedRFFTArray) = PaddedRFFTArray(a)
-parent(a::PaddedRFFTArray) = a.parent
 
 # For ESTIMATE plans, FFTW allows one to pass NULL for the array pointer,
 # since it is not written to.  Hence, it is convenient to create an
@@ -569,13 +543,7 @@ unsafe_execute!(plan::r2rFFTWPlan{T},
 # Compute dims and howmany for FFTW guru planner
 function dims_howmany(X::StridedArray, Y::StridedArray,
                       sz::Array{Int,1}, region)
-    reg = if isa(region, Int)
-        region:region
-    elseif isa(region, Tuple)
-        Int[region...]
-    else
-        Int.(region)
-    end
+    reg = Int[region...]
     if length(unique(reg)) < length(reg)
         throw(ArgumentError("each dimension can be transformed at most once"))
     end
@@ -788,27 +756,16 @@ function *(p::cFFTWPlan{T,K,true}, x::StridedArray{T}) where {T,K}
     return x
 end
 
-# rfft/brfft and planned variants.  No in-place version for now.
-
-for f in (:rfft!,)
-    pf = Symbol("plan_", f)
-    @eval begin
-        $f(x::AbstractArray) = (x = PaddedRFFTArray(x); $pf(x) * x)
-        $f(x::AbstractArray, region) = (x = PaddedRFFTArray(x); $pf(x, region) * x)
-        $f(x::PaddedRFFTArray) = $pf(x) * x
-        $f(x::PaddedRFFTArray, region) = $pf(x, region) * x
-        $pf(x::PaddedRFFTArray; kws...) = $pf(x, 1:ndims(x); kws...)
-    end
-end
+# rfft/brfft and planned variants.  No in-place version of rfft for now.
 
 for f in (:brfft!, :irfft!)
     pf = Symbol("plan_", f)
     @eval begin
-        $f(x::AbstractArray, d::Integer) = $f(PaddedRFFTArray(x),d)
+        $f(x::AbstractArray, d::Integer) = $f(PaddedRFFTArray(x), d)
         $f(x::AbstractArray, d::Integer, region) = $f(PaddedRFFTArray(x), d, region)
         $f(x::PaddedRFFTArray, d::Integer) = $pf(x, d) * x
         $f(x::PaddedRFFTArray, d::Integer, region) = $pf(x, d, region) * x
-        $pf(x::PaddedRFFTArray, d::Integer; kws...) = $pf(x, d, 1:ndims(x);kws...)
+        $pf(x::PaddedRFFTArray, d::Integer; kws...) = $pf(x, d, 1:ndims(x); kws...)
     end
 end
 
@@ -845,12 +802,6 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
         plan_rfft(X::StridedArray{$Tr};kws...)=plan_rfft(X,1:ndims(X);kws...)
         plan_brfft(X::StridedArray{$Tr};kws...)=plan_brfft(X,1:ndims(X);kws...)
 
-        function plan_rfft!(X::PaddedRFFTArray{$Tr,N}, region;
-                            flags::Integer=ESTIMATE,
-                            timelimit::Real=NO_TIMELIMIT) where N
-            rFFTWPlan{$Tr,$FORWARD,true,N}(real_view(X), complex_view(X), region, flags, timelimit)
-        end
-
         function plan_brfft!(X::PaddedRFFTArray{$Tr,N}, d::Integer, region;
             flags::Integer=ESTIMATE,
             timelimit::Real=NO_TIMELIMIT) where N
@@ -858,7 +809,6 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             return rFFTWPlan{$Tc,$BACKWARD,true,N}(complex_view(X), real_view(X), region, flags, timelimit)
         end
 
-        plan_rfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_rfft!(X,1:ndims(X);kws...)
         plan_brfft!(X::PaddedRFFTArray{$Tr};kws...)=plan_brfft!(X,1:ndims(X);kws...)
 
         function plan_irfft!(x::PaddedRFFTArray{$Tr,N}, d::Integer, region; kws...) where N
@@ -917,14 +867,6 @@ for (Tr,Tc) in ((:Float32,:(Complex{Float32})),(:Float64,:(Complex{Float64})))
             return y
         end
 
-        *(p::rFFTWPlan{$Tr,$FORWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N = 
-            (mul!(complex_view(f), p, real_view(f)); complex_view(f))
-
-        function \(p::rFFTWPlan{$Tr,$FORWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N
-            isdefined(p, :pinv) || (p.pinv = plan_irfft!(f, p.region))
-            return p.pinv * f
-        end
-
         *(p::rFFTWPlan{$Tc,$BACKWARD,true,N}, f::PaddedRFFTArray{$Tr,N}) where N = 
             (mul!(real_view(f), p, complex_view(f)); real_view(f))
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7785bcb..8c5b95f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -529,50 +529,32 @@ end
     end
 end
 
-let a = rand(Float64,(8,4,4)), b = PaddedRFFTArray(a), c = copy(b)
+begin
     @testset "PaddedRFFTArray creation" begin
-      @test a == FFTW.real_view(b)
-      @test c == b
-      @test c.r == b.r
+        a = rand(Complex{Float64}, (8, 4, 4))
+        b = FFTW.PaddedRFFTArray(a)
+        c = copy(b)
+        @test a == FFTW.complex_view(b)
+        @test c == b
+        @test FFTW.real_view(c) == FFTW.real_view(b)
+        @test FFTW.complex_view(c) == FFTW.complex_view(b)
     end
     
-    @testset "rfft! and irfft!" begin
-      @test rfft(a) ≈ rfft!(b) 
-      @test a ≈ irfft!(b, size(a, 1))
-      @test rfft(a, 1:2) ≈ rfft!(b, 1:2) 
-      @test a ≈ irfft!(b, size(a, 1), 1:2)
-      # @test rfft(a, (1,3)) ≈ rfft!(b, (1,3)) fails 
-      # @test a ≈ irfft!(b, size(a, 1), (1,3))
-      
-      p = plan_rfft!(c)
-      @test p*c ≈ rfft!(b)
-      @test p\c ≈ irfft!(b, size(a, 1))
-    
-      a = rand(Float64, (9,4,4))
-      b = PaddedRFFTArray(a)
-      @test a == FFTW.real_view(b)
-      @test rfft(a) ≈ rfft!(b) 
-      @test a ≈ irfft!(b, size(a, 1))
-      @test rfft(a, 1:2) ≈ rfft!(b, 1:2) 
-      @test a ≈ irfft!(b, size(a, 1), 1:2)
-      # @test rfft(a, (1,3)) ≈ rfft!(b, (1,3)) # fails?
-      # @test a ≈ irfft!(b, size(a, 1), (1,3))
+    @testset "irfft!" begin
+        a = rand(Float64, (8, 4, 4))
+        c = rfft(a)
+        d = copy(c)
+        e = irfft!(d, size(a, 1))
+        @test a ≈ e
+        @test irfft(c, size(a, 1)) ≈ e
+        @test d === parent(parent(e))
+        @test d != c
+        @test irfft(c, size(a, 1), 1:2) ≈ irfft!(copy(c), size(a, 1), 1:2)
     end
     
     @testset "brfft!" begin
       a = rand(Float64,(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)
+      b = rfft(a)
+      @test (brfft!(b, size(a, 1)) ./ 16) ≈ a
     end
-end #let block
+end

From 11926ce76ad2b41c0948e03f1c1f7804e2f4f7a3 Mon Sep 17 00:00:00 2001
From: Nicholas Bauer <nicholasbauer@outlook.com>
Date: Wed, 29 Sep 2021 11:06:49 -0400
Subject: [PATCH 4/4] Restore accidental lines

---
 src/fft.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/fft.jl b/src/fft.jl
index 8a6ac13..80f03e2 100644
--- a/src/fft.jl
+++ b/src/fft.jl
@@ -610,7 +610,7 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
                                                      Y::StridedArray{$Tc,N},
                                                      region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
-        region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
+        region = circshift(Int[region...],-1) # FFTW halves last dim
         unsafe_set_timelimit($Tr, timelimit)
         dims, howmany = dims_howmany(X, Y, [size(X)...], region)
         plan = ccall(($(string(fftw,"_plan_guru64_dft_r2c")),$lib[]),
@@ -630,7 +630,7 @@ for (Tr,Tc,fftw,lib) in ((:Float64,:(Complex{Float64}),"fftw",:libfftw3),
                                                       Y::StridedArray{$Tr,N},
                                                       region, flags::Integer, timelimit::Real) where {inplace,N}
         R = isa(region, Tuple) ? region : copy(region)
-        region = isa(region, AbstractArray) ? circshift(region, -1) : region # FFTW halves last dim
+        region = circshift(Int[region...],-1) # FFTW halves last dim
         unsafe_set_timelimit($Tr, timelimit)
         dims, howmany = dims_howmany(X, Y, [size(Y)...], region)
         plan = ccall(($(string(fftw,"_plan_guru64_dft_c2r")),$lib[]),