[SYMSCAE]

Summary

Speeds up PDFCalculator, DebyePDFCalculator and bond generation for periodic and crystal structures. Most of the gain comes from SIMD calculation kernels, batched bond generation, and dropping uBLAS in favor of fixed-size vector and matrix types.

Changes per commit

1. Fixed-size linear algebra and hashing

• R3::Vector and R3::Matrix are now std::array types.
• Adds diffpy/hash.hpp which uses the same mixing scheme as Boost.
• Removes the dependencies on Boost.Functional/Hash, Boost.StringAlgo and Boost.uBLAS.
• R3::inverse, R3::floor, R3::mxvecproduct and the Lattice coordinate and matrix methods now return by value instead of a reference to a shared static. This makes them thread-safe for future multithreaded use.

2. SIMD kernels for the Gaussian and Debye loops

• Adds kernels for AVX-512, AVX2 and SSE2 on x86-64, and NEON on arm64 based on xsimd and selects the best kernel the CPU supports on the fly with manual override option available.
• All kernels are compiled with IEEE math (-fno-fast-math -ffp-contract=off) so they give consistent results. The SIMD files are built with -fno-lto so their ISA-specific inline functions can't leak into the rest of the library.
• New public method PeakProfile::addRDFPeak. New protected hook GaussianProfile::kernelShape lets derived profiles use the kernels.

3. Batched periodic bond generation

• PeriodicStructureBondGenerator stores the unit cell positions and lattice vectors as separate cartesian arrays and finds candidate bonds in chunks.
• CrystalStructureAdapter reuses symmetry-expanded positions and expands again only when atoms, the lattice or the symmetry operations change.

4. Build: LTO and -fno-semantic-interposition

• Release and RelWithDebInfo builds of the shared library use LTO by default.
• The library is compiled with -fno-semantic-interposition when the compiler supports it.

5. Faster FFT

• fftgtof and fftftog compute the sine transform with a half-length complex FFT instead of a 4N-point odd extension. GSL wavetables are cached per thread for each power-of-two length.

Benchmarks

[SYMSCAE]

@sbillinge PR ready for review. CI failed due to codecov not enabled.

[sbillinge]

this is great. Here is a Fable review:

I read all the source changes in PR #76 and found no correctness bugs, but I could not build it or run the tests.

Things to raise with the author

1. Scope. The PR has five themes in about 4,200 lines: the R3/hash rewrite, SIMD kernels, batched bond generation, LTO flags, and the FFT rewrite. The five commits are already cleanly separable, so I'd ask for a split into separate PRs peeling off the breaking R3 change.

2. Breaking API and ABI for downstream. The news file documents these, but they need a decision from you:

   • R3::Vector and R3::Matrix no longer accept uBLAS expressions, and R3::row(A, 0) = v no longer compiles.
   • Several R3 and Lattice functions now return by value.
   • Protected members of the periodic and crystal bond generators are removed, and PeriodicStructureBondGenerator is non-copyable.
   • PeakProfile gains a virtual method and BaseBondGenerator gains members.

3. SIMD silently off without xsimd. LIBDIFFPY_SIMD=AUTO falls back to the scalar kernel with only a STATUS message. xsimd is added to requirements/build.txt with no minimum version, and the conda-forge feedstock would need it too, otherwise users get the slow build without knowing. [simon says, I think this is ok, but let's remember to make sure it is in the build requirements of the meta.yaml]

4. Instruction-set leakage risk. The SIMD files rely on -O2 inlining plus a renamed xsimd namespace to stop AVX-512/AVX2 copies of inline standard-library functions being shared with the rest of the library. That is a mitigation, not a guarantee. The #error guard for unoptimized builds only covers GCC/Clang, so an MSVC Debug build with /arch:AVX2 is unprotected. A CI check with nm for weak symbols in the kernel objects would close this.

5. SIMD falls back to scalar for far, narrow peaks. In CalculationKernelSIMD.ipp, gaussianWindow uses the scalar kernel when dist > ~450 * fwhm. For fwhm of 0.1 Å that is r above roughly 45 Å, so large-rmax calculations lose much of the speedup. Ask whether the benchmarks cover that case; the PR's benchmark is an image only. [@sbillinge asks why is this? We often want to compute beyond 45 A-1, and the speed of this calculator will only encourage more of that....is there a reason for that design choice?]

Minor

• CalculationKernelSSE2.cpp header comment says "NEON calculation kernel".
• R3linalg.ipp is a new file but carries "(c) 2009 … File coded by: Pavol Juhas", and the new files mix LICENSE.rst and LICENSE_DANSE.rst headers.
• fftgtof keeps per-thread FFT buffers for every power-of-two length for the life of the thread. Probably fine, but worth a comment.
• -fno-semantic-interposition and LTO are on by default; this is in the news, just be aware for packaging.

What I checked and found sound

• FFT rewrite (PDFUtils.cpp): I worked through the half-length transform algebra and it reproduces the old odd-extension result, including the small-N cases. Negative rmin now throws where it previously indexed out of bounds.
• Bond batching (PeriodicStructureAdapter.cpp): bond order matches the old nested loops for both periodic and crystal generators. The row-narrowing geometry is conservative, the exact distance check is still applied to every candidate, and buffer sizes are sufficient.
• Crystal symmetry caching: invalidation covers in-place atom and lattice edits as well as symmetry-operation and precision changes.
• Legacy serialization of R3::Vector and R3::Matrix: the layout matches the old uBLAS one, and there is a byte-level test for it.
• Evaluator fix after exceptions (PQEvaluator.cpp): correct.
• Kernels: the SIMD Debye early-stop and tail handling match the scalar semantics, and the cropped-Gaussian peak height is consistent.

[sbillinge]

Is this ready to merge? Can anyone see any issues? @cadenmyers13 please can we build an env with this version of libdiffpy and srreal and run all the examples?

[cadenmyers13]

@sbillinge @SYMSCAE Wow this is amazing!!! Ran this on the diffpy.cmi examples and everything looked good. The examples used to take ~20-25 min to run. This only took ~6 min

[SYMSCAE]

@sbillinge Investigating on msvc builds which the linker oddly expose vectorized stl functions out of scope. Should be ready within an hour or so.