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
/ input
of 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 usedfftw_execute_dft_r2c
, and real-output (c2r) DFT plans should usedfftw_execute_dft_c2r
. The various r2r plans should usedfftw_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 theFFTW_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)
-
reporter -
reporter Pass the arrays to use in transform_x/inverse_x
See issue #214 for further discussion
→ <<cset 25b16aca4789>>
-
reporter See PR #674 for a temporary solution.
- Log in to comment
Having a type which wraps all this up may make it a little easier to make use of
fftw_malloc
as this involves a few steps. The main downside I can see is that it may end up with us being more likely to have to copy data unless we’re careful and it’s likely to require some wider code changes.An alternative solution is to remove public visibility from
transform_x
as it isn’t actually used in other modules. We can then simply ensurexxf
is aligned (and private).