fix(FFTXlib): add FFTW_UNALIGNED to cached FFTW plans to match new-array execute

Maurice Hanischrequested to merge
MauriceDHanisch/q-e:fix-fftw-unaligned-plan-execute into develop

Summary

The FFTW3 backend (FFTXlib/src/fft_scalar.FFTW3.f90) builds cached plans on a freshly allocated scratch array (f_test, ALLOCATE(..., mold=...)) and later executes them via fftw_execute_dft on the caller's array f, which can have a different SIMD alignment. FFTW's new-array execute requires the execute array to share the planning array's alignment; otherwise the behavior is undefined.

With an FFTW built for AVX-512 (e.g. 3.3.10 DYNAMIC_ARCH), the cached plan selects an alignment-assuming codelet and silently computes a wrong transform when f is not 64-byte aligned. In pw.x this corrupts the wavefunction FFT (cfft3ds), then the overlap matrix S, making the Cholesky in cdiaghg/rdiaghg fail with S matrix not positive definite on systems that converge with less aggressive FFTW builds (e.g. a plain SSE2/AVX FFTW 3.3.8).

cfft3ds additionally reuses one cached plan across many 16-byte-shifted offsets f(i:)/f(ii:), so no allocation-alignment scheme can satisfy every execute site; the plans must be alignment-agnostic. This MR adds FFTW_UNALIGNED to all cached plan creations in cft_1z, cft_2xy, cfft3d, and cfft3ds (16 sites).

Why FFTW_UNALIGNED rather than forcing alignment

A single cached plan is executed both on different caller arrays (different alignments) and, within cfft3ds, on many starts 16 bytes apart (f(i:), f(ii:)). No allocation alignment makes all of those simultaneously 64-byte aligned, so the plans are inherently used in an alignment-agnostic way. FFTW_UNALIGNED describes exactly how the code already runs.

Reproducer

A single Mg atom in an 8 Å vacuum box, gamma point, 40/160 Ry, with pw.x linked against an AVX-512 FFTW (e.g. spack fftw@3.3.10 with DYNAMIC_ARCH):

&CONTROL
  calculation='scf'
  prefix='m'
  pseudo_dir='.'
/
&SYSTEM
  ibrav=0, nat=1, ntyp=1
  ecutwfc=40, ecutrho=160
  occupations='smearing', smearing='mv', degauss=0.01
/
&ELECTRONS
  conv_thr=1.0d-8
/
ATOMIC_SPECIES
Mg 24.305 Mg_ONCV_PBE-1.0.upf
CELL_PARAMETERS angstrom
 8.0 0.0 0.0
 0.0 8.0 0.0
 0.0 0.0 8.0
ATOMIC_POSITIONS angstrom
Mg 0.0 0.0 0.0
K_POINTS gamma

Before this patch the run aborts at the first diagonalization:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine cdiaghg (41):
     S matrix not positive definite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The same input and binary converge when linked against a plain SSE2/AVX FFTW (3.3.8), which is what made this look machine-dependent.

Verification

With the patch, the failing SCF converges and the total energy matches a known-good FFTW build to all printed digits. On a separate MgO test the energy is bit-identical to an MKL/DFTI build of the same input. The only cost is a small FFT slowdown from disabling the aligned codelet (single-digit percent on realistic systems; larger on FFT-bound cases), which is the price of correct results on alignment-aggressive FFTW builds.

Merge request reports