# \`plan\_rfft\` together used with \`mul!\` and \`@view\` complains about wrong-strides-output
**URL:** https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215
**Category:** Signal and Image Processing
**Tags:** fftw, views
**Created:** [March 30, 2021, 11:05am UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215 "2021-03-30T11:05:10Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 11:05am UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/1 "2021-03-30T11:05:10Z")
Hey,
initially I wanted to have `rfft!` but since it doesn’t exist, I tried the workaround mentioned [in this discourse discussion](https://discourse.julialang.org/t/in-place-rfft-with-padded-arrays/3340).
However, `mul!` complains when used with a `plan_rfft` generated matrix:
```julia
using LinearAlgebra, FFTW
function main()
# array we want to process
x = randn((60,40))
y = rfft(randn((80,60)))
# just create a view storage
y2 = view(y, 1:31, 1:40)
p = plan_rfft(x)
@time y2 .= p * x; # this still allocates memory because of p * x
@time y2 .= p * x;
@time mul!(y2, p, x); # this should be memory allocation free
@time mul!(y2, p, x);
end
```
Output then is:
```julia
julia> main()
0.000020 seconds (2 allocations: 19.453 KiB)
0.000012 seconds (2 allocations: 19.453 KiB)
ERROR: ArgumentError: FFTW plan applied to wrong-strides output
Stacktrace:
[1] assert_applicable(p::FFTW.rFFTWPlan{Float64, -1, false, 2, UnitRange{Int64}}, X::Matrix{Float64}, Y::SubArray{ComplexF64, 2, Matrix{ComplexF64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false})
@ FFTW ~/.julia/packages/FFTW/G3lSO/src/fft.jl:435
[2] mul!
@ ~/.julia/packages/FFTW/G3lSO/src/fft.jl:781 [inlined]
[3] macro expansion
@ ./timing.jl:210 [inlined]
[4] main()
@ Main ~/julia/FourierTools.jl/src/rfft_issue.jl:17
[5] top-level scope
@ REPL[4]:1
```
Apparently FFTW.jl seems not to correctly accept the view as output.
Am I doing something wrong or what is the issue here?
Thanks,
Felix
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 3:58pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/2 "2021-03-30T15:58:45Z")
The output is non-contiguous and apparently FFTW cannot handle that.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 4:09pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/3 "2021-03-30T16:09:48Z")
Strange, that’s true. So `views` are somehow supported within `mul!` but only in very special cases:
Might that be easy fixed within FFTW.jl because usual Julia code works with such views?
The following works, where I cropped instead the second dimension instead of the first one (`31 == size(y, 1)`):
```julia
function main()
# array we want to process
x = randn((60,40))
y = rfft(randn((60,50)))
# just create a view storage which crops second dimension:
y2 = view(y, 1:31, 1:40)
p = plan_rfft(x, flags=FFTW.UNALIGNED)
@time y2 .= p * x;
@time mul!(y2, p, x);
end
```
There have been some [issues fixed](https://github.com/JuliaMath/FFTW.jl/issues/67) in the past regarding views and rfft but during assignment.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 4:17pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/4 "2021-03-30T16:17:55Z")
The issue is that you are taking a corner out of a 2D Array. If you take a contiguous chunk out of a 1D Vector it should be fine
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 4:20pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/5 "2021-03-30T16:20:52Z")
Memory is linear, see this graphic
So if you just take a corner from the array (with incomplete columns), there are gaps that FFTW would have to jump over.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 4:21pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/6 "2021-03-30T16:21:28Z")
Yes, but I actually want that.
I mean, is it just `mul!` or FFTW.jl? because a normal statement assigning an array to a new, even scattered view array just works. So it must be due to the C wrapper, right?
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 4:22pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/7 "2021-03-30T16:22:46Z")
Yes, `mul!` in this case is just calling another method from the C-library FFTW which is not set up to handle non-contiguous output.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 4:24pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/8 "2021-03-30T16:24:44Z")
Keep in mind you can just pre-allocate two arrays: one for the unpadded FFTW output and another one with padding.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 4:27pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/9 "2021-03-30T16:27:17Z")
The issue is, I’m calling several 1D rffts for different dimensions. Apparently in place rfft does not exist for single arbitrary dimensions.
That’s why I wanted to allocate one large array covering at least the size in each dimension needed to accommodate each single result.
For each dimension I simply would have needed a view halving the rfft dimension.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:29pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/10 "2021-03-30T18:29:38Z")
I’m not sure I understand correctly—are performing the FFT only along one axis?
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 6:34pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/11 "2021-03-30T18:34:43Z")
Yeah. It would be nice to prevent allocation of `arr_ft` every time, because x could be simply in-place overwritten. For the inverse it works, because the resulting size is always the same:
See below some part:
```julia
for d in dims
p = plan_rfft(arr, d)
arr_ft = p * arr
arr_ft .*= ϕ # do something with arr_ft
mul!(arr, inv(p), arr_ft)
end
```
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:45pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/12 "2021-03-30T18:45:03Z")
Ok, so firstly that is not what your code currently does. Your plan `p` performs a 2D FFT — because it was set up with a 2D input array.
If you want to perform a 1D FFT, you need to use a Vector to create the plan `p`. Then you could write a for-loop to go through your data column by column (you can use views for this) and apply `p` to each.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:46pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/13 "2021-03-30T18:46:49Z")
also, I’m pretty sure you dont need to use inv(p) — since that probably allocates a new plan for the inverse transform. There should be an `ldiv!()` method which works like `mul!()` but performs and inverse transform.
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 6:47pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/14 "2021-03-30T18:47:23Z")
I’m not sure if you are correct. I believe the `dims` argument for `fft, rfft` is intended for exactly that case:
`rfft` performs a `length(dims)`-dimensional Fourier transform over the `dims` dimensions of the N-dimensional array.
```julia
rfft(A [, dims])
Multidimensional FFT of a real array A, exploiting the fact that the transform has conjugate symmetry in order to save
roughly half the computational time and storage costs compared with fft. If A has size (n_1, ..., n_d), the result has
size (div(n_1,2)+1, ..., n_d).
The optional dims argument specifies an iterable subset of one or more dimensions of A to transform, similar to fft.
Instead of (roughly) halving the first dimension of A in the result, the dims[1] dimension is (roughly) halved in the
same way.
```
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:49pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/15 "2021-03-30T18:49:26Z")
But you are not using `dims` in the code you posted
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 6:49pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/16 "2021-03-30T18:49:33Z")
The issue with` ldiv!(A, B)` is that, both `A` and `B` have different types than `inv(P) * B`, so can’t work?
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 6:50pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/17 "2021-03-30T18:50:10Z")
`rfft(arr, d)` with `d` being an integer does that.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:50pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/18 "2021-03-30T18:50:22Z")
```julia
ldiv!(output, p, input)
```
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [March 30, 2021, 6:52pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/19 "2021-03-30T18:52:24Z")
where is the integer in your code?
---
### Author:  [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)
#### Post date: [March 30, 2021, 6:54pm UTC](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215/20 "2021-03-30T18:54:16Z")
Manual too long for `ldiv!` couldn’t see that 😂
`for d in ntuple(identity, N)` is what I’m using. I’m looping over the dimensions of the array.
`lmul!` doesn’t bring any speedup, I believe internally `inv` and `ldiv` are probably doing exactly the same thing
[Next page](https://discourse.julialang.org/t/plan-rfft-together-used-with-mul-and-view-complains-about-wrong-strides-output/58215.md?page=2)