perf: lightning-fast vectorized PDF calculators

.github/workflows/matrix-and-codecov.yml

@@ -76,14 +76,14 @@ jobs:
           ctest --test-dir build\release --output-on-failure --no-tests=error
 
       - name: Generate coverage report
-        if: runner.os == 'Linux'
+        if: matrix.coverage == 'ON'
         run: >-
           gcovr --gcov-executable "$CONDA_PREFIX/bin/gcov"
           --root . --filter 'src/diffpy/'
           --xml coverage.xml
 
       - name: Save coverage report
-        if: runner.os == 'Linux'
+        if: matrix.coverage == 'ON'
         uses: actions/upload-artifact@v4
         with:
           name: coverage-report
@@ -92,7 +92,7 @@ jobs:
 
       - name: Upload coverage to Codecov
         # Fork push builds have no upstream Codecov token; keep their artifact.
-        if: runner.os == 'Linux' && github.repository == 'diffpy/libdiffpy'
+        if: matrix.coverage == 'ON' && github.repository == 'diffpy/libdiffpy'
         uses: codecov/codecov-action@v5
         with:
           files: coverage.xml

CMakeLists.txt

@@ -35,6 +35,8 @@ endif()
 option(LIBDIFFPY_BUILD_OBJECT_LIB "Build an OBJECT library for embedding" OFF)
 option(LIBDIFFPY_INSTALL "Install libdiffpy and its CMake package" ON)
 option(LIBDIFFPY_ENABLE_OBJCRYST "Build with the ObjCryst structure adapter" ON)
+option(LIBDIFFPY_LTO
+  "Link the optimized shared library with link-time optimization" ON)
 option(LIBDIFFPY_BUILD_TESTS "Build the Catch2 unit tests" OFF)
 option(LIBDIFFPY_TEST_INSTALLED
   "Build only the unit tests and link them with the installed libdiffpy" OFF)
@@ -294,6 +296,27 @@ endif()
 if(MSVC)
   # Required for the Boost.Serialization code, including embedded builds.
   target_compile_options(diffpy PRIVATE /bigobj)
+else()
+  # Let GCC inline calls between the exported functions of the library.
+  include(CheckCXXCompilerFlag)
+  check_cxx_compiler_flag(-fno-semantic-interposition
+    LIBDIFFPY_HAS_NO_SEMANTIC_INTERPOSITION)
+  if(LIBDIFFPY_HAS_NO_SEMANTIC_INTERPOSITION)
+    target_compile_options(diffpy PRIVATE -fno-semantic-interposition)
+  endif()
+  # LTO inlines the short per-pair calls between source files.  Static and
+  # OBJECT libraries skip it, their objects would need the same compiler.
+  if(LIBDIFFPY_LTO AND _libdiffpy_library_type STREQUAL "SHARED_LIBRARY")
+    include(CheckIPOSupported)
+    check_ipo_supported(RESULT _ipo_supported OUTPUT _ipo_error LANGUAGES CXX)
+    if(_ipo_supported)
+      set_target_properties(diffpy PROPERTIES
+        INTERPROCEDURAL_OPTIMIZATION_RELEASE ON
+        INTERPROCEDURAL_OPTIMIZATION_RELWITHDEBINFO ON)
+    else()
+      message(STATUS "Link-time optimization is not supported, skipping it.")
+    endif()
+  endif()
 endif()
 target_compile_features(diffpy PUBLIC cxx_std_23)
 target_include_directories(diffpy PUBLIC
@@ -319,6 +342,7 @@ target_link_libraries(diffpy
 if(LIBDIFFPY_HAS_OBJCRYST)
   target_link_libraries(diffpy PUBLIC ObjCryst::ObjCryst)
 endif()
+include(CalculationKernels)
 
 # Installation ----------------------------------------------------------------
 

cmake/CalculationKernels.cmake

@@ -0,0 +1,107 @@
+# Calculation kernels.  The scalar kernel is always built.  With LIBDIFFPY_SIMD
+# and xsimd, all SIMD kernels of the target architecture are built too and
+# the library selects the best one the CPU supports at runtime, see
+# src/diffpy/srreal/kernels/CalculationKernelDispatch.cpp.
+set(LIBDIFFPY_SIMD "AUTO" CACHE STRING
+  "SIMD calculation kernels: AUTO (all kernels of the target architecture \
+when xsimd is found) or OFF (scalar kernel only)")
+set_property(CACHE LIBDIFFPY_SIMD PROPERTY STRINGS AUTO OFF)
+
+# Target architecture.  A universal macOS build matches no architecture
+# below and gets only the scalar kernel.
+set(_kernel_arch "${CMAKE_SYSTEM_PROCESSOR}")
+if(MSVC AND CMAKE_CXX_COMPILER_ARCHITECTURE_ID)
+  set(_kernel_arch "${CMAKE_CXX_COMPILER_ARCHITECTURE_ID}")
+endif()
+if(APPLE AND CMAKE_OSX_ARCHITECTURES)
+  set(_kernel_arch "${CMAKE_OSX_ARCHITECTURES}")
+endif()
+string(TOLOWER "${_kernel_arch}" _kernel_arch)
+if(_kernel_arch MATCHES "^(x86_64|amd64|x64)$")
+  set(_kernel_arch "x86_64")
+elseif(_kernel_arch MATCHES "^(aarch64|arm64)$")
+  set(_kernel_arch "arm64")
+else()
+  set(_kernel_arch "")
+endif()
+
+# The SIMD kernels need xsimd.
+set(_kernel_simd FALSE)
+if(LIBDIFFPY_SIMD AND _kernel_arch)
+  if(NOT TARGET xsimd)
+    find_package(xsimd CONFIG QUIET)
+  endif()
+  if(TARGET xsimd)
+    set(_kernel_simd TRUE)
+  else()
+    message(STATUS "xsimd not found, building only the scalar calculation "
+      "kernel.  Install xsimd for the SIMD kernels.")
+  endif()
+endif()
+
+# The main source glob is not recursive, add the kernels explicitly.
+set(_kernel_dir "${PROJECT_SOURCE_DIR}/src/diffpy/srreal/kernels")
+target_sources(diffpy PRIVATE
+  "${_kernel_dir}/CalculationKernelScalar.cpp"
+  "${_kernel_dir}/CalculationKernelDispatch.cpp")
+set(_kernel_names scalar)
+
+# All kernels need IEEE arithmetic to give consistent results.
+if(MSVC)
+  set(_kernel_math_options /fp:precise)
+else()
+  set(_kernel_math_options -fno-fast-math -ffp-contract=off)
+endif()
+set_source_files_properties("${_kernel_dir}/CalculationKernelScalar.cpp"
+  PROPERTIES COMPILE_OPTIONS "${_kernel_math_options}")
+
+# Inline functions compiled for the SIMD instruction sets must not be shared
+# with the rest of the library, see CalculationKernelSIMD.ipp.  -O2 inlines
+# them also in unoptimized builds and -fno-lto keeps LTO from sharing them.
+set(_simd_isolation_options)
+if(_kernel_simd AND NOT MSVC)
+  include(CheckCXXCompilerFlag)
+  check_cxx_compiler_flag(-fno-lto LIBDIFFPY_HAS_NO_LTO)
+  if(LIBDIFFPY_HAS_NO_LTO)
+    list(APPEND _simd_isolation_options -fno-lto)
+  endif()
+  list(APPEND _simd_isolation_options
+    "$<$<NOT:$<CONFIG:Release,RelWithDebInfo>>:-O2>")
+endif()
+
+# Add the SIMD kernel in CalculationKernel<name>.cpp, compiled with the
+# instruction set flags in ARGN.  CalculationKernelDispatch.cpp uses the
+# kernel when LIBDIFFPY_KERNEL_<name> is defined.
+function(_libdiffpy_add_simd_kernel name)
+  set(_source "${_kernel_dir}/CalculationKernel${name}.cpp")
+  set(_options ${ARGN} ${_kernel_math_options} ${_simd_isolation_options})
+  set_source_files_properties("${_source}"
+    PROPERTIES COMPILE_OPTIONS "${_options}")
+  target_sources(diffpy PRIVATE "${_source}")
+  target_compile_definitions(diffpy PRIVATE LIBDIFFPY_KERNEL_${name}=1)
+  string(TOLOWER "${name}" _lower)
+  set(_kernel_names ${_kernel_names} ${_lower} PARENT_SCOPE)
+endfunction()
+
+if(_kernel_simd AND _kernel_arch STREQUAL "x86_64")
+  _libdiffpy_add_simd_kernel(SSE2)
+  if(MSVC)
+    _libdiffpy_add_simd_kernel(AVX2 /arch:AVX2)
+    _libdiffpy_add_simd_kernel(AVX512 /arch:AVX512)
+  else()
+    _libdiffpy_add_simd_kernel(AVX2 -mavx2 -mfma)
+    _libdiffpy_add_simd_kernel(AVX512
+      -mavx512f -mavx512cd -mavx512dq -mfma)
+  endif()
+elseif(_kernel_simd AND _kernel_arch STREQUAL "arm64")
+  _libdiffpy_add_simd_kernel(NEON)
+endif()
+list(JOIN _kernel_names ", " _kernel_names)
+message(STATUS "Calculation kernels: ${_kernel_names}")
+
+if(_kernel_simd)
+  target_compile_definitions(diffpy PRIVATE LIBDIFFPY_SIMD_KERNELS=1)
+  # Include path only, keeps xsimd out of the static library's link interface.
+  target_include_directories(diffpy SYSTEM PRIVATE
+    $<TARGET_PROPERTY:xsimd,INTERFACE_INCLUDE_DIRECTORIES>)
+endif()

docs/source/development.rst

@@ -45,7 +45,24 @@ To use the same options for every build, set them in a ``CMakeUserPresets.json``
 file in the repository root. To use another generator, configure a separate
 directory with ``cmake -S . -B <build-dir> -G "<generator>"``. To disable the
 ObjCryst adapter, pass ``-DLIBDIFFPY_ENABLE_OBJCRYST=OFF``. For gprof
-instrumentation, use ``-DLIBDIFFPY_PROFILE=ON``.
+instrumentation, use ``-DLIBDIFFPY_PROFILE=ON``. To build the optimized
+shared library without link-time optimization, pass ``-DLIBDIFFPY_LTO=OFF``.
+
+Calculation kernels
+-------------------
+
+``PDFCalculator`` and ``DebyePDFCalculator`` run their inner loops in a
+calculation kernel. The scalar kernel is always built. The CMake option
+``LIBDIFFPY_SIMD`` controls the SIMD kernels:
+
+* ``AUTO`` - the default, build all SIMD kernels of the target architecture
+  when `xsimd <https://github.com/xtensor-stack/xsimd>`_ is found, otherwise
+  only the scalar kernel
+* ``OFF`` - build only the scalar kernel
+
+The SIMD kernels are ``avx512``, ``avx2`` and ``sse2`` on x86-64 and ``neon``
+on arm64. On first use the library picks the best kernel that the CPU
+supports, for example ``avx2`` on a CPU without AVX-512.
 
 Embed libdiffpy in another CMake project
 ----------------------------------------

news/faster-pdf.rst

@@ -0,0 +1,59 @@
+**Added:**
+
+* ``diffpy/hash.hpp`` with ``hash_combine``, ``hash_range``, ``hash_values``
+  and ``PairHash``.
+* SIMD calculation kernels for ``PDFCalculator`` and ``DebyePDFCalculator``:
+  AVX-512, AVX2 and SSE2 on x86-64, NEON on arm64. They are built when the
+  header-only xsimd library is found, unless the ``LIBDIFFPY_SIMD`` CMake
+  option is ``OFF``. The best kernel the CPU supports is used at runtime;
+  ``availableCalculationKernels``, ``getCalculationKernel`` and
+  ``setCalculationKernel`` select it by name.
+* ``PeakProfile::addRDFPeak`` to add a peak to an RDF array, and protected
+  ``GaussianProfile::kernelShape``, which derived profiles that keep
+  ``operator()`` can override to use the calculation kernels.
+
+**Changed:**
+
+* ``R3::Vector`` and ``R3::Matrix`` are fixed-size ``std::array`` types,
+  zero-initialized by default. Existing archives and pickles still load.
+* ``R3::row`` and ``R3::column`` return copies, so code such as
+  ``R3::row(A, 0) = v`` no longer compiles; use the new ``R3::setrow`` and
+  ``R3::setcolumn``.
+* ``R3::inverse``, ``R3::floor``, ``R3::mxvecproduct`` and the ``Lattice``
+  methods ``cartesian``, ``fractional``, ``ucvCartesian``, ``ucvFractional``,
+  ``cartesianMatrix`` and ``fractionalMatrix`` return by value instead of a
+  reference to a shared static result, so later calls no longer overwrite
+  earlier results. They are thread-safe, as are ``Lattice::distance`` and
+  ``Lattice::ucMaxDiagonal``.
+* Faster Gaussian peak evaluation and bond generation for periodic and
+  crystal structures. ``CrystalStructureAdapter`` expands symmetry faster
+  and only after its atoms, lattice or symmetry operations change.
+* Faster ``fftgtof`` and ``fftftog`` using cached mixed-radix complex FFTs.
+* Optimized builds of the shared library use link-time optimization unless
+  configured with ``-DLIBDIFFPY_LTO=OFF``. The library is compiled with
+  ``-fno-semantic-interposition`` when supported, so ``LD_PRELOAD`` may no
+  longer replace calls inside the library.
+
+**Removed:**
+
+* Dependencies on Boost.Functional/Hash, Boost.StringAlgo and Boost.uBLAS.
+  ``R3::Vector`` and ``R3::Matrix`` no longer accept uBLAS expressions.
+* Protected ``updater1`` method and ``msphere``, ``mrcsphere``, ``msymidx``,
+  ``mcstructure`` and ``mpuc1`` members of ``PeriodicStructureBondGenerator``
+  and ``CrystalStructureBondGenerator``, which no longer override
+  ``iterateSymmetry``.
+
+**Fixed:**
+
+* ``R3::inverse`` throws ``invalid_argument`` for a singular matrix instead
+  of aborting the program.
+* ``Lattice::setLatPar`` and ``Lattice::setLatBase`` leave the lattice
+  unchanged when they throw for degenerate parameters or base vectors.
+* ``CroppedGaussianProfile`` peaks of zero width no longer produce NaN in
+  ``PDFCalculator`` results.
+* Better accuracy for narrow Gaussian peaks at large pair distances.
+* Calculators give the correct result on the first evaluation after
+  a ``CrystalStructureAdapter`` atom moves in place to or from a special
+  position.
+* After an evaluation fails with an exception, the ``OPTIMIZED`` evaluator
+  recalculates fully instead of updating the incomplete result.

requirements/build.txt

@@ -1,3 +1,4 @@
 cxx-compiler
 cmake>=3.22
 ninja
+xsimd

src/diffpy/hash.hpp

@@ -0,0 +1,86 @@
+/*****************************************************************************
+ *
+ * libdiffpy         by DANSE Diffraction group
+ *                   Simon J. L. Billinge
+ *                   (c) 2026 The Trustees of Columbia University
+ *                   in the City of New York.  All rights reserved.
+ *
+ * See AUTHORS.rst for a list of people who contributed.
+ * See LICENSE_DANSE.rst for license information.
+ *
+ ******************************************************************************
+ *
+ * Hash utilities that replace Boost.Functional/Hash:
+ *
+ *   diffpy::hash_of      -- hash_value(v) found by ADL, else std::hash<T>
+ *   diffpy::hash_combine -- mix the hash of a value into a seed
+ *   diffpy::hash_range   -- hash of a sequence of values
+ *   diffpy::hash_values  -- hash of several values
+ *   diffpy::PairHash     -- hash function object for std::pair
+ *
+ *****************************************************************************/
+
+#ifndef HASH_HPP_INCLUDED
+#define HASH_HPP_INCLUDED
+
+#include <concepts>
+#include <cstddef>
+#include <functional>
+#include <utility>
+
+namespace diffpy {
+
+namespace detail {
+
+/// True when hash_value(v) can be found by argument-dependent lookup.
+template <class T>
+concept HasHashValue = requires(const T& v) {
+  { hash_value(v) } -> std::convertible_to<std::size_t>;
+};
+
+}  // namespace detail
+
+/// Hash of v from hash_value(v) when available, otherwise std::hash<T>.
+template <class T>
+std::size_t hash_of(const T& v) {
+  if constexpr (detail::HasHashValue<T>) {
+    return hash_value(v);
+  } else {
+    return std::hash<T>{}(v);
+  }
+}
+
+/// Mix the hash of v into seed, same scheme as boost::hash_combine.
+template <class T>
+void hash_combine(std::size_t& seed, const T& v) {
+  const std::size_t magic = static_cast<std::size_t>(0x9e3779b97f4a7c15ull);
+  seed ^= hash_of(v) + magic + (seed << 6) + (seed >> 2);
+}
+
+/// Order-dependent hash of the values in [first, last).
+template <class It>
+std::size_t hash_range(It first, It last) {
+  std::size_t seed = 0;
+  for (; first != last; ++first) hash_combine(seed, *first);
+  return seed;
+}
+
+/// Order-dependent hash of the arguments, same as hash_combine of each.
+template <class... Ts>
+std::size_t hash_values(const Ts&... values) {
+  std::size_t seed = 0;
+  (hash_combine(seed, values), ...);
+  return seed;
+}
+
+/// Hash function object for std::pair, for use in unordered containers.
+struct PairHash {
+  template <class A, class B>
+  std::size_t operator()(const std::pair<A, B>& p) const {
+    return hash_values(p.first, p.second);
+  }
+};
+
+}  // namespace diffpy
+
+#endif  // HASH_HPP_INCLUDED

src/diffpy/srreal/AtomicStructureAdapter.cpp

@@ -19,8 +19,8 @@
 
 #include <cassert>
 #include <algorithm>
-#include <boost/functional/hash.hpp>
 
+#include <diffpy/hash.hpp>
 #include <diffpy/serialization.ipp>
 #include <diffpy/srreal/AtomicStructureAdapter.hpp>
 #include <diffpy/srreal/StructureDifference.hpp>
@@ -69,13 +69,8 @@ bool operator==(const Atom& a0, const Atom& a1) {
 bool operator!=(const Atom& a0, const Atom& a1) { return !(a0 == a1); }
 
 size_t hash_value(const Atom& a) {
-  size_t seed = 0;
-  boost::hash_combine(seed, a.atomtype);
-  boost::hash_combine(seed, a.xyz_cartn);
-  boost::hash_combine(seed, a.occupancy);
-  boost::hash_combine(seed, a.anisotropy);
-  boost::hash_combine(seed, a.uij_cartn);
-  return seed;
+  return diffpy::hash_values(a.atomtype, a.xyz_cartn, a.occupancy, a.anisotropy,
+                             a.uij_cartn);
 }
 
 //////////////////////////////////////////////////////////////////////////////