What FFT methods should I use?

[brianreinhold]

I have a spectral transform model that every time step has to

1. transform spectral components to grid point data

2. Do nonlinear computations in grid point space

3. transform back to spectral space.

The spectral components are all real, using eigenfunctions

streamfunction(k,l) = a(k,l)cosknxsinly + b(k,l)sinknxsinly + cosly

Derivatives are taken for the Jacobian so I will also have to transform
d/dy = a(k,l)cosknxcosly + b(k,l)sinknxcosly + sinly

The domain is periodic in x and bounded in y with range y=0 to y=pi.

The model is written but used the old F77 FFT991 FFT routines. I wanted to update to F90+ at least.

I started with FFTW3 but ran into problems so I downloaded MKL and started with the DFTI methods. However, in spite of the fact all my data is real, I had to allocate a complex array to go from spectral to grid point space. I have many arrays in this project and they are all real and are allocated and constructed using pointers.

Then I saw the TT methods in the partial differential equations helpers. That is exactly what I am integrating so I thought this was good. But as far as I could tell I had to revert to using the DFTI routines for at least one set of these transforms (likely spectral to grid data) which meant allocating complex arrays again.

So I gave up and went back to FFTW3 but could only get a precompiled version for Windows for 3.3.5 which is aways behind 3.3.10. Building them myself seemed painfully complicated so I went for the FFTW3 wrappers in the MKL. As far as I can tell, this appears to be my best solution if I don't want to reconstruct my arrays to be complex (with wasted space).

Seems funny I would have to do this. Maybe I misunderstood the TT options and they can apply to my case using all real arrays ...?

[Fengrui]

Thank you for posting in the forum!

However, in spite of the fact all my data is real, I had to allocate a complex array to go from spectral to grid point space.

>> Actually you don't have to allocate a complex array. In ${MKLROOT}/include/mkl_dfti.h, the compute API is declared as "DFTI_EXTERN MKL_LONG DftiComputeForward(DFTI_DESCRIPTOR_HANDLE, void*, ...);", the arrays got passed as void pointers. So, for example, if using double* a to store the complex data in MKL_Complex16 *b, a[0] takes b[0].real, a[1] takes b[0].imag, a[2] takes b[1].real, ...