Caution regarding gs2_transforms

In attempting to write a few simple tests for gs2_transforms I hit a couple of what seemed like very confusing errors. I’ve now understood the origin of these and will record this here for others and to allow discussion about possible changes.

Ultimately the issues boil down to the fact that the transform_x5d and inverse_x5d methods take xxf as an array to store the output / inputof the transforms. However these are not passed to the fftw_execute call, which instead just takes the fft plan handle.

This is because with fftw3 one may choose to provide the actual array used to store these at the time of the plan creation and then these are automatically used when the plan is executed. It is possible to used advanced methods (e.g. fftw_execute_dft )which do take the input and output arrays and this is the approach taken in the y transforms*. This has no impact for the commonly used transform2 and inverse2 methods as xxf is only used to temporarily store intermediary results. However, as we expose transform_x for use outside of gs2_transforms this could cause problems for people in the future unaware of this restriction.

Indeed the FFTW documentation notes (https://www.fftw.org/fftw3_doc/FFTW-Execution-in-Fortran.html)

There are various workarounds to this, but the safest and simplest thing is to not use dfftw_execute in Fortran.

Furthermore, even using the methods which take the arrays is not without possible issue:

There are a few things to be careful of, however:

You must use the correct type of execute function, matching the way the plan was created. Complex DFT plans should use dfftw_execute_dft, Real-input (r2c) DFT plans should use use dfftw_execute_dft_r2c, and real-output (c2r) DFT plans should use dfftw_execute_dft_c2r. The various r2r plans should use dfftw_execute_r2r.

You should normally pass the same input/output arrays that were used when creating the plan. This is always safe.

If you pass different input/output arrays compared to those used when creating the plan, you must abide by all the restrictions of the new-array execute functions (see New-array Execute Functions). The most difficult of these, in Fortran, is the requirement that the new arrays have the same alignment as the original arrays, because there seems to be no way in legacy Fortran to obtain guaranteed-aligned arrays (analogous to fftw_malloc in C). You can, of course, use the FFTW_UNALIGNED flag when creating the plan, in which case the plan does not depend on the alignment, but this may sacrifice substantial performance on architectures (like x86) with SIMD instructions (see SIMD alignment and fftw_malloc).

We currently use FFTW_UNALIGNED when creating the plans so our use of these approaches in the y transforms should be safe, but potentially less performant than it could be. See https://www.fftw.org/fftw3_doc/New_002darray-Execute-Functions.html

‌

I believe a bit of restructuring to all of this might help. In particular wrapping up the plan, arrays and transform methods into a small derived type may provide a nice way to ensure consistency and access any performance gain from alignment.

Comments (3)