From 73981a61e71387eac48b7f70193251f44e13658d Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Tue, 10 Nov 2026 21:48:16 -0800 Subject: [PATCH 1/6] refactor: replace Boost hash and uBLAS with fixed-size types --- news/faster-pdf.rst | 30 ++ src/diffpy/hash.hpp | 86 +++++ src/diffpy/srreal/AtomicStructureAdapter.cpp | 11 +- src/diffpy/srreal/BVParam.cpp | 9 +- src/diffpy/srreal/BVSCalculator.cpp | 4 +- src/diffpy/srreal/Lattice.cpp | 118 +++---- src/diffpy/srreal/Lattice.hpp | 64 ++-- src/diffpy/srreal/PairQuantity.hpp | 7 +- src/diffpy/srreal/R3linalg.cpp | 114 ++++--- src/diffpy/srreal/R3linalg.hpp | 339 +++++++++++-------- src/diffpy/srreal/R3linalg.ipp | 102 ++++++ src/tests/TestHash.cpp | 68 ++++ src/tests/TestLattice.cpp | 33 ++ src/tests/TestR3linalg.cpp | 198 +++++++++++ src/tests/test_helpers.cpp | 26 +- 15 files changed, 909 insertions(+), 300 deletions(-) create mode 100644 news/faster-pdf.rst create mode 100644 src/diffpy/hash.hpp create mode 100644 src/diffpy/srreal/R3linalg.ipp create mode 100644 src/tests/TestHash.cpp diff --git a/news/faster-pdf.rst b/news/faster-pdf.rst new file mode 100644 index 0000000..c9f32d1 --- /dev/null +++ b/news/faster-pdf.rst @@ -0,0 +1,30 @@ +**Added:** + +* ``diffpy/hash.hpp`` with ``hash_combine``, ``hash_range``, ``hash_values`` + and ``PairHash``. + +**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``. + +**Removed:** + +* Dependency on Boost.Functional/Hash, Boost.StringAlgo and Boost.uBLAS. + ``R3::Vector`` and ``R3::Matrix`` no longer accept uBLAS expressions. + +**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. diff --git a/src/diffpy/hash.hpp b/src/diffpy/hash.hpp new file mode 100644 index 0000000..e64ef67 --- /dev/null +++ b/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 + * 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 +#include +#include +#include + +namespace diffpy { + +namespace detail { + +/// True when hash_value(v) can be found by argument-dependent lookup. +template +concept HasHashValue = requires(const T& v) { + { hash_value(v) } -> std::convertible_to; +}; + +} // namespace detail + +/// Hash of v from hash_value(v) when available, otherwise std::hash. +template +std::size_t hash_of(const T& v) { + if constexpr (detail::HasHashValue) { + return hash_value(v); + } else { + return std::hash{}(v); + } +} + +/// Mix the hash of v into seed, same scheme as boost::hash_combine. +template +void hash_combine(std::size_t& seed, const T& v) { + const std::size_t magic = static_cast(0x9e3779b97f4a7c15ull); + seed ^= hash_of(v) + magic + (seed << 6) + (seed >> 2); +} + +/// Order-dependent hash of the values in [first, last). +template +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 +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 + std::size_t operator()(const std::pair& p) const { + return hash_values(p.first, p.second); + } +}; + +} // namespace diffpy + +#endif // HASH_HPP_INCLUDED diff --git a/src/diffpy/srreal/AtomicStructureAdapter.cpp b/src/diffpy/srreal/AtomicStructureAdapter.cpp index 967bdd0..dec63d4 100644 --- a/src/diffpy/srreal/AtomicStructureAdapter.cpp +++ b/src/diffpy/srreal/AtomicStructureAdapter.cpp @@ -19,8 +19,8 @@ #include #include -#include +#include #include #include #include @@ -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); } ////////////////////////////////////////////////////////////////////////////// diff --git a/src/diffpy/srreal/BVParam.cpp b/src/diffpy/srreal/BVParam.cpp index 90ae2eb..15bec30 100644 --- a/src/diffpy/srreal/BVParam.cpp +++ b/src/diffpy/srreal/BVParam.cpp @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include #include @@ -92,12 +92,7 @@ void BVParam::setFromCifLine(const std::string& cifline) { // class BVParam::BondHash --------------------------------------------------- size_t BVParam::BondHash::operator()(const BVParam& bp) const { - size_t seed = 0; - boost::hash_combine(seed, bp.matom0); - boost::hash_combine(seed, bp.mvalence0); - boost::hash_combine(seed, bp.matom1); - boost::hash_combine(seed, bp.mvalence1); - return seed; + return diffpy::hash_values(bp.matom0, bp.mvalence0, bp.matom1, bp.mvalence1); } // class BVParam::BondEqual -------------------------------------------------- diff --git a/src/diffpy/srreal/BVSCalculator.cpp b/src/diffpy/srreal/BVSCalculator.cpp index c4ec1e1..63228d0 100644 --- a/src/diffpy/srreal/BVSCalculator.cpp +++ b/src/diffpy/srreal/BVSCalculator.cpp @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -164,8 +165,7 @@ void BVSCalculator::cacheStructureData() { double BVSCalculator::rmaxFromPrecision(double eps) const { const BVParametersTable& bvtb = *(this->getBVParamTable()); // build a set of unique (symbol, valence) pairs - typedef std::unordered_set, - boost::hash > > + typedef std::unordered_set, diffpy::PairHash> SymbolValenceSet; SymbolValenceSet allsymvals; assert(int(mstructure_cache.baresymbols.size()) == this->countSites()); diff --git a/src/diffpy/srreal/Lattice.cpp b/src/diffpy/srreal/Lattice.cpp index 9cf9e93..447a8d7 100644 --- a/src/diffpy/srreal/Lattice.cpp +++ b/src/diffpy/srreal/Lattice.cpp @@ -17,7 +17,8 @@ * *****************************************************************************/ -#include +#include +#include #include #include @@ -33,19 +34,35 @@ using namespace diffpy::srreal; Lattice::Lattice() { mbaserot = R3::identity(); - this->setLatPar(1.0, 1.0, 1.0, 90.0, 90.0, 90.0); + this->assignLatPar(1.0, 1.0, 1.0, 90.0, 90.0, 90.0); } Lattice::Lattice(double a0, double b0, double c0, double alpha0, double beta0, double gamma0) { mbaserot = R3::identity(); - this->setLatPar(a0, b0, c0, alpha0, beta0, gamma0); + this->assignLatPar(a0, b0, c0, alpha0, beta0, gamma0); } // Public Methods ------------------------------------------------------------ +// The setters leave this lattice unchanged when they throw. + void Lattice::setLatPar(double a0, double b0, double c0, double alpha0, double beta0, double gamma0) { + Lattice L(*this); + L.assignLatPar(a0, b0, c0, alpha0, beta0, gamma0); + *this = L; +} + +void Lattice::setLatBase(const R3::Vector& va0, const R3::Vector& vb0, + const R3::Vector& vc0) { + Lattice L(*this); + L.assignLatBase(va0, vb0, vc0); + *this = L; +} + +void Lattice::assignLatPar(double a0, double b0, double c0, double alpha0, + double beta0, double gamma0) { using diffpy::mathutils::acosd; using diffpy::mathutils::cosd; using diffpy::mathutils::sind; @@ -89,20 +106,20 @@ void Lattice::setLatPar(double a0, double b0, double c0, double alpha0, mvar = R3::column(mrecbase, 0); mvbr = R3::column(mrecbase, 1); mvcr = R3::column(mrecbase, 2); - R3::row(mnormbase, 0) = R3::row(mbase, 0) * mar; - R3::row(mnormbase, 1) = R3::row(mbase, 1) * mbr; - R3::row(mnormbase, 2) = R3::row(mbase, 2) * mcr; - R3::column(mrecnormbase, 0) = R3::column(mrecbase, 0) / mar; - R3::column(mrecnormbase, 1) = R3::column(mrecbase, 1) / mbr; - R3::column(mrecnormbase, 2) = R3::column(mrecbase, 2) / mcr; + R3::setrow(mnormbase, 0, R3::row(mbase, 0) * mar); + R3::setrow(mnormbase, 1, R3::row(mbase, 1) * mbr); + R3::setrow(mnormbase, 2, R3::row(mbase, 2) * mcr); + R3::setcolumn(mrecnormbase, 0, R3::column(mrecbase, 0) / mar); + R3::setcolumn(mrecnormbase, 1, R3::column(mrecbase, 1) / mbr); + R3::setcolumn(mrecnormbase, 2, R3::column(mrecbase, 2) / mcr); } -void Lattice::setLatBase(const R3::Vector& va0, const R3::Vector& vb0, - const R3::Vector& vc0) { +void Lattice::assignLatBase(const R3::Vector& va0, const R3::Vector& vb0, + const R3::Vector& vc0) { using diffpy::mathutils::acosd; - R3::row(mbase, 0) = va0; - R3::row(mbase, 1) = vb0; - R3::row(mbase, 2) = vc0; + R3::setrow(mbase, 0, va0); + R3::setrow(mbase, 1, vb0); + R3::setrow(mbase, 2, vc0); double detbase = R3::determinant(mbase); if (fabs(detbase) < 1.0e-8) { throw invalid_argument("base vectors are degenerate"); @@ -148,12 +165,12 @@ void Lattice::setLatBase(const R3::Vector& va0, const R3::Vector& vb0, mvbr = R3::column(mrecbase, 1); mvcr = R3::column(mrecbase, 2); // bases normalized to unit reciprocal vectors - R3::row(mnormbase, 0) = R3::row(mbase, 0) * mar; - R3::row(mnormbase, 1) = R3::row(mbase, 1) * mbr; - R3::row(mnormbase, 2) = R3::row(mbase, 2) * mcr; - R3::column(mrecnormbase, 0) = R3::column(mrecbase, 0) / mar; - R3::column(mrecnormbase, 1) = R3::column(mrecbase, 1) / mbr; - R3::column(mrecnormbase, 2) = R3::column(mrecbase, 2) / mcr; + R3::setrow(mnormbase, 0, R3::row(mbase, 0) * mar); + R3::setrow(mnormbase, 1, R3::row(mbase, 1) * mbr); + R3::setrow(mnormbase, 2, R3::row(mbase, 2) * mcr); + R3::setcolumn(mrecnormbase, 0, R3::column(mrecbase, 0) / mar); + R3::setcolumn(mrecnormbase, 1, R3::column(mrecbase, 1) / mbr); + R3::setcolumn(mrecnormbase, 2, R3::column(mrecbase, 2) / mcr); this->updateMetrics(); } @@ -169,67 +186,52 @@ double Lattice::volume() const { return rv; } -const R3::Vector& Lattice::cartesian(const R3::Vector& lv) const { - static R3::Vector res; - res = R3::mxvecproduct(lv, mbase); - return res; +R3::Vector Lattice::cartesian(const R3::Vector& lv) const { + return R3::mxvecproduct(lv, mbase); } -const R3::Vector& Lattice::fractional(const R3::Vector& cv) const { - static R3::Vector res; - res = R3::mxvecproduct(cv, mrecbase); - return res; +R3::Vector Lattice::fractional(const R3::Vector& cv) const { + return R3::mxvecproduct(cv, mrecbase); } -const R3::Vector& Lattice::ucvCartesian(const R3::Vector& cv) const { - static R3::Vector res; - res = cartesian(ucvFractional(fractional(cv))); - return res; +R3::Vector Lattice::ucvCartesian(const R3::Vector& cv) const { + return cartesian(ucvFractional(fractional(cv))); } -const R3::Vector& Lattice::ucvFractional(const R3::Vector& lv) const { +R3::Vector Lattice::ucvFractional(const R3::Vector& lv) const { using mathutils::eps_eq; - static R3::Vector res; - res = lv - floor(lv); + R3::Vector res = lv - floor(lv); if (eps_eq(res[0], 1.0)) res[0] = 0.0; if (eps_eq(res[1], 1.0)) res[1] = 0.0; if (eps_eq(res[2], 1.0)) res[2] = 0.0; return res; } -const R3::Matrix& Lattice::cartesianMatrix(const R3::Matrix& Ml) const { - static R3::Matrix res0, res1; - res0 = prod(Ml, mnormbase); - res1 = prod(R3::trans(mnormbase), res0); - return res1; +R3::Matrix Lattice::cartesianMatrix(const R3::Matrix& Ml) const { + return R3::prod(R3::trans(mnormbase), R3::prod(Ml, mnormbase)); } -const R3::Matrix& Lattice::fractionalMatrix(const R3::Matrix& Mc) const { - static R3::Matrix res0, res1; - res0 = prod(Mc, mrecnormbase); - res1 = prod(R3::trans(mrecnormbase), res0); - return res1; +R3::Matrix Lattice::fractionalMatrix(const R3::Matrix& Mc) const { + return R3::prod(R3::trans(mrecnormbase), R3::prod(Mc, mrecnormbase)); } const R3::Vector& Lattice::ucMaxDiagonal() const { - static list ucdiagonals; - if (ucdiagonals.empty()) { - ucdiagonals.push_back(R3::Vector(+1, +1, +1)); - ucdiagonals.push_back(R3::Vector(-1, +1, +1)); - ucdiagonals.push_back(R3::Vector(+1, -1, +1)); - ucdiagonals.push_back(R3::Vector(+1, +1, -1)); - } + static constexpr array ucdiagonals = {{ + {+1, +1, +1}, + {-1, +1, +1}, + {+1, -1, +1}, + {+1, +1, -1}, + }}; double maxnorm = -1; - list::iterator ucd; - list::iterator maxucd = ucdiagonals.end(); - for (ucd = ucdiagonals.begin(); ucd != ucdiagonals.end(); ++ucd) { - double normucd = this->norm(*ucd); + const R3::Vector* maxucd = nullptr; + for (const R3::Vector& ucd : ucdiagonals) { + double normucd = this->norm(ucd); if (normucd > maxnorm) { maxnorm = normucd; - maxucd = ucd; + maxucd = &ucd; } } - assert(maxucd != ucdiagonals.end()); + assert(maxucd); return *maxucd; } diff --git a/src/diffpy/srreal/Lattice.hpp b/src/diffpy/srreal/Lattice.hpp index 54616d7..d075869 100644 --- a/src/diffpy/srreal/Lattice.hpp +++ b/src/diffpy/srreal/Lattice.hpp @@ -108,20 +108,20 @@ class DLL_EXPORT Lattice { template double anglerad(const V& u, const V& v) const; // conversion of coordinates and tensors - const R3::Vector& cartesian(const R3::Vector& lv) const; + R3::Vector cartesian(const R3::Vector& lv) const; template - const R3::Vector& cartesian(const V& lv) const; - const R3::Vector& fractional(const R3::Vector& cv) const; + R3::Vector cartesian(const V& lv) const; + R3::Vector fractional(const R3::Vector& cv) const; template - const R3::Vector& fractional(const V& cv) const; - const R3::Vector& ucvCartesian(const R3::Vector& cv) const; + R3::Vector fractional(const V& cv) const; + R3::Vector ucvCartesian(const R3::Vector& cv) const; template - const R3::Vector& ucvCartesian(const V& cv) const; - const R3::Vector& ucvFractional(const R3::Vector& lv) const; + R3::Vector ucvCartesian(const V& cv) const; + R3::Vector ucvFractional(const R3::Vector& lv) const; template - const R3::Vector& ucvFractional(const V& lv) const; - const R3::Matrix& cartesianMatrix(const R3::Matrix& Ml) const; - const R3::Matrix& fractionalMatrix(const R3::Matrix& Mc) const; + R3::Vector ucvFractional(const V& lv) const; + R3::Matrix cartesianMatrix(const R3::Matrix& Ml) const; + R3::Matrix fractionalMatrix(const R3::Matrix& Mc) const; // largest cell diagonal in fractional coordinates const R3::Vector& ucMaxDiagonal() const; double ucMaxDiagonalLength() const; @@ -132,6 +132,11 @@ class DLL_EXPORT Lattice { private: // methods + /// setLatPar and setLatBase in place, not exception safe + void assignLatPar(double a, double b, double c, double alpha, double beta, + double gamma); + void assignLatBase(const R3::Vector& va, const R3::Vector& vb, + const R3::Vector& vc); void updateMetrics(); void updateStandardBase(); @@ -191,7 +196,7 @@ class DLL_EXPORT Lattice { // Template Constructor ------------------------------------------------------ template -Lattice::Lattice(const V& va0, const V& vb0, const V& vc0) { +Lattice::Lattice(const V& va0, const V& vb0, const V& vc0) : Lattice() { this->setLatBase(va0, vb0, vc0); } @@ -298,10 +303,7 @@ double Lattice::norm(const V& u) const { template double Lattice::distance(const V& u, const V& v) const { - static R3::Vector duv; - duv[0] = u[0] - v[0]; - duv[1] = u[1] - v[1]; - duv[2] = u[2] - v[2]; + const R3::Vector duv(u[0] - v[0], u[1] - v[1], u[2] - v[2]); return this->norm(duv); } @@ -319,39 +321,23 @@ double Lattice::anglerad(const V& u, const V& v) const { } template -const R3::Vector& Lattice::cartesian(const V& lv) const { - static R3::Vector lvcopy; - lvcopy[0] = lv[0]; - lvcopy[1] = lv[1]; - lvcopy[2] = lv[2]; - return this->cartesian(lvcopy); +R3::Vector Lattice::cartesian(const V& lv) const { + return this->cartesian(R3::Vector(lv[0], lv[1], lv[2])); } template -const R3::Vector& Lattice::fractional(const V& cv) const { - static R3::Vector cvcopy; - cvcopy[0] = cv[0]; - cvcopy[1] = cv[1]; - cvcopy[2] = cv[2]; - return this->fractional(cvcopy); +R3::Vector Lattice::fractional(const V& cv) const { + return this->fractional(R3::Vector(cv[0], cv[1], cv[2])); } template -const R3::Vector& Lattice::ucvCartesian(const V& cv) const { - static R3::Vector cvcopy; - cvcopy[0] = cv[0]; - cvcopy[1] = cv[1]; - cvcopy[2] = cv[2]; - return this->ucvCartesian(cvcopy); +R3::Vector Lattice::ucvCartesian(const V& cv) const { + return this->ucvCartesian(R3::Vector(cv[0], cv[1], cv[2])); } template -const R3::Vector& Lattice::ucvFractional(const V& cv) const { - static R3::Vector cvcopy; - cvcopy[0] = cv[0]; - cvcopy[1] = cv[1]; - cvcopy[2] = cv[2]; - return ucvFractional(cvcopy); +R3::Vector Lattice::ucvFractional(const V& lv) const { + return this->ucvFractional(R3::Vector(lv[0], lv[1], lv[2])); } } // namespace srreal diff --git a/src/diffpy/srreal/PairQuantity.hpp b/src/diffpy/srreal/PairQuantity.hpp index 988c51d..ee8fb51 100644 --- a/src/diffpy/srreal/PairQuantity.hpp +++ b/src/diffpy/srreal/PairQuantity.hpp @@ -24,8 +24,8 @@ #include #include #include -#include +#include #include #include #include @@ -104,11 +104,10 @@ class DLL_EXPORT PairQuantity : public diffpy::Attributes { virtual void restorePartialValue(); // data - typedef std::unordered_set, - boost::hash > > + typedef std::unordered_set, diffpy::PairHash> PairMaskStorage; typedef std::unordered_map, bool, - boost::hash > > + diffpy::PairHash> TypeMaskStorage; QuantityType mvalue; StructureAdapterPtr mstructure; diff --git a/src/diffpy/srreal/R3linalg.cpp b/src/diffpy/srreal/R3linalg.cpp index 8ba3043..7fb27ee 100644 --- a/src/diffpy/srreal/R3linalg.cpp +++ b/src/diffpy/srreal/R3linalg.cpp @@ -16,9 +16,12 @@ * *****************************************************************************/ -#include -#include -#include +#include +#include +#include +#include + +#include #include namespace diffpy { @@ -30,56 +33,90 @@ DLL_EXPORT const int Ndim = 3; // Functions const Matrix& identity() { - static Matrix mx = ublas::identity_matrix(Ndim); + static const Matrix mx(1, 0, 0, 0, 1, 0, 0, 0, 1); return mx; } const Matrix& zeromatrix() { - static Matrix mx = ublas::zero_matrix(Ndim, Ndim); + static const Matrix mx; return mx; } -double determinant(const Matrix& A) { - gsl_matrix* gA = gsl_matrix_alloc(Ndim, Ndim); - for (int i = 0; i != Ndim; ++i) { - for (int j = 0; j != Ndim; ++j) { - gsl_matrix_set(gA, i, j, A(i, j)); +namespace { + +/// LU decomposition with partial pivoting, P * A = L * U. Unlike cofactor +/// expansion it stays accurate for nearly singular matrices. +struct LUDecomposition { + /// U on and above the diagonal, unit lower-triangular L below it. + double lu[3][3]; + /// Row i of P * A is row perm[i] of A. + int perm[3] = {0, 1, 2}; + /// Odd number of row exchanges, which flips the sign of determinant. + bool odd = false; + /// Some pivot is exactly zero. + bool singular = false; + + explicit LUDecomposition(const Matrix& A) { + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) lu[i][j] = A(i, j); + } + for (int k = 0; k < 3; ++k) { + int p = k; + for (int i = k + 1; i < 3; ++i) { + if (std::fabs(lu[i][k]) > std::fabs(lu[p][k])) p = i; + } + if (p != k) { + std::swap(lu[p], lu[k]); + std::swap(perm[p], perm[k]); + odd = !odd; + } + // zero pivot means the rest of the column is zero as well + if (lu[k][k] == 0.0) { + singular = true; + continue; + } + for (int i = k + 1; i < 3; ++i) { + const double f = (lu[i][k] /= lu[k][k]); + for (int j = k + 1; j < 3; ++j) lu[i][j] -= f * lu[k][j]; + } } } - gsl_permutation* gP = gsl_permutation_alloc(Ndim); - int signum; - gsl_linalg_LU_decomp(gA, gP, &signum); - double det = gsl_linalg_LU_det(gA, signum); - gsl_permutation_free(gP); - gsl_matrix_free(gA); - return det; +}; + +} // namespace + +double determinant(const Matrix& A) { + const LUDecomposition d(A); + const double det = d.lu[0][0] * d.lu[1][1] * d.lu[2][2]; + return d.odd ? -det : det; } -const Matrix& inverse(const Matrix& A) { - static Matrix B; - gsl_matrix* gA = gsl_matrix_alloc(Ndim, Ndim); - for (int i = 0; i != Ndim; ++i) { - for (int j = 0; j != Ndim; ++j) { - gsl_matrix_set(gA, i, j, A(i, j)); +Matrix inverse(const Matrix& A) { + const LUDecomposition d(A); + if (d.singular) throw std::invalid_argument("singular matrix"); + Matrix B; + // solve L * U * x = P * e_j for every column j of the identity matrix + for (int j = 0; j < 3; ++j) { + double x[3]; + for (int i = 0; i < 3; ++i) { + x[i] = (d.perm[i] == j) ? 1.0 : 0.0; + for (int k = 0; k < i; ++k) x[i] -= d.lu[i][k] * x[k]; + } + for (int i = 2; i >= 0; --i) { + for (int k = i + 1; k < 3; ++k) x[i] -= d.lu[i][k] * x[k]; + x[i] /= d.lu[i][i]; + B(i, j) = x[i]; } } - gsl_permutation* gP = gsl_permutation_alloc(Ndim); - int signum; - gsl_linalg_LU_decomp(gA, gP, &signum); - double* bdata = &(B.data()[0]); - gsl_matrix_view gB = gsl_matrix_view_array(bdata, Ndim, Ndim); - gsl_linalg_LU_invert(gA, gP, &gB.matrix); - gsl_permutation_free(gP); - gsl_matrix_free(gA); return B; } size_t hash_value(const Vector& v) { - return boost::hash_range(v.begin(), v.end()); + return diffpy::hash_range(v.begin(), v.end()); } size_t hash_value(const Matrix& A) { - return boost::hash_range(A.data().begin(), A.data().end()); + return diffpy::hash_range(A.data().begin(), A.data().end()); } } // namespace R3 @@ -92,9 +129,7 @@ namespace mathutils { template <> bool EpsilonLess::operator()( const srreal::R3::Matrix& A, const srreal::R3::Matrix& B) const { - bool rv = std::lexicographical_compare( - A.data().begin(), A.data().end(), B.data().begin(), B.data().end(), *this); - return rv; + return std::ranges::lexicographical_compare(A.data(), B.data(), *this); } // EpsilonEqual specializations ---------------------------------------------- @@ -102,16 +137,13 @@ bool EpsilonLess::operator()( template <> bool EpsilonEqual::operator()( const srreal::R3::Vector& u, const srreal::R3::Vector& v) const { - bool rv = std::equal(u.begin(), u.end(), v.begin(), *this); - return rv; + return std::ranges::equal(u, v, *this); } template <> bool EpsilonEqual::operator()( const srreal::R3::Matrix& A, const srreal::R3::Matrix& B) const { - bool rv = - std::equal(A.data().begin(), A.data().end(), B.data().begin(), *this); - return rv; + return std::ranges::equal(A.data(), B.data(), *this); } } // namespace mathutils diff --git a/src/diffpy/srreal/R3linalg.hpp b/src/diffpy/srreal/R3linalg.hpp index 5d51b80..01728df 100644 --- a/src/diffpy/srreal/R3linalg.hpp +++ b/src/diffpy/srreal/R3linalg.hpp @@ -20,10 +20,10 @@ #define R3LINALG_HPP_INCLUDED #include -#include -#include -#include -#include +#include +#include +#include +#include #include #include @@ -32,209 +32,277 @@ namespace diffpy { namespace srreal { namespace R3 { -// Declarations -------------------------------------------------------------- - -namespace ublas = boost::numeric::ublas; -using ublas::column; -using ublas::prod; -using ublas::row; -using ublas::trans; - -// Constants +// Constants ----------------------------------------------------------------- extern DLL_EXPORT const int Ndim; using ::diffpy::mathutils::SQRT_DOUBLE_EPS; -const ublas::zero_vector zerovector(Ndim); -// Classes - -class DLL_EXPORT Vector - : public ublas::vector > { - typedef ublas::vector > BaseVector; +// Classes ------------------------------------------------------------------- +class Vector { public: - // constructors - Vector() : BaseVector(3) {} + typedef double value_type; + typedef std::size_t size_type; + typedef std::array array_type; + typedef array_type::iterator iterator; + typedef array_type::const_iterator const_iterator; - Vector(const double& x, const double& y, const double& z) : BaseVector(3) { - Vector::iterator xi = this->begin(); - *(xi++) = x; - *(xi++) = y; - *(xi++) = z; + // constructors + constexpr Vector() noexcept : mdata{} {} + constexpr Vector(double x, double y, double z) noexcept : mdata{x, y, z} {} + constexpr Vector(const Vector&) noexcept = default; + constexpr Vector(Vector&&) noexcept = default; + + // Assignments are restricted to lvalues so that writes to a temporary, + // such as row(A, 0) = v, fail to compile instead of being lost. + constexpr Vector& operator=(const Vector&) & noexcept = default; + constexpr Vector& operator=(Vector&&) & noexcept = default; + + // Element writes need an lvalue too, row(A, 0)[1] = x does not compile. + constexpr double& operator[](size_type i) & noexcept { return mdata[i]; } + constexpr const double& operator[](size_type i) const& noexcept { + return mdata[i]; + } + constexpr double& operator()(size_type i) & noexcept { return mdata[i]; } + constexpr const double& operator()(size_type i) const& noexcept { + return mdata[i]; + } + static constexpr size_type size() noexcept { return 3; } + constexpr array_type& data() & noexcept { return mdata; } + constexpr const array_type& data() const& noexcept { return mdata; } + + // iterators, writable only for lvalues + constexpr iterator begin() & noexcept { return mdata.begin(); } + constexpr iterator end() & noexcept { return mdata.end(); } + constexpr const_iterator begin() const& noexcept { return mdata.begin(); } + constexpr const_iterator end() const& noexcept { return mdata.end(); } + + // exact element-wise equality, != is rewritten from it + constexpr bool operator==(const Vector&) const noexcept = default; + + // in-place arithmetic + constexpr Vector& operator+=(const Vector& v) & noexcept { + for (size_type k = 0; k < 3; ++k) mdata[k] += v.mdata[k]; + return *this; } - template - Vector(const ublas::vector_expression& r) : BaseVector(r) {} + constexpr Vector& operator-=(const Vector& v) & noexcept { + for (size_type k = 0; k < 3; ++k) mdata[k] -= v.mdata[k]; + return *this; + } - template - void operator=(const ublas::vector_expression& r) { - BaseVector::operator=(r); + constexpr Vector& operator*=(double s) & noexcept { + for (double& x : mdata) x *= s; + return *this; } - template - void operator=(const BaseVector& r) { - BaseVector::operator=(r); + constexpr Vector& operator/=(double s) & noexcept { + for (double& x : mdata) x /= s; + return *this; } private: + array_type mdata; + // serialization friend class boost::serialization::access; template - void serialize(Archive& ar, const unsigned int version) { - using boost::serialization::base_object; - ar& base_object(*this); - } + void serialize(Archive& ar, const unsigned int version); }; -class DLL_EXPORT Matrix - : public ublas::matrix > { - typedef ublas::matrix > - BaseMatrix; - +class Matrix { public: + typedef double value_type; + typedef std::size_t size_type; + typedef std::array array_type; + // constructors - Matrix() : BaseMatrix(3, 3) {} - - Matrix(const double& x0, const double& x1, const double& x2, const double& x3, - const double& x4, const double& x5, const double& x6, const double& x7, - const double& x8) - : BaseMatrix(3, 3) { - Matrix::array_type::iterator xi = this->data().begin(); - *(xi++) = x0; - *(xi++) = x1; - *(xi++) = x2; - *(xi++) = x3; - *(xi++) = x4; - *(xi++) = x5; - *(xi++) = x6; - *(xi++) = x7; - *(xi++) = x8; + constexpr Matrix() noexcept : mdata{} {} + constexpr Matrix(double x0, double x1, double x2, double x3, double x4, + double x5, double x6, double x7, double x8) noexcept + : mdata{x0, x1, x2, x3, x4, x5, x6, x7, x8} {} + constexpr Matrix(const Matrix&) noexcept = default; + constexpr Matrix(Matrix&&) noexcept = default; + + // assignments are restricted to lvalues, see Vector + constexpr Matrix& operator=(const Matrix&) & noexcept = default; + constexpr Matrix& operator=(Matrix&&) & noexcept = default; + + // row-major element access, writable only for lvalues + constexpr double& operator()(size_type i, size_type j) & noexcept { + return mdata[3 * i + j]; + } + constexpr const double& operator()(size_type i, size_type j) const& noexcept { + return mdata[3 * i + j]; + } + static constexpr size_type size1() noexcept { return 3; } + static constexpr size_type size2() noexcept { return 3; } + constexpr array_type& data() & noexcept { return mdata; } + constexpr const array_type& data() const& noexcept { return mdata; } + + // exact element-wise equality, != is rewritten from it + constexpr bool operator==(const Matrix&) const noexcept = default; + + // in-place arithmetic + constexpr Matrix& operator+=(const Matrix& A) & noexcept { + for (size_type k = 0; k < 9; ++k) mdata[k] += A.mdata[k]; + return *this; } - template - Matrix(const ublas::matrix_expression& r) : BaseMatrix(r) {} + constexpr Matrix& operator-=(const Matrix& A) & noexcept { + for (size_type k = 0; k < 9; ++k) mdata[k] -= A.mdata[k]; + return *this; + } - template - void operator=(const ublas::matrix_expression& r) { - BaseMatrix::operator=(r); + constexpr Matrix& operator*=(double s) & noexcept { + for (double& x : mdata) x *= s; + return *this; } - template - void operator=(const BaseMatrix& r) { - BaseMatrix::operator=(r); + constexpr Matrix& operator/=(double s) & noexcept { + for (double& x : mdata) x /= s; + return *this; } private: + array_type mdata; + // serialization friend class boost::serialization::access; template - void serialize(Archive& ar, const unsigned int version) { - using boost::serialization::base_object; - ar& base_object(*this); - } + void serialize(Archive& ar, const unsigned int version); }; -// Functions +inline constexpr Vector zerovector{}; + +// Functions ----------------------------------------------------------------- DLL_EXPORT const Matrix& identity(); DLL_EXPORT const Matrix& zeromatrix(); DLL_EXPORT double determinant(const Matrix& A); -DLL_EXPORT const Matrix& inverse(const Matrix& A); +DLL_EXPORT Matrix inverse(const Matrix& A); -const Vector& floor(const Vector&); -template -double norm(const V&); -template -double distance(const V& u, const V& v); -template -double dot(const V& u, const V& v); -template -Vector cross(const V& u, const V& v); -template -const Vector& mxvecproduct(const Matrix&, const V&); -template -const Vector& mxvecproduct(const V&, const Matrix&); +// Vector arithmetic --------------------------------------------------------- -// Equality ------------------------------------------------------------------ +constexpr Vector operator-(const Vector& u) noexcept { + return Vector(-u[0], -u[1], -u[2]); +} + +constexpr Vector operator+(Vector u, const Vector& v) noexcept { + return u += v; +} -inline bool operator==(const Vector& u, const Vector& v) { - bool rv = std::equal(u.begin(), u.end(), v.begin()); - return rv; +constexpr Vector operator-(Vector u, const Vector& v) noexcept { + return u -= v; } -inline bool operator!=(const Vector& u, const Vector& v) { return !(u == v); } +constexpr Vector operator*(Vector u, double s) noexcept { return u *= s; } +constexpr Vector operator*(double s, Vector u) noexcept { return u *= s; } +constexpr Vector operator/(Vector u, double s) noexcept { return u /= s; } -inline bool operator==(const Matrix& A, const Matrix& B) { - bool rv = (&A == &B) || - std::equal(A.data().begin(), A.data().end(), B.data().begin()); - return rv; +// Matrix arithmetic --------------------------------------------------------- + +constexpr Matrix operator-(Matrix A) noexcept { return A *= -1.0; } + +constexpr Matrix operator+(Matrix A, const Matrix& B) noexcept { + return A += B; +} + +constexpr Matrix operator-(Matrix A, const Matrix& B) noexcept { + return A -= B; +} + +constexpr Matrix operator*(Matrix A, double s) noexcept { return A *= s; } +constexpr Matrix operator*(double s, Matrix A) noexcept { return A *= s; } +constexpr Matrix operator/(Matrix A, double s) noexcept { return A /= s; } + +/// Matrix product A * B. +constexpr Matrix prod(const Matrix& A, const Matrix& B) noexcept { + Matrix C; + for (std::size_t i = 0; i < 3; ++i) { + for (std::size_t j = 0; j < 3; ++j) { + C(i, j) = A(i, 0) * B(0, j) + A(i, 1) * B(1, j) + A(i, 2) * B(2, j); + } + } + return C; +} + +constexpr Matrix trans(const Matrix& A) noexcept { + return Matrix(A(0, 0), A(1, 0), A(2, 0), A(0, 1), A(1, 1), A(2, 1), A(0, 2), + A(1, 2), A(2, 2)); +} + +/// Copy of the i-th row of A. +constexpr Vector row(const Matrix& A, std::size_t i) noexcept { + return Vector(A(i, 0), A(i, 1), A(i, 2)); +} + +/// Copy of the j-th column of A. +constexpr Vector column(const Matrix& A, std::size_t j) noexcept { + return Vector(A(0, j), A(1, j), A(2, j)); } -inline bool operator!=(const Matrix& A, const Matrix& B) { return !(A == B); } +constexpr void setrow(Matrix& A, std::size_t i, const Vector& v) noexcept { + A(i, 0) = v[0]; + A(i, 1) = v[1]; + A(i, 2) = v[2]; +} + +constexpr void setcolumn(Matrix& A, std::size_t j, const Vector& v) noexcept { + A(0, j) = v[0]; + A(1, j) = v[1]; + A(2, j) = v[2]; +} // Hashing ------------------------------------------------------------------- DLL_EXPORT size_t hash_value(const Vector& v); DLL_EXPORT size_t hash_value(const Matrix& A); -// Inlined functions --------------------------------------------------------- +// Vector functions ---------------------------------------------------------- -inline const Vector& floor(const Vector& v) { - static Vector res; - Vector::const_iterator xi = v.begin(); - Vector::iterator xo = res.begin(); - for (; xi != v.end(); ++xi, ++xo) *xo = std::floor(*xi); - return res; +inline Vector floor(const Vector& v) { + return Vector(std::floor(v[0]), std::floor(v[1]), std::floor(v[2])); } -// Template functions -------------------------------------------------------- - template -double norm(const V& u) { - return sqrt(R3::dot(u, u)); +double dot(const V& u, const V& v) { + return (u[0] * v[0] + u[1] * v[1] + u[2] * v[2]); } template -double distance(const V& u, const V& v) { - static R3::Vector duv; - duv[0] = u[0] - v[0]; - duv[1] = u[1] - v[1]; - duv[2] = u[2] - v[2]; - return R3::norm(duv); +double norm(const V& u) { + return std::sqrt(R3::dot(u, u)); } template -double dot(const V& u, const V& v) { - return (u[0] * v[0] + u[1] * v[1] + u[2] * v[2]); +double distance(const V& u, const V& v) { + const double dx = u[0] - v[0]; + const double dy = u[1] - v[1]; + const double dz = u[2] - v[2]; + return std::sqrt(dx * dx + dy * dy + dz * dz); } template Vector cross(const V& u, const V& v) { - Vector res; - res[0] = u[1] * v[2] - u[2] * v[1]; - res[1] = u[2] * v[0] - u[0] * v[2]; - res[2] = u[0] * v[1] - u[1] * v[0]; - return res; + return Vector(u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], + u[0] * v[1] - u[1] * v[0]); } +/// Matrix-vector product M * u. template -const Vector& mxvecproduct(const Matrix& M, const V& u) { - static Vector res; - res[0] = M(0, 0) * u[0] + M(0, 1) * u[1] + M(0, 2) * u[2]; - res[1] = M(1, 0) * u[0] + M(1, 1) * u[1] + M(1, 2) * u[2]; - res[2] = M(2, 0) * u[0] + M(2, 1) * u[1] + M(2, 2) * u[2]; - return res; +Vector mxvecproduct(const Matrix& M, const V& u) { + return Vector(M(0, 0) * u[0] + M(0, 1) * u[1] + M(0, 2) * u[2], + M(1, 0) * u[0] + M(1, 1) * u[1] + M(1, 2) * u[2], + M(2, 0) * u[0] + M(2, 1) * u[1] + M(2, 2) * u[2]); } +/// Vector-matrix product u * M. template -const Vector& mxvecproduct(const V& u, const Matrix& M) { - static Vector res; - res[0] = u[0] * M(0, 0) + u[1] * M(1, 0) + u[2] * M(2, 0); - res[1] = u[0] * M(0, 1) + u[1] * M(1, 1) + u[2] * M(2, 1); - res[2] = u[0] * M(0, 2) + u[1] * M(1, 2) + u[2] * M(2, 2); - return res; +Vector mxvecproduct(const V& u, const Matrix& M) { + return Vector(u[0] * M(0, 0) + u[1] * M(1, 0) + u[2] * M(2, 0), + u[0] * M(0, 1) + u[1] * M(1, 1) + u[2] * M(2, 1), + u[0] * M(0, 2) + u[1] * M(1, 2) + u[2] * M(2, 2)); } } // namespace R3 @@ -261,4 +329,7 @@ EpsilonEqual::operator()( } // namespace mathutils } // namespace diffpy +// Implementation of serialization +#include + #endif // R3LINALG_HPP_INCLUDED diff --git a/src/diffpy/srreal/R3linalg.ipp b/src/diffpy/srreal/R3linalg.ipp new file mode 100644 index 0000000..068000b --- /dev/null +++ b/src/diffpy/srreal/R3linalg.ipp @@ -0,0 +1,102 @@ +/***************************************************************************** + * + * libdiffpy by DANSE Diffraction group + * Simon J. L. Billinge + * (c) 2009 The Trustees of Columbia University + * in the City of New York. All rights reserved. + * + * File coded by: Pavol Juhas + * + * See AUTHORS.rst for a list of people who contributed. + * See LICENSE_DANSE.rst for license information. + * + ****************************************************************************** + * + * R3linalg -- serialization of the Vector and Matrix types. + * + *****************************************************************************/ + +#ifndef R3LINALG_IPP_INCLUDED +#define R3LINALG_IPP_INCLUDED + +#include +#include +#include +#include + +namespace diffpy { +namespace srreal { +namespace R3 { +namespace detail { + +// These wrappers keep the archive layout of the former ublas base types, so +// older archives and pickles still load. Each stands for one ublas type and +// must stay a distinct class with default serialization traits, because the +// archive writes class information once per type. + +inline void check_legacy_size(bool valid) { + if (!valid) { + throw boost::archive::archive_exception( + boost::archive::archive_exception::array_size_too_short); + } +} + +/// Stands for ublas::bounded_array. +template +struct LegacyBoundedArray { + std::array& a; + + template + void serialize(Archive& ar, const unsigned int) { + boost::serialization::collection_size_type s(N); + ar& boost::serialization::make_nvp("size", s); + check_legacy_size(s == N); + ar& boost::serialization::make_array(a.data(), N); + } +}; + +/// Stands for ublas::vector >. +struct LegacyVectorBase { + std::array& a; + + template + void serialize(Archive& ar, const unsigned int) { + LegacyBoundedArray<3> data{a}; + ar& boost::serialization::make_nvp("data", data); + } +}; + +/// Stands for ublas::matrix >. +struct LegacyMatrixBase { + std::array& a; + + template + void serialize(Archive& ar, const unsigned int) { + boost::serialization::collection_size_type s1(3), s2(3); + ar& boost::serialization::make_nvp("size1", s1); + ar& boost::serialization::make_nvp("size2", s2); + check_legacy_size(s1 == 3 && s2 == 3); + LegacyBoundedArray<9> data{a}; + ar& boost::serialization::make_nvp("data", data); + } +}; + +} // namespace detail + +template +void Vector::serialize(Archive& ar, const unsigned int) { + detail::LegacyVectorBase base{mdata}; + ar & base; +} + +template +void Matrix::serialize(Archive& ar, const unsigned int) { + detail::LegacyMatrixBase base{mdata}; + ar & base; +} + +} // namespace R3 +} // namespace srreal +} // namespace diffpy + +#endif // R3LINALG_IPP_INCLUDED diff --git a/src/tests/TestHash.cpp b/src/tests/TestHash.cpp new file mode 100644 index 0000000..89594c5 --- /dev/null +++ b/src/tests/TestHash.cpp @@ -0,0 +1,68 @@ +/***************************************************************************** + * + * 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. + * + ****************************************************************************** + * + * TestHash -- unit tests for the hash utilities in diffpy/hash.hpp + * + *****************************************************************************/ + +#include +#include + +#include + +#include +#include +#include + +using namespace std; +using namespace diffpy::srreal; + +TEST_CASE("TestHash::test_vector", "[hash]") { + R3::Vector v0(1.0, 2.0, 3.0); + R3::Vector v1(1.0, 2.0, 3.0); + CHECK(hash_value(v0) == hash_value(v1)); + v1[2] = 4.0; + CHECK(hash_value(v0) != hash_value(v1)); + // negative zero compares equal, so it must hash equal + R3::Vector z0(0.0, 0.0, 0.0); + R3::Vector z1(-0.0, 0.0, -0.0); + REQUIRE(z0 == z1); + CHECK(hash_value(z0) == hash_value(z1)); +} + +TEST_CASE("TestHash::test_atom", "[hash]") { + Atom a0; + a0.atomtype = "C"; + a0.xyz_cartn = R3::Vector(1, 2, 3); + Atom a1 = a0; + CHECK(hash_value(a0) == hash_value(a1)); + a1.uij_cartn(0, 0) = 0.01; + CHECK(hash_value(a0) != hash_value(a1)); + // hash_of finds the ADL hash_value overload + CHECK(diffpy::hash_of(a0) == hash_value(a0)); +} + +TEST_CASE("TestHash::test_hash_values", "[hash]") { + const R3::Vector v(1.0, 2.0, 3.0); + size_t seed = 0; + diffpy::hash_combine(seed, string("C")); + diffpy::hash_combine(seed, v); + diffpy::hash_combine(seed, 7); + CHECK(diffpy::hash_values(string("C"), v, 7) == seed); + CHECK(diffpy::hash_values(1, 2) != diffpy::hash_values(2, 1)); + CHECK(diffpy::hash_values() == 0); +} + +TEST_CASE("TestHash::test_pair_hash", "[hash]") { + diffpy::PairHash h; + CHECK(h(make_pair(1, 2)) != h(make_pair(2, 1))); +} diff --git a/src/tests/TestLattice.cpp b/src/tests/TestLattice.cpp index b304031..55dc115 100644 --- a/src/tests/TestLattice.cpp +++ b/src/tests/TestLattice.cpp @@ -16,6 +16,8 @@ * *****************************************************************************/ +#include + #include #include @@ -117,6 +119,37 @@ TEST_CASE_METHOD(TestLatticeFixture, "TestLattice::test_setLatBase", CHECK(allclose(recbase_check, lattice->recbase())); } +TEST_CASE_METHOD(TestLatticeFixture, "TestLattice::test_invalid_unchanged", + "[Lattice]") { + // A degenerate lattice throws and leaves the lattice as it was. + lattice->setLatPar(1.0, 2.0, 3.0, 80, 90, 100); + const Lattice saved = *lattice; + const R3::Vector v(0.3, -1.7, 2.9); + auto unchanged = [&]() { + return lattice->a() == saved.a() && lattice->b() == saved.b() && + lattice->c() == saved.c() && lattice->alpha() == saved.alpha() && + lattice->beta() == saved.beta() && + lattice->gamma() == saved.gamma() && lattice->ar() == saved.ar() && + lattice->base() == saved.base() && + lattice->recbase() == saved.recbase() && + lattice->baserot() == saved.baserot() && + lattice->metrics() == saved.metrics() && + lattice->fractional(v) == saved.fractional(v) && + lattice->cartesian(v) == saved.cartesian(v); + }; + REQUIRE(unchanged()); + CHECK_THROWS_AS(lattice->setLatPar(0.0, 4.0, 4.0, 90, 90, 90), + std::invalid_argument); + CHECK(unchanged()); + const R3::Vector va(1.0, 0.0, 0.0); + const R3::Vector vc(0.0, 0.0, 1.0); + CHECK_THROWS_AS(lattice->setLatBase(va, va, vc), std::invalid_argument); + CHECK(unchanged()); + CHECK_THROWS_AS(lattice->setLatBase(va, vc, R3::Vector(0.0, 1.0, 0.0)), + std::invalid_argument); + CHECK(unchanged()); +} + TEST_CASE_METHOD(TestLatticeFixture, "TestLattice::test_dist", "[Lattice]") { R3::Vector va(1.0, 2.0, 2.0); R3::Vector vb(0.0, 0.0, 0.0); diff --git a/src/tests/TestR3linalg.cpp b/src/tests/TestR3linalg.cpp index 53307e9..87219d7 100644 --- a/src/tests/TestR3linalg.cpp +++ b/src/tests/TestR3linalg.cpp @@ -16,12 +16,22 @@ * *****************************************************************************/ +#include +#include +#include +#include +#include +#include + +#include +#include #include #include #include using Catch::Matchers::WithinAbs; +using Catch::Matchers::WithinRel; using namespace std; using namespace diffpy::srreal; @@ -47,6 +57,11 @@ TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_determinant", double detA = 549; CHECK(detA == R3::determinant(A1)); CHECK(-detA == R3::determinant(A2)); + // singular matrices with dependent rows and with a zero column + R3::Matrix S1(1, 2, 3, 2, 4, 6, 1, 0, 1); + R3::Matrix S2(0, 1, 2, 0, 3, 4, 0, 5, 6); + CHECK(0.0 == R3::determinant(S1)); + CHECK(0.0 == R3::determinant(S2)); } TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_inverse", @@ -57,6 +72,28 @@ TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_inverse", -0.24089965988852, 1.32489393619466, -0.15249839300481, 0.68609361943181, 0.02249568627913, 0.41178393851247); CHECK(allclose(invA, R3::inverse(A))); + // the argument may be the matrix returned from the previous call + CHECK(allclose(A, R3::inverse(R3::inverse(A)))); + R3::Matrix AinvA = R3::prod(A, R3::inverse(A)); + CHECK(allclose(R3::identity(), AinvA)); + // singular matrices with dependent rows and with a zero column + R3::Matrix S1(1, 2, 3, 2, 4, 6, 1, 0, 1); + R3::Matrix S2(0, 1, 2, 0, 3, 4, 0, 5, 6); + CHECK_THROWS_AS(R3::inverse(S1), invalid_argument); + CHECK_THROWS_AS(R3::inverse(S2), invalid_argument); +} + +TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_nearly_singular", + "[R3linalg]") { + // Determinant d * d is lost in round-off when evaluated from cofactors. + const double x = 1 + 1e-8; + const double d = x - 1; + R3::Matrix A(1, 1, 1, 1, x, 1, 1, 1, x); + CHECK_THAT(R3::determinant(A), WithinRel(d * d, 1e-12)); + // exact inverse scaled by d + R3::Matrix dinvA0(2 + d, -1, -1, -1, 1, 0, -1, 0, 1); + R3::Matrix dinvA = d * R3::inverse(A); + CHECK(allclose(dinvA0, dinvA)); } TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_transpose", @@ -102,4 +139,165 @@ TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_product", CHECK(allclose(prod_vM, R3::mxvecproduct(v, M))); } +TEST_CASE("TestR3linalg::test_default_zero", "[R3linalg]") { + R3::Vector v; + R3::Matrix M; + CHECK(R3::Vector(0, 0, 0) == v); + CHECK(R3::zeromatrix() == M); + CHECK(R3::Vector(0, 0, 0) == R3::zerovector); +} + +TEST_CASE("TestR3linalg::test_arithmetic", "[R3linalg]") { + R3::Vector u(1, 2, 3); + R3::Vector v(4, 5, 7); + CHECK(R3::Vector(5, 7, 10) == u + v); + CHECK(R3::Vector(-3, -3, -4) == u - v); + CHECK(R3::Vector(-1, -2, -3) == -u); + CHECK(R3::Vector(2, 4, 6) == u * 2.0); + CHECK(R3::Vector(2, 4, 6) == 2.0 * u); + CHECK(R3::Vector(0.5, 1, 1.5) == u / 2.0); + u += v; + CHECK(R3::Vector(5, 7, 10) == u); + u -= v; + CHECK(R3::Vector(1, 2, 3) == u); + R3::Matrix A(1, 2, 3, 4, 5, 6, 7, 8, 9); + CHECK(R3::Matrix(2, 4, 6, 8, 10, 12, 14, 16, 18) == A + A); + CHECK(R3::zeromatrix() == A - A); + CHECK(A * 2.0 == 2.0 * A); + CHECK(A == (A * 4.0) / 4.0); + CHECK(R3::zeromatrix() == A + (-A)); +} + +TEST_CASE("TestR3linalg::test_rows_columns", "[R3linalg]") { + R3::Matrix A(1, 2, 3, 4, 5, 6, 7, 8, 9); + CHECK(R3::Vector(4, 5, 6) == R3::row(A, 1)); + CHECK(R3::Vector(3, 6, 9) == R3::column(A, 2)); + R3::setrow(A, 0, R3::Vector(-1, -2, -3)); + CHECK(R3::Matrix(-1, -2, -3, 4, 5, 6, 7, 8, 9) == A); + R3::setcolumn(A, 1, R3::Vector(0, 0, 0)); + CHECK(R3::Matrix(-1, 0, -3, 4, 0, 6, 7, 0, 9) == A); +} + +namespace { + +// T&& is an rvalue when T is not a reference type +template +concept ScalableInPlace = requires(T&& x) { std::forward(x) *= 2.0; }; + +template +concept VectorElementWritable = requires(T&& x) { + std::forward(x)[1] = 2.0; + std::forward(x)(1) += 2.0; + *std::forward(x).begin() = 2.0; + std::forward(x).data()[1] = 2.0; +}; + +template +concept MatrixElementWritable = requires(T&& x) { + std::forward(x)(0, 1) = 2.0; + std::forward(x).data()[1] = 2.0; +}; + +} // namespace + +// Writes to temporaries such as row(A, 0) = v or row(A, 0)[1] = x must not +// compile, because they would silently leave A unchanged. +static_assert(std::is_assignable_v); +static_assert(!std::is_assignable_v); +static_assert(ScalableInPlace); +static_assert(!ScalableInPlace); +static_assert(VectorElementWritable); +static_assert(!VectorElementWritable); +static_assert(std::is_assignable_v); +static_assert(!std::is_assignable_v); +static_assert(ScalableInPlace); +static_assert(!ScalableInPlace); +static_assert(MatrixElementWritable); +static_assert(!MatrixElementWritable); +static_assert(std::is_trivially_copyable_v); +static_assert(std::is_trivially_copyable_v); + +TEST_CASE("TestR3linalg::test_prod", "[R3linalg]") { + R3::Matrix A(1, 2, 3, 4, 5, 6, 7, 8, 9); + R3::Matrix B(9, 8, 7, 6, 5, 4, 3, 2, 1); + R3::Matrix AB(30, 24, 18, 84, 69, 54, 138, 114, 90); + CHECK(AB == R3::prod(A, B)); + CHECK(A == R3::prod(A, R3::identity())); +} + +TEST_CASE_METHOD(TestR3linalgFixture, "TestR3linalg::test_inverse_pivoting", + "[R3linalg]") { + // zero at the first pivot requires row exchanges + R3::Matrix P(0, 1, 0, 0, 0, 1, 1, 0, 0); + CHECK(1.0 == R3::determinant(P)); + CHECK(R3::trans(P) == R3::inverse(P)); + R3::Matrix A(0, 2, 1, 3, 0, 4, 5, 6, 0); + CHECK_THAT(R3::determinant(A), WithinAbs(58.0, precision)); + CHECK(allclose(R3::identity(), R3::prod(A, R3::inverse(A)))); +} + +namespace { + +// Binary archive of a Vector(1, -2.5, 3.25), Matrix(1, 2, ..., 8, 9.5) and +// the same two objects again, as written by the boost::numeric::ublas based +// implementation. The archive uses no_header, little-endian byte order and +// 8-byte collection sizes. +const char* const LEGACY_ARCHIVE_HEX = + "0000000000000000000000000000000300000000000000000000000000f03f00" + "000000000004c00000000000000a400000000000000000000003000000000000" + "00030000000000000000000000000900000000000000000000000000f03f0000" + "0000000000400000000000000840000000000000104000000000000014400000" + "0000000018400000000000001c40000000000000204000000000000023400300" + "000000000000000000000000f03f00000000000004c00000000000000a400300" + "00000000000003000000000000000900000000000000000000000000f03f0000" + "0000000000400000000000000840000000000000104000000000000014400000" + "0000000018400000000000001c4000000000000020400000000000002340"; + +std::string fromhex(const std::string& hex) { + std::string rv; + for (size_t i = 0; i < hex.size(); i += 2) { + rv += char(std::stoi(hex.substr(i, 2), nullptr, 16)); + } + return rv; +} + +} // namespace + +TEST_CASE("TestR3linalg::test_legacy_serialization", "[R3linalg]") { + if (std::endian::native != std::endian::little || sizeof(size_t) != 8) { + SKIP("legacy archive bytes assume little-endian 64-bit platform"); + } + const R3::Vector v(1.0, -2.5, 3.25); + const R3::Matrix M(1, 2, 3, 4, 5, 6, 7, 8, 9.5); + using boost::archive::no_header; + const std::string legacy = fromhex(LEGACY_ARCHIVE_HEX); + // saving produces the same bytes as the ublas based implementation + std::ostringstream os; + { + boost::archive::binary_oarchive oa(os, no_header); + oa << v << M << v << M; + } + CHECK(legacy == os.str()); + // loading of the legacy archive + std::istringstream is(legacy); + { + boost::archive::binary_iarchive ia(is, no_header); + R3::Vector v1, v2; + R3::Matrix M1, M2; + ia >> v1 >> M1 >> v2 >> M2; + CHECK(v == v1); + CHECK(M == M1); + CHECK(v == v2); + CHECK(M == M2); + } + // archive with a wrong Vector size cannot be loaded + std::string bad = legacy; + REQUIRE(3 == bad[15]); + bad[15] = 4; + std::istringstream isbad(bad); + boost::archive::binary_iarchive iabad(isbad, no_header); + R3::Vector vbad; + CHECK_THROWS_AS(iabad >> vbad, boost::archive::archive_exception); +} + // End of file diff --git a/src/tests/test_helpers.cpp b/src/tests/test_helpers.cpp index c11da3b..33c37ba 100644 --- a/src/tests/test_helpers.cpp +++ b/src/tests/test_helpers.cpp @@ -19,15 +19,27 @@ #include #include #include -#include -#include -#include +#include +#include +#include #include #include #include "test_helpers.hpp" +namespace { + +char lower_char(char c) { + return static_cast(std::tolower(static_cast(c))); +} + +bool iequals(std::string_view a, std::string_view b) { + return std::ranges::equal(a, b, {}, lower_char, lower_char); +} + +} // namespace + // Resolve DIFFPYTESTSDIRPATH ------------------------------------------------ #ifndef DIFFPYTESTSDIRPATH @@ -71,12 +83,12 @@ diffpy::srreal::StructureAdapterPtr loadTestPeriodicStructure( string fileformat; // process structure file header while (fp >> line) { - if (line.wcount() >= 2 && boost::iequals(line.words[0], "format")) { + if (line.wcount() >= 2 && iequals(line.words[0], "format")) { fileformat = line.words[1]; - assert(boost::iequals(fileformat, "pdffit")); + assert(iequals(fileformat, "pdffit")); } if (line.wcount() == 7 && line.words[0] == "cell") { - boost::replace_all(line.line, ",", " "); + std::ranges::replace(line.line, ',', ' '); string skip; double a, b, c, alpha, beta, gamma; istringstream fpwords(line.line); @@ -101,7 +113,7 @@ diffpy::srreal::StructureAdapterPtr loadTestPeriodicStructure( Uc(2, 1) = Uc(1, 2); if (!fp) break; // convert string symbol to standard case - boost::to_lower(ai.atomtype); + std::ranges::transform(ai.atomtype, ai.atomtype.begin(), lower_char); string::iterator tt = ai.atomtype.begin(); for (; tt != ai.atomtype.end(); ++tt) { if (isalpha(*tt)) { From 28ea6aee167d4f0421939d56f33d8941864fe11b Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Tue, 10 Nov 2026 21:48:40 -0800 Subject: [PATCH 2/6] perf: vectorize Gaussian and Debye calculations --- .github/workflows/matrix-and-codecov.yml | 6 +- CMakeLists.txt | 2 + cmake/CalculationKernels.cmake | 107 ++++ docs/source/development.rst | 16 + news/faster-pdf.rst | 15 +- requirements/build.txt | 1 + src/diffpy/srreal/BaseDebyeSum.cpp | 42 +- src/diffpy/srreal/BaseDebyeSum.hpp | 3 +- src/diffpy/srreal/CalculationKernel.hpp | 43 ++ src/diffpy/srreal/CroppedGaussianProfile.cpp | 8 + src/diffpy/srreal/CroppedGaussianProfile.hpp | 4 + src/diffpy/srreal/GaussianProfile.cpp | 34 ++ src/diffpy/srreal/GaussianProfile.hpp | 12 + src/diffpy/srreal/PDFCalculator.cpp | 11 +- src/diffpy/srreal/PeakProfile.cpp | 14 + src/diffpy/srreal/PeakProfile.hpp | 6 + .../srreal/kernels/CalculationKernelAVX2.cpp | 32 ++ .../kernels/CalculationKernelAVX512.cpp | 32 ++ .../kernels/CalculationKernelDispatch.cpp | 113 ++++ .../srreal/kernels/CalculationKernelImpl.hpp | 111 ++++ .../srreal/kernels/CalculationKernelNEON.cpp | 32 ++ .../srreal/kernels/CalculationKernelSIMD.ipp | 283 ++++++++++ .../srreal/kernels/CalculationKernelSSE2.cpp | 32 ++ .../kernels/CalculationKernelScalar.cpp | 75 +++ src/tests/CMakeLists.txt | 6 + src/tests/TestCalculationKernel.cpp | 487 ++++++++++++++++++ src/tests/test_helpers.cpp | 20 + src/tests/test_helpers.hpp | 36 ++ 28 files changed, 1547 insertions(+), 36 deletions(-) create mode 100644 cmake/CalculationKernels.cmake create mode 100644 src/diffpy/srreal/CalculationKernel.hpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelAVX2.cpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelAVX512.cpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelDispatch.cpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelImpl.hpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelNEON.cpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelSSE2.cpp create mode 100644 src/diffpy/srreal/kernels/CalculationKernelScalar.cpp create mode 100644 src/tests/TestCalculationKernel.cpp diff --git a/.github/workflows/matrix-and-codecov.yml b/.github/workflows/matrix-and-codecov.yml index 8fcbcdd..a4a6d37 100644 --- a/.github/workflows/matrix-and-codecov.yml +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1a7f3fa..0564a61 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -320,6 +320,8 @@ if(LIBDIFFPY_HAS_OBJCRYST) target_link_libraries(diffpy PUBLIC ObjCryst::ObjCryst) endif() +include(CalculationKernels) + # Installation ---------------------------------------------------------------- if(LIBDIFFPY_INSTALL) diff --git a/cmake/CalculationKernels.cmake b/cmake/CalculationKernels.cmake new file mode 100644 index 0000000..f473993 --- /dev/null +++ b/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 + "$<$>:-O2>") +endif() + +# Add the SIMD kernel in CalculationKernel.cpp, compiled with the +# instruction set flags in ARGN. CalculationKernelDispatch.cpp uses the +# kernel when LIBDIFFPY_KERNEL_ 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 + $) +endif() diff --git a/docs/source/development.rst b/docs/source/development.rst index 49ecc45..41db078 100644 --- a/docs/source/development.rst +++ b/docs/source/development.rst @@ -47,6 +47,22 @@ directory with ``cmake -S . -B -G ""``. To disable the ObjCryst adapter, pass ``-DLIBDIFFPY_ENABLE_OBJCRYST=OFF``. For gprof instrumentation, use ``-DLIBDIFFPY_PROFILE=ON``. +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 `_ 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 ---------------------------------------- diff --git a/news/faster-pdf.rst b/news/faster-pdf.rst index c9f32d1..10b98bf 100644 --- a/news/faster-pdf.rst +++ b/news/faster-pdf.rst @@ -2,6 +2,15 @@ * ``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:** @@ -16,10 +25,11 @@ 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. **Removed:** -* Dependency on Boost.Functional/Hash, Boost.StringAlgo and Boost.uBLAS. +* Dependencies on Boost.Functional/Hash, Boost.StringAlgo and Boost.uBLAS. ``R3::Vector`` and ``R3::Matrix`` no longer accept uBLAS expressions. **Fixed:** @@ -28,3 +38,6 @@ 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. diff --git a/requirements/build.txt b/requirements/build.txt index b3fb0da..abe058a 100644 --- a/requirements/build.txt +++ b/requirements/build.txt @@ -1,3 +1,4 @@ cxx-compiler cmake>=3.22 ninja +xsimd diff --git a/src/diffpy/srreal/BaseDebyeSum.cpp b/src/diffpy/srreal/BaseDebyeSum.cpp index 9572f43..ea32084 100644 --- a/src/diffpy/srreal/BaseDebyeSum.cpp +++ b/src/diffpy/srreal/BaseDebyeSum.cpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -141,19 +142,27 @@ void BaseDebyeSum::addPairContribution(const BaseBondGenerator& bnds, // calculate sigma parameter for the Debye-Waller dampign Gaussian const double fwhm = this->getPeakWidthModel()->calculate(bnds); const double fwhmtosigma = 1.0 / (2 * sqrt(2 * M_LN2)); - const double dwsigma = fwhmtosigma * fwhm; const int nqpts = pdfutils_qmaxSteps(this); - const int smscale = summationscale * bnds.multiplicity(); - const double& sineprec = this->getDebyePrecision(); - for (int kq = pdfutils_qminSteps(this); kq < nqpts; ++kq) { - const double q = kq * this->getQstep(); - const double dwscale = exp(-0.5 * pow(dwsigma * q, 2)); - const double sinescale = smscale * dwscale * - this->sfSiteAtkQ(bnds.site0(), kq) * - this->sfSiteAtkQ(bnds.site1(), kq) / dist; - if (eps_eq(0.0, sinescale, sineprec)) break; - mvalue[kq] += sinescale * sin(q * dist); - } + const auto& sftypeatkq = mstructure_cache.sftypeatkq; + const auto& typeofsite = mstructure_cache.typeofsite; + const auto& sf0 = sftypeatkq[typeofsite[bnds.site0()]]; + const auto& sf1 = sftypeatkq[typeofsite[bnds.site1()]]; + assert(nqpts <= int(mvalue.size())); + assert(nqpts <= int(sf0.size())); + assert(nqpts <= int(sf1.size())); + + kernels::activeCalculationKernel().debyeSum({ + .out = mvalue.data(), + .kqlo = pdfutils_qminSteps(this), + .kqhi = nqpts, + .qstep = this->getQstep(), + .dist = dist, + .dwsigma = fwhmtosigma * fwhm, + .scale = double(summationscale * bnds.multiplicity()), + .prec = this->getDebyePrecision(), + .sf0 = sf0.data(), + .sf1 = sf1.data(), + }); } void BaseDebyeSum::stashPartialValue() { mdbsumstash = this->value(); } @@ -170,15 +179,6 @@ double BaseDebyeSum::sfSiteAtQ(int siteidx, const double& q) const { // Private Methods ----------------------------------------------------------- -double BaseDebyeSum::sfSiteAtkQ(int siteidx, int kq) const { - assert(0 <= siteidx && siteidx < int(mstructure_cache.typeofsite.size())); - int typeidx = mstructure_cache.typeofsite[siteidx]; - assert(typeidx < int(mstructure_cache.sftypeatkq.size())); - const QuantityType& sfarray = mstructure_cache.sftypeatkq[typeidx]; - assert(0 <= kq && kq < int(sfarray.size())); - return sfarray[kq]; -} - double BaseDebyeSum::sfAverageAtkQ(int kq) const { assert(kq < int(mstructure_cache.sfaverageatkq.size())); return mstructure_cache.sfaverageatkq[kq]; diff --git a/src/diffpy/srreal/BaseDebyeSum.hpp b/src/diffpy/srreal/BaseDebyeSum.hpp index d290168..b5cd11a 100644 --- a/src/diffpy/srreal/BaseDebyeSum.hpp +++ b/src/diffpy/srreal/BaseDebyeSum.hpp @@ -71,9 +71,8 @@ class DLL_EXPORT BaseDebyeSum : public PairQuantity, private: // methods - /// cache structure factors data for a quick access during summation - double sfSiteAtkQ(int siteidx, int kq) const; double sfAverageAtkQ(int kq) const; + /// cache structure factors data for a quick access during summation void cacheStructureData(); // data diff --git a/src/diffpy/srreal/CalculationKernel.hpp b/src/diffpy/srreal/CalculationKernel.hpp new file mode 100644 index 0000000..35902a1 --- /dev/null +++ b/src/diffpy/srreal/CalculationKernel.hpp @@ -0,0 +1,43 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Runtime selection of the PDF calculation kernel + * + *****************************************************************************/ + +#ifndef CALCULATIONKERNEL_HPP_INCLUDED +#define CALCULATIONKERNEL_HPP_INCLUDED + +#include +#include + +#include + +namespace diffpy { +namespace srreal { + +/// Names of the kernels that this build and CPU can run, from the most to +/// the least preferred: "avx512", "avx2" and "sse2" on x86-64, "neon" on +/// arm64, and "scalar", which is last and always available. +DLL_EXPORT std::vector availableCalculationKernels(); + +/// Name of the kernel in use, by default the most preferred one. +DLL_EXPORT std::string getCalculationKernel(); + +/// Use the named kernel for all calculators. +/// Throw std::invalid_argument when the kernel is not available. +/// Do not call while a calculation runs in another thread. +DLL_EXPORT void setCalculationKernel(const std::string& name); + +} // namespace srreal +} // namespace diffpy + +#endif // CALCULATIONKERNEL_HPP_INCLUDED diff --git a/src/diffpy/srreal/CroppedGaussianProfile.cpp b/src/diffpy/srreal/CroppedGaussianProfile.cpp index 84d142a..b2990fe 100644 --- a/src/diffpy/srreal/CroppedGaussianProfile.cpp +++ b/src/diffpy/srreal/CroppedGaussianProfile.cpp @@ -19,6 +19,7 @@ *****************************************************************************/ #include +#include #include @@ -73,6 +74,13 @@ void CroppedGaussianProfile::setPrecision(double eps) { } } +// Protected Methods --------------------------------------------------------- + +GaussianProfile::KernelShape CroppedGaussianProfile::kernelShape() const { + const bool exact = (typeid(*this) == typeid(CroppedGaussianProfile)); + return exact ? KernelShape::CroppedGaussian : KernelShape::None; +} + // Registration -------------------------------------------------------------- bool reg_CroppedGaussianProfile = CroppedGaussianProfile().registerThisType(); diff --git a/src/diffpy/srreal/CroppedGaussianProfile.hpp b/src/diffpy/srreal/CroppedGaussianProfile.hpp index b393565..88fcdb5 100644 --- a/src/diffpy/srreal/CroppedGaussianProfile.hpp +++ b/src/diffpy/srreal/CroppedGaussianProfile.hpp @@ -40,6 +40,10 @@ class DLL_EXPORT CroppedGaussianProfile : public GaussianProfile { double operator()(double x, double fwhm) const; void setPrecision(double eps); + protected: + // methods + KernelShape kernelShape() const override; + private: // data double mscale; diff --git a/src/diffpy/srreal/GaussianProfile.cpp b/src/diffpy/srreal/GaussianProfile.cpp index 1773c2b..d8fa02a 100644 --- a/src/diffpy/srreal/GaussianProfile.cpp +++ b/src/diffpy/srreal/GaussianProfile.cpp @@ -18,8 +18,10 @@ *****************************************************************************/ #include +#include #include +#include #include #include @@ -70,6 +72,31 @@ double GaussianProfile::xboundhi(double fwhm) const { return rv; } +void GaussianProfile::addRDFPeak(double* out, int ilo, int ihi, int rlosteps, + double rstep, double dist, double fwhm, + double scale) const { + const KernelShape shape = this->kernelShape(); + if (shape == KernelShape::None) { + this->PeakProfile::addRDFPeak(out, ilo, ihi, rlosteps, rstep, dist, fwhm, + scale); + return; + } + if (fwhm <= 0.0) return; + kernels::activeCalculationKernel().gaussianWindow({ + .out = out, + .ilo = ilo, + .ihi = ihi, + .rlosteps = rlosteps, + .rstep = rstep, + .dist = dist, + .fwhm = fwhm, + .scale = scale, + .peakheight = (*this)(0.0, fwhm), + .cropped = (shape == KernelShape::CroppedGaussian), + .halfboundrel = mhalfboundrel, + }); +} + void GaussianProfile::setPrecision(double eps) { // correct any settings below DOUBLE_EPS double eps1 = max(eps, DOUBLE_EPS); @@ -80,6 +107,13 @@ void GaussianProfile::setPrecision(double eps) { mhalfboundrel = 0.0; } +// Protected Methods --------------------------------------------------------- + +GaussianProfile::KernelShape GaussianProfile::kernelShape() const { + const bool exact = (typeid(*this) == typeid(GaussianProfile)); + return exact ? KernelShape::Gaussian : KernelShape::None; +} + // Registration -------------------------------------------------------------- bool reg_GaussianProfile = GaussianProfile().registerThisType(); diff --git a/src/diffpy/srreal/GaussianProfile.hpp b/src/diffpy/srreal/GaussianProfile.hpp index 84ef7d3..c0cd598 100644 --- a/src/diffpy/srreal/GaussianProfile.hpp +++ b/src/diffpy/srreal/GaussianProfile.hpp @@ -39,12 +39,24 @@ class DLL_EXPORT GaussianProfile : public PeakProfile { double operator()(double x, double fwhm) const; double xboundlo(double fwhm) const; double xboundhi(double fwhm) const; + void addRDFPeak(double* out, int ilo, int ihi, int rlosteps, double rstep, + double dist, double fwhm, double scale) const override; void setPrecision(double eps); protected: // data double mhalfboundrel; + /// Profiles that addRDFPeak can evaluate with the calculation kernel. + /// CroppedGaussian is zero for |x| >= mhalfboundrel * fwhm. + enum class KernelShape { None, Gaussian, CroppedGaussian }; + + // methods + /// Shape of operator() for the kernel in addRDFPeak, None for the generic + /// loop. It is None for derived classes, which may change operator(), + /// unless they override this. + virtual KernelShape kernelShape() const; + private: // serialization friend class boost::serialization::access; diff --git a/src/diffpy/srreal/PDFCalculator.cpp b/src/diffpy/srreal/PDFCalculator.cpp index 144f553..d464349 100644 --- a/src/diffpy/srreal/PDFCalculator.cpp +++ b/src/diffpy/srreal/PDFCalculator.cpp @@ -433,15 +433,8 @@ void PDFCalculator::addPairContribution(const BaseBondGenerator& bnds, int ilast = min(this->countCalcPoints(), this->calcIndex(xhi) + 1); assert(ilast <= int(mvalue.size())); assert(eps_gt(dist, 0.0)); - for (; i < ilast; ++i) { - double x = (this->rcalcloSteps() + i) * this->getRstep() - dist; - double y = pkf(x, fwhm); - // Contributions in G(r) need to be normalized by pair distance, - // not by r as done in PDFfit or PDFfit2. Here we rescale RDF - // in such way that division by r will give a correct result. - double yrdf = y * (x / dist + 1); - mvalue[i] += peakscale * yrdf; - } + pkf.addRDFPeak(mvalue.data(), i, ilast, this->rcalcloSteps(), + this->getRstep(), dist, fwhm, peakscale); } void PDFCalculator::stashPartialValue() { diff --git a/src/diffpy/srreal/PeakProfile.cpp b/src/diffpy/srreal/PeakProfile.cpp index 27a8db1..2cd72c5 100644 --- a/src/diffpy/srreal/PeakProfile.cpp +++ b/src/diffpy/srreal/PeakProfile.cpp @@ -41,6 +41,20 @@ PeakProfile::PeakProfile() : mprecision(0.0) { // Public Methods ------------------------------------------------------------ +void PeakProfile::addRDFPeak(double* out, int ilo, int ihi, int rlosteps, + double rstep, double dist, double fwhm, + double scale) const { + for (int i = ilo; i < ihi; ++i) { + const double x = (rlosteps + i) * rstep - dist; + const double y = (*this)(x, fwhm); + // Contributions in G(r) need to be normalized by pair distance, + // not by r as done in PDFfit or PDFfit2. Here we rescale RDF + // in such way that division by r will give a correct result. + const double yrdf = y * (x / dist + 1); + out[i] += scale * yrdf; + } +} + void PeakProfile::setPrecision(double eps) { if (mprecision != eps) mticker.click(); mprecision = eps; diff --git a/src/diffpy/srreal/PeakProfile.hpp b/src/diffpy/srreal/PeakProfile.hpp index 0d0e20a..a25accd 100644 --- a/src/diffpy/srreal/PeakProfile.hpp +++ b/src/diffpy/srreal/PeakProfile.hpp @@ -47,6 +47,12 @@ class DLL_EXPORT PeakProfile : public diffpy::Attributes, virtual double operator()(double x, double fwhm) const = 0; virtual double xboundlo(double fwhm) const = 0; virtual double xboundhi(double fwhm) const = 0; + /// Add the peak of a pair at distance dist to an RDF on the r-grid + /// (rlosteps + i) * rstep. For i in [ilo, ihi) with x = r - dist: + /// out[i] += scale * (*this)(x, fwhm) * (x / dist + 1) + virtual void addRDFPeak(double* out, int ilo, int ihi, int rlosteps, + double rstep, double dist, double fwhm, + double scale) const; virtual void setPrecision(double eps); const double& getPrecision() const; virtual eventticker::EventTicker& ticker() const { return mticker; } diff --git a/src/diffpy/srreal/kernels/CalculationKernelAVX2.cpp b/src/diffpy/srreal/kernels/CalculationKernelAVX2.cpp new file mode 100644 index 0000000..4a53d21 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelAVX2.cpp @@ -0,0 +1,32 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * AVX2+FMA calculation kernel. Include nothing else here. + * + *****************************************************************************/ + +// xsimd namespace of this file, see CalculationKernelSIMD.ipp +#define LIBDIFFPY_XSIMD_NAMESPACE xsimd_avx2 +#include "CalculationKernelSIMD.ipp" + +namespace diffpy { +namespace srreal { +namespace kernels { + +const CalculationKernel& avx2CalculationKernel() { + static constexpr CalculationKernel kernel = + makeSIMDKernel>("avx2"); + return kernel; +} + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelAVX512.cpp b/src/diffpy/srreal/kernels/CalculationKernelAVX512.cpp new file mode 100644 index 0000000..5ed9d54 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelAVX512.cpp @@ -0,0 +1,32 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * AVX-512 calculation kernel. Include nothing else here. + * + *****************************************************************************/ + +// xsimd namespace of this file, see CalculationKernelSIMD.ipp +#define LIBDIFFPY_XSIMD_NAMESPACE xsimd_avx512 +#include "CalculationKernelSIMD.ipp" + +namespace diffpy { +namespace srreal { +namespace kernels { + +const CalculationKernel& avx512CalculationKernel() { + static constexpr CalculationKernel kernel = + makeSIMDKernel("avx512"); + return kernel; +} + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelDispatch.cpp b/src/diffpy/srreal/kernels/CalculationKernelDispatch.cpp new file mode 100644 index 0000000..7a94175 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelDispatch.cpp @@ -0,0 +1,113 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Runtime selection of the calculation kernel + * + *****************************************************************************/ + +#include +#include +#include +#include +#include + +#if defined(LIBDIFFPY_SIMD_KERNELS) +#include +#endif + +#include +#include + +namespace diffpy { +namespace srreal { + +using kernels::CalculationKernel; + +namespace { + +/// Kernels of this build that the CPU can run, from the most to the least +/// preferred. CMake defines LIBDIFFPY_KERNEL_ for each SIMD kernel it +/// compiles. The scalar kernel is last and always available. +const std::vector& availableKernels() { + static const std::vector available = [] { + std::vector rv; +#if defined(LIBDIFFPY_SIMD_KERNELS) + const auto cpu = xsimd::available_architectures(); +#endif +#if defined(LIBDIFFPY_KERNEL_AVX512) + // The AVX-512 file may also use AVX2 and FMA instructions. + if (cpu.avx512f && cpu.avx512cd && cpu.avx512dq && cpu.fma3_avx2) { + rv.push_back(&kernels::avx512CalculationKernel()); + } +#endif +#if defined(LIBDIFFPY_KERNEL_AVX2) + if (cpu.fma3_avx2) rv.push_back(&kernels::avx2CalculationKernel()); +#endif +#if defined(LIBDIFFPY_KERNEL_SSE2) + if (cpu.sse2) rv.push_back(&kernels::sse2CalculationKernel()); +#endif +#if defined(LIBDIFFPY_KERNEL_NEON) + if (cpu.neon64) rv.push_back(&kernels::neonCalculationKernel()); +#endif + rv.push_back(&kernels::scalarCalculationKernel()); + return rv; + }(); + return available; +} + +/// The active kernel, the most preferred one is selected on first use. +std::atomic& activeKernelSlot() { + static std::atomic slot{availableKernels().front()}; + return slot; +} + +} // namespace + +// Public functions ---------------------------------------------------------- + +std::vector availableCalculationKernels() { + std::vector rv; + for (const CalculationKernel* k : availableKernels()) + rv.emplace_back(k->name); + return rv; +} + +std::string getCalculationKernel() { + return kernels::activeCalculationKernel().name; +} + +void setCalculationKernel(const std::string& name) { + const auto& available = availableKernels(); + const auto it = std::ranges::find_if( + available, [&](const CalculationKernel* k) { return name == k->name; }); + if (it == available.end()) { + std::string names; + for (const CalculationKernel* k : available) { + if (!names.empty()) names += ", "; + names += k->name; + } + throw std::invalid_argument("Calculation kernel '" + name + + "' is not available, use one of: " + names + + "."); + } + activeKernelSlot().store(*it, std::memory_order_relaxed); +} + +namespace kernels { + +const CalculationKernel& activeCalculationKernel() { + return *activeKernelSlot().load(std::memory_order_relaxed); +} + +} // namespace kernels + +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp b/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp new file mode 100644 index 0000000..5a32f94 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp @@ -0,0 +1,111 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Internal interface of the calculation kernels, not installed. + * + * No standard library includes, see CalculationKernelSIMD.ipp. + * + *****************************************************************************/ + +#ifndef CALCULATIONKERNELIMPL_HPP_INCLUDED +#define CALCULATIONKERNELIMPL_HPP_INCLUDED + +#include + +namespace diffpy { +namespace srreal { +namespace kernels { + +/// Add one Gaussian peak to the PDFCalculator RDF array. +/// +/// For each i in [ilo, ihi) with x = (rlosteps + i) * rstep - dist: +/// out[i] += scale * G(x, fwhm) * (x / dist + 1) +/// where G(x, fwhm) = peakheight * exp(-4 * ln(2) * (x / fwhm)^2), or zero +/// for |x / fwhm| >= halfboundrel when cropped. peakheight is the profile +/// value at x = 0. Nothing is added for fwhm <= 0. +struct GaussianWindowArgs { + double* out; + int ilo; + int ihi; + int rlosteps; + double rstep; + double dist; + double fwhm; + double scale; + double peakheight; + bool cropped; + double halfboundrel; +}; + +/// Add one pair to the BaseDebyeSum F(Q) array. +/// +/// For kq = kqlo, kqlo + 1, ... while kq < kqhi: +/// q = kq * qstep +/// t = scale * exp(-0.5 * (dwsigma * q)^2) * sf0[kq] * sf1[kq] / dist +/// stop at the first kq where |t| <= prec +/// out[kq] += t * sin(q * dist) +struct DebyeSumArgs { + double* out; + int kqlo; + int kqhi; + double qstep; + double dist; + double dwsigma; + double scale; + double prec; + const double* sf0; + const double* sf1; +}; + +/// Select candidate bonds from n points by their squared distance. +/// +/// For each i in [0, n) with +/// d2 = (x[i] + ox)^2 + (y[i] + oy)^2 + (z[i] + oz)^2 +/// append i to idx when d2lo <= d2 <= d2hi or d2 is NaN. +/// Return the count. idx must hold n values. +/// Rounding of d2 differs between kernels. +struct BondCandidatesArgs { + const double* x; + const double* y; + const double* z; + int n; + double ox; + double oy; + double oz; + double d2lo; + double d2hi; + int* idx; +}; + +struct CalculationKernel { + /// Name for setCalculationKernel, a string literal. + const char* name; + void (*gaussianWindow)(const GaussianWindowArgs&); + void (*debyeSum)(const DebyeSumArgs&); + int (*bondCandidates)(const BondCandidatesArgs&); +}; + +/// The active kernel, exported for the unit tests. +DLL_EXPORT const CalculationKernel& activeCalculationKernel(); + +// Kernel tables per instruction set. The SIMD ones exist only when CMake +// compiles them, see cmake/CalculationKernels.cmake. +const CalculationKernel& scalarCalculationKernel(); +const CalculationKernel& sse2CalculationKernel(); +const CalculationKernel& avx2CalculationKernel(); +const CalculationKernel& avx512CalculationKernel(); +const CalculationKernel& neonCalculationKernel(); + +} // namespace kernels +} // namespace srreal +} // namespace diffpy + +#endif // CALCULATIONKERNELIMPL_HPP_INCLUDED diff --git a/src/diffpy/srreal/kernels/CalculationKernelNEON.cpp b/src/diffpy/srreal/kernels/CalculationKernelNEON.cpp new file mode 100644 index 0000000..e65a4c5 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelNEON.cpp @@ -0,0 +1,32 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * NEON calculation kernel. Include nothing else here. + * + *****************************************************************************/ + +// xsimd namespace of this file, see CalculationKernelSIMD.ipp +#define LIBDIFFPY_XSIMD_NAMESPACE xsimd_neon +#include "CalculationKernelSIMD.ipp" + +namespace diffpy { +namespace srreal { +namespace kernels { + +const CalculationKernel& neonCalculationKernel() { + static constexpr CalculationKernel kernel = + makeSIMDKernel("neon"); + return kernel; +} + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp b/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp new file mode 100644 index 0000000..bf2e5c4 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp @@ -0,0 +1,283 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * SIMD calculation kernel shared by the SSE2, AVX2, AVX-512 and NEON files. + * + * On the uniform grid exp and sin are advanced by recurrences instead of + * being called per point: + * g = exp(-b x^2): g(x + h) = g(x) rho(x), rho(x + h) = rho(x) exp(-2 b h^2) + * sin, cos of x+d: s' = s cos(d) + c sin(d), c' = c cos(d) - s sin(d) + * xsimd::exp and xsimd::sincos reseed them every RESEED_STEPS vector steps + * to limit the rounding error. Windows too far from the peak for the + * recurrence go to the scalar kernel. + * + * The linker keeps one copy of each inline function for the whole library, + * so a copy compiled here could run AVX-512 code on a CPU without it. To + * avoid that, each file renames the xsimd namespace to + * LIBDIFFPY_XSIMD_NAMESPACE and is compiled with optimization and without + * link-time optimization, which inlines the few standard library functions + * used here. For the same reason use only C math functions. + * + *****************************************************************************/ + +#if defined(__FAST_MATH__) || defined(_M_FP_FAST) +#error "The xsimd math functions give wrong results with fast-math." +#endif + +#if (defined(__GNUC__) || defined(__clang__)) && !defined(__OPTIMIZE__) +#error "Compile the SIMD kernels with optimization, see the comment above." +#endif + +#ifndef LIBDIFFPY_XSIMD_NAMESPACE +#error "Define LIBDIFFPY_XSIMD_NAMESPACE, see the comment above." +#endif + +#include +#include +#include + +#include + +// The xsimd headers and all code below use the renamed namespace. +#define xsimd LIBDIFFPY_XSIMD_NAMESPACE +#include + +namespace diffpy { +namespace srreal { +namespace kernels { + +namespace { + +/// Vector steps in one block of the recurrences. +constexpr int RESEED_STEPS = 32; + +/// Extra vector steps allowed in the last block. +constexpr int RESEED_SLACK = 8; + +/// Limit of b x^2 at the seeds and of |log(rho)|, which keeps the seeds +/// normal and rho finite. exp(-SAFE_EXPONENT) is about 1e-261. +constexpr double SAFE_EXPONENT = 600.0; + +/// True when xsimd::fma of Arch is a fused multiply-add. SSE2 has none and +/// xsimd evaluates u * v + w there. +template +constexpr bool FUSED_FMA = true; +template <> +constexpr bool FUSED_FMA = false; + +template +struct SIMDKernel { + using B = xsimd::batch; + using BB = xsimd::batch_bool; + static constexpr int W = static_cast(B::size); + + /// The batch {0, 1, ..., W - 1}. + static B laneIndex() { + double v[W]; + for (int j = 0; j < W; ++j) v[j] = j; + return B::load_unaligned(v); + } + + /// End of the block that starts at i, the last block ends at hi. + static int blockEnd(int i, int hi) { + const bool last = (hi - i <= (RESEED_STEPS + RESEED_SLACK) * W); + return last ? hi : i + RESEED_STEPS * W; + } + + /// Load the first n lanes from p, the other lanes are zero. + static B load(const double* p, int n) { + if (n >= W) return B::load_unaligned(p); + double v[W] = {}; + for (int j = 0; j < n; ++j) v[j] = p[j]; + return B::load_unaligned(v); + } + + /// Add the first n lanes of u * v to out. + static void addProduct(double* out, const B& u, const B& v, int n) { + if (n >= W) { + xsimd::fma(u, v, B::load_unaligned(out)).store_unaligned(out); + return; + } + double uv[W]; + (u * v).store_unaligned(uv); + for (int j = 0; j < n; ++j) out[j] += uv[j]; + } + + /// u * v + w for one value, rounded as in the lanes of xsimd::fma. + static double mulAdd(double u, double v, double w) { + return FUSED_FMA ? fma(u, v, w) : u * v + w; + } + + static void gaussianWindow(const GaussianWindowArgs& args) { + // Vector stores may alias anything, work on a copy of the arguments. + const GaussianWindowArgs a = args; + if (a.fwhm <= 0.0) return; + // exp(-b x^2) and its ratio over one vector step. + const double b = 4 * M_LN2 / (a.fwhm * a.fwhm); + const double wh = W * a.rstep; + // Use the scalar kernel for NaN, for windows too far from the peak, where + // b * x^2 or 2 * b * wh * |x| reach SAFE_EXPONENT, and for narrow peaks + // far out, where the grid rounding (~eps * dist) is too coarse. + const double xlo = (a.rlosteps + a.ilo) * a.rstep - a.dist; + const double xhi = (a.rlosteps + a.ihi) * a.rstep - a.dist; + const double xfar = (-xlo > xhi) ? -xlo : xhi; + const double reach = (xfar > 2 * wh) ? xfar : 2 * wh; + const bool recurrence = (b * xfar * reach < SAFE_EXPONENT) && + (DBL_EPSILON * a.dist <= 1e-13 * a.fwhm); + if (!recurrence) { + scalarCalculationKernel().gaussianWindow(args); + return; + } + const B negb(-b); + const B twowh(2 * wh); + const B whwh(wh * wh); + const B ratiostep(exp(-2 * b * wh * wh)); + const B vw(1.0 * W); + const B vrstep(a.rstep); + const B vdist(a.dist); + const B vinvdist(1.0 / a.dist); + const B vnorm(a.scale * a.peakheight); + const B vfwhm(a.fwhm); + const B vhalfbound(a.halfboundrel); + const B one(1.0); + const B zero(0.0); + // Grid steps of the lanes, whole numbers are exact in double. + B steps = B(double(a.rlosteps + a.ilo)) + laneIndex(); + // x / dist + 1 is r / dist, which grows by a constant per vector step. + const B vrscale(a.rstep / a.dist); + const B vrscalestep(wh / a.dist); + for (int i = a.ilo; i < a.ihi;) { + const int iend = blockEnd(i, a.ihi); + const B xseed = steps * vrstep - vdist; + B g = vnorm * xsimd::exp(negb * xseed * xseed); + B ratio = xsimd::exp(negb * (twowh * xseed + whwh)); + if (a.cropped) { + for (; i < iend; i += W) { + // Same x as the scalar kernel, so cropping agrees. + const B x = steps * vrstep - vdist; + const B y = + xsimd::select(xsimd::abs(x / vfwhm) >= vhalfbound, zero, g); + addProduct(a.out + i, y, xsimd::fma(x, vinvdist, one), iend - i); + g *= ratio; + ratio *= ratiostep; + steps += vw; + } + } else { + B rscaled = steps * vrscale; + const int istart = i; + for (; i < iend; i += W) { + addProduct(a.out + i, g, rscaled, iend - i); + g *= ratio; + ratio *= ratiostep; + rscaled += vrscalestep; + } + steps += B(double(i - istart)); + } + } + } + + static void debyeSum(const DebyeSumArgs& a) { + // Debye-Waller factor exp(-alpha kq^2) and its ratio over W steps. + const double dwstep = a.dwsigma * a.qstep; + const double alpha = 0.5 * dwstep * dwstep; + const B negalpha(-alpha); + const B twow(2.0 * W); + const B ww(1.0 * W * W); + const B ratiostep(exp(-2 * alpha * W * W)); + // Rotation of (sin, cos) of q * dist over W steps. + const double wangle = W * a.qstep * a.dist; + const B cosw(cos(wangle)); + const B sinw(sin(wangle)); + const B vqstep(a.qstep); + const B vdist(a.dist); + const B vdwsigma(a.dwsigma); + const B vsinescale(a.scale / a.dist); + const B vprec(a.prec); + const B neghalf(-0.5); + const B lanes = laneIndex(); + for (int kq = a.kqlo; kq < a.kqhi;) { + const int kqend = blockEnd(kq, a.kqhi); + const B k = B(double(kq)) + lanes; + const B q = k * vqstep; + const B dwq = vdwsigma * q; + // g includes the sine scale. + B g = vsinescale * xsimd::exp(neghalf * dwq * dwq); + B ratio = xsimd::exp(negalpha * (twow * k + ww)); + auto [s, c] = xsimd::sincos(q * vdist); + for (; kq < kqend; kq += W) { + const int n = kqend - kq; + // Unused lanes have zero t and do not reach the output. + B t = g * load(a.sf0 + kq, n) * load(a.sf1 + kq, n); + const auto small = xsimd::abs(t) <= vprec; + const bool last = xsimd::any(small); + if (last) { + // Keep the lanes before the first small term. + const uint64_t m = small.mask(); + t = xsimd::select(BB::from_mask((m ^ (m - 1)) >> 1), t, B(0.0)); + } + addProduct(a.out + kq, t, s, n); + if (last) return; + g *= ratio; + ratio *= ratiostep; + const B snext = xsimd::fma(s, cosw, c * sinw); + c = xsimd::fms(c, cosw, s * sinw); + s = snext; + } + } + } + + static int bondCandidates(const BondCandidatesArgs& a) { + const B vox(a.ox); + const B voy(a.oy); + const B voz(a.oz); + const B vlo(a.d2lo); + const B vhi(a.d2hi); + int m = 0; + int i = 0; + for (; i + W <= a.n; i += W) { + const B dx = B::load_unaligned(a.x + i) + vox; + const B dy = B::load_unaligned(a.y + i) + voy; + const B dz = B::load_unaligned(a.z + i) + voz; + const B d2 = xsimd::fma(dz, dz, xsimd::fma(dy, dy, dx * dx)); + // NaN fails both comparisons and stays a candidate. + const uint64_t bits = (!((d2 < vlo) || (d2 > vhi))).mask(); + if (!bits) continue; + for (int j = 0; j < W; ++j) { + a.idx[m] = i + j; + m += int((bits >> j) & 1); + } + } + // The last points go one by one with the rounding of the vector lanes, + // a partial vector load would be much slower. + for (; i < a.n; ++i) { + const double dx = a.x[i] + a.ox; + const double dy = a.y[i] + a.oy; + const double dz = a.z[i] + a.oz; + const double d2 = mulAdd(dz, dz, mulAdd(dy, dy, dx * dx)); + a.idx[m] = i; + m += int(!(d2 < a.d2lo || d2 > a.d2hi)); + } + return m; + } +}; + +/// Kernel table for the Arch instruction set. +template +constexpr CalculationKernel makeSIMDKernel(const char* name) { + using K = SIMDKernel; + return {name, K::gaussianWindow, K::debyeSum, K::bondCandidates}; +} + +} // namespace + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelSSE2.cpp b/src/diffpy/srreal/kernels/CalculationKernelSSE2.cpp new file mode 100644 index 0000000..3697aa1 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelSSE2.cpp @@ -0,0 +1,32 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * NEON calculation kernel. Include nothing else here. + * + *****************************************************************************/ + +// xsimd namespace of this file, see CalculationKernelSIMD.ipp +#define LIBDIFFPY_XSIMD_NAMESPACE xsimd_sse2 +#include "CalculationKernelSIMD.ipp" + +namespace diffpy { +namespace srreal { +namespace kernels { + +const CalculationKernel& sse2CalculationKernel() { + static constexpr CalculationKernel kernel = + makeSIMDKernel("sse2"); + return kernel; +} + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp b/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp new file mode 100644 index 0000000..5849f12 --- /dev/null +++ b/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp @@ -0,0 +1,75 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Scalar calculation kernel, the reference and fallback implementation + * + *****************************************************************************/ + +#include +#include + +#include + +namespace diffpy { +namespace srreal { +namespace kernels { + +namespace { + +void gaussianWindowScalar(const GaussianWindowArgs& a) { + if (a.fwhm <= 0.0) return; + for (int i = a.ilo; i < a.ihi; ++i) { + const double x = (a.rlosteps + i) * a.rstep - a.dist; + const double xrel = x / a.fwhm; + if (a.cropped && fabs(xrel) >= a.halfboundrel) continue; + const double y = a.peakheight * exp(-4 * std::numbers::ln2 * xrel * xrel); + // Normalized by dist instead of r, see PeakProfile::addRDFPeak. + const double yrdf = y * (x / a.dist + 1); + a.out[i] += a.scale * yrdf; + } +} + +void debyeSumScalar(const DebyeSumArgs& a) { + for (int kq = a.kqlo; kq < a.kqhi; ++kq) { + const double q = kq * a.qstep; + const double dwq = a.dwsigma * q; + const double dwscale = exp(-0.5 * dwq * dwq); + const double sinescale = a.scale * dwscale * a.sf0[kq] * a.sf1[kq] / a.dist; + if (fabs(sinescale) <= a.prec) break; + a.out[kq] += sinescale * sin(q * a.dist); + } +} + +int bondCandidatesScalar(const BondCandidatesArgs& a) { + int m = 0; + for (int i = 0; i < a.n; ++i) { + const double dx = a.x[i] + a.ox; + const double dy = a.y[i] + a.oy; + const double dz = a.z[i] + a.oz; + const double d2 = dx * dx + dy * dy + dz * dz; + a.idx[m] = i; + // NaN fails both comparisons and stays a candidate. + m += int(!(d2 < a.d2lo || d2 > a.d2hi)); + } + return m; +} + +} // namespace + +const CalculationKernel& scalarCalculationKernel() { + static constexpr CalculationKernel kernel = { + "scalar", gaussianWindowScalar, debyeSumScalar, bondCandidatesScalar}; + return kernel; +} + +} // namespace kernels +} // namespace srreal +} // namespace diffpy diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 158fecb..5f879eb 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -50,7 +50,13 @@ add_executable(alltests ${_suites} ${_helpers}) set_target_properties(alltests PROPERTIES CXX_EXTENSIONS OFF) if(MSVC) target_compile_options(alltests PRIVATE /bigobj) + set(_test_math_options /fp:precise) +else() + set(_test_math_options -fno-fast-math) endif() +# Numerical checks must recognize NaN and infinity even in Release builds. +set_property(SOURCE TestCalculationKernel.cpp test_helpers.cpp + APPEND PROPERTY COMPILE_OPTIONS "${_test_math_options}") target_include_directories(alltests PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(alltests PRIVATE libdiffpy::diffpy Catch2::Catch2WithMain) # test_helpers.cpp finds the testdata directory from this macro. diff --git a/src/tests/TestCalculationKernel.cpp b/src/tests/TestCalculationKernel.cpp new file mode 100644 index 0000000..faa7fc8 --- /dev/null +++ b/src/tests/TestCalculationKernel.cpp @@ -0,0 +1,487 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Tests of the calculation kernels against the scalar kernel + * + *****************************************************************************/ + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "test_helpers.hpp" + +// The internal kernel header is missing when testing an installed library. +#if __has_include() +#include +#define LIBDIFFPY_TEST_KERNEL_INTERNALS 1 +#endif + +using namespace diffpy::srreal; + +namespace { + +const double KERNEL_TOLERANCE = 1e-12; + +/// One calculation that is repeated with every kernel. +struct KernelCase { + std::string name; + std::string calculator; // "pdf" or "debye" + std::string structure; // .stru file in testdata + std::string profile; // peak profile, PDFCalculator only + std::vector> attributes; +}; + +const std::vector KERNEL_CASES = { + {"pdf_LiTaO3", "pdf", "LiTaO3.stru", "gaussian", {{"rmax", 15.0}}}, + {"pdf_LiTaO3_qmax", + "pdf", + "LiTaO3.stru", + "gaussian", + {{"rmin", 1.0}, {"qmax", 25.0}, {"delta2", 1.0}}}, + {"pdf_K2Bi8Se13_cropped", + "pdf", + "alpha_K2Bi8Se13.stru", + "croppedgaussian", + {{"rmax", 12.0}}}, + {"debye_LiTaO3", + "debye", + "LiTaO3.stru", + "", + {{"qmin", 1.0}, {"qmax", 20.0}, {"debyeprecision", 1e-8}}}, +}; + +template +void setAttributes(Calc& calc, const KernelCase& kc) { + for (const auto& [name, value] : kc.attributes) { + calc.setDoubleAttr(name, value); + } +} + +/// PDF and the second result array (RDF or F) for the active kernel. +std::vector computeCase(const KernelCase& kc) { + StructureAdapterPtr stru = loadTestPeriodicStructure(kc.structure); + if (kc.calculator == "pdf") { + PDFCalculator calc; + calc.setPeakProfileByType(kc.profile); + setAttributes(calc, kc); + calc.eval(stru); + return {calc.getPDF(), calc.getRDF()}; + } + DebyePDFCalculator calc; + setAttributes(calc, kc); + calc.eval(stru); + return {calc.getPDF(), calc.getF()}; +} + +} // namespace + +TEST_CASE("TestCalculationKernel::test_nonfinite_error", + "[CalculationKernel]") { + const double inf = std::numeric_limits::infinity(); + const double nan = std::numeric_limits::quiet_NaN(); + const QuantityType reference(2, 1.0); + for (double invalid : {nan, inf, -inf}) { + QuantityType actual(reference); + actual[0] = invalid; + CHECK(maxRelativeError(actual, reference) == inf); + CHECK(maxRelativeError(reference, actual) == inf); + CHECK(maxRelativeError(actual, actual) == inf); + } +} + +TEST_CASE("TestCalculationKernel::test_narrow_peaks", "[CalculationKernel]") { + KernelGuard guard; + for (double distance : {2.123, 100.123, 1000.123}) { + auto stru = std::make_shared(); + Atom atom; + atom.atomtype = "C"; + stru->append(atom); + atom.xyz_cartn[0] = distance; + stru->append(atom); + for (const char* profile : {"gaussian", "croppedgaussian"}) { + for (double width : {0.001, 0.01, 0.1}) { + PDFCalculator calc; + calc.setEvaluatorType(BASIC); + calc.setRmin(distance - 0.25); + calc.setRmax(distance + 0.25); + calc.setRstep(0.001); + calc.setPeakProfileByType(profile); + calc.setDoubleAttr("peakprecision", 1e-16); + calc.setPeakWidthModelByType("constant"); + calc.setDoubleAttr("width", width); + setCalculationKernel("scalar"); + calc.eval(stru); + const QuantityType reference = calc.getPDF(); + INFO("profile " << profile << ", distance " << distance << ", width " + << width); + testAllKernels([&] { + calc.eval(stru); + CHECK(maxRelativeError(calc.getPDF(), reference) <= KERNEL_TOLERANCE); + }); + } + } + } +} + +TEST_CASE("TestCalculationKernel::test_wide_windows", "[CalculationKernel]") { + // addRDFPeak over windows that reach far beyond the peak, where the + // Gaussian underflows to zero, must match the scalar kernel. + KernelGuard guard; + const int npts = 2000; + const double rstep = 0.01; + const double dist = 10.0; + const double fwhm = 0.05; + GaussianProfile gaussian; + CroppedGaussianProfile cropped; + const std::pair windows[] = {{0, npts}, {900, npts}, {0, 1010}, + {995, 1005}, {1500, npts}, {0, 600}}; + const PeakProfile* profiles[] = {&gaussian, &cropped}; + for (const PeakProfile* pk : profiles) { + for (const auto& [ilo, ihi] : windows) { + setCalculationKernel("scalar"); + QuantityType reference(npts, 0.0); + pk->addRDFPeak(reference.data(), ilo, ihi, 0, rstep, dist, fwhm, 1.0); + if (ilo <= 1000 && 1000 < ihi) { + // x is zero at the point 1000, where y equals the peak height + CHECK(reference[1000] == (*pk)(0.0, fwhm)); + } + INFO("profile " << pk->type() << ", window [" << ilo << ", " << ihi + << ")"); + testAllKernels([&] { + QuantityType actual(npts, 0.0); + pk->addRDFPeak(actual.data(), ilo, ihi, 0, rstep, dist, fwhm, 1.0); + // relative to the peak height, the reference may be all zero + CHECK(maxRelativeError(actual, reference, (*pk)(0.0, fwhm)) <= + KERNEL_TOLERANCE); + }); + } + } +} + +namespace { + +/// Gaussian profile with a changed operator(), which must not use the kernel. +class DoubledGaussianProfile : public GaussianProfile { + public: + double operator()(double x, double fwhm) const override { + return 2 * this->GaussianProfile::operator()(x, fwhm); + } +}; + +/// Cropped profile with a changed operator(), which must not use the kernel. +class DoubledCroppedGaussianProfile : public CroppedGaussianProfile { + public: + double operator()(double x, double fwhm) const override { + return 2 * this->CroppedGaussianProfile::operator()(x, fwhm); + } +}; + +/// Gaussian profile that keeps operator() and asks for the kernel. +class KernelGaussianProfile : public GaussianProfile { + protected: + KernelShape kernelShape() const override { return KernelShape::Gaussian; } +}; + +QuantityType rdfPeak(const PeakProfile& pk) { + QuantityType rv(300, 0.0); + pk.addRDFPeak(rv.data(), 50, 250, 900, 0.01, 10.5, 0.37, 1.5); + return rv; +} + +} // namespace + +TEST_CASE("TestCalculationKernel::test_derived_profiles", + "[CalculationKernel]") { + testAllKernels([] { + const QuantityType gaussian = rdfPeak(GaussianProfile()); + const QuantityType cropped = rdfPeak(CroppedGaussianProfile()); + QuantityType doubled(gaussian); + for (double& y : doubled) y *= 2; + QuantityType doubledcropped(cropped); + for (double& y : doubledcropped) y *= 2; + CHECK(maxRelativeError(rdfPeak(DoubledGaussianProfile()), doubled) <= + KERNEL_TOLERANCE); + CHECK(maxRelativeError(rdfPeak(DoubledCroppedGaussianProfile()), + doubledcropped) <= KERNEL_TOLERANCE); + CHECK(rdfPeak(KernelGaussianProfile()) == gaussian); + }); +} + +TEST_CASE("TestCalculationKernel::test_available", "[CalculationKernel]") { + const std::vector names = availableCalculationKernels(); + const auto has = [&](const char* name) { + return std::ranges::find(names, name) != names.end(); + }; + REQUIRE(!names.empty()); + // The scalar kernel is the last resort. + CHECK(names.back() == "scalar"); + const std::vector known = {"avx512", "avx2", "sse2", "neon", + "scalar"}; + for (const std::string& kname : names) { + INFO("kernel " << kname); + CHECK(std::ranges::count(names, kname) == 1); + CHECK(std::ranges::find(known, kname) != known.end()); + } + // The x86-64 kernels are built together and wider ones need newer CPUs. + if (has("avx512")) { + CHECK(has("avx2")); + } + if (has("avx2")) { + CHECK(has("sse2")); + } + CHECK(!(has("sse2") && has("neon"))); +} + +TEST_CASE("TestCalculationKernel::test_selection", "[CalculationKernel]") { + KernelGuard guard; + const std::vector names = availableCalculationKernels(); + CHECK(std::ranges::find(names, getCalculationKernel()) != names.end()); + for (const std::string& kname : names) { + setCalculationKernel(kname); + CHECK(kname == getCalculationKernel()); + } + // Names are lowercase and kernels of other CPUs are not available. + const std::string active = getCalculationKernel(); + for (const char* invalid : {"", "sse", "avx", "AVX2", "Scalar", "sve"}) { + INFO("kernel '" << invalid << "'"); + CHECK_THROWS_AS(setCalculationKernel(invalid), std::invalid_argument); + CHECK(active == getCalculationKernel()); + } + if (std::ranges::find(names, "sse2") != names.end()) { + CHECK_THROWS_AS(setCalculationKernel("neon"), std::invalid_argument); + } else { + CHECK_THROWS_AS(setCalculationKernel("sse2"), std::invalid_argument); + } +} + +TEST_CASE("TestCalculationKernel::test_match_scalar", "[CalculationKernel]") { + KernelGuard guard; + for (const KernelCase& kc : KERNEL_CASES) { + setCalculationKernel("scalar"); + const std::vector reference = computeCase(kc); + INFO("case " << kc.name); + testAllKernels([&] { + const std::vector actual = computeCase(kc); + REQUIRE(actual.size() == reference.size()); + for (size_t i = 0; i < actual.size(); ++i) { + INFO("array " << i); + CHECK(maxRelativeError(actual[i], reference[i]) <= KERNEL_TOLERANCE); + } + }); + } +} + +#ifdef LIBDIFFPY_TEST_KERNEL_INTERNALS +namespace { + +/// True when both arrays keep the fill value at the same indices. +bool sameUntouched(const QuantityType& actual, const QuantityType& reference, + double fill) { + if (actual.size() != reference.size()) return false; + for (size_t i = 0; i < actual.size(); ++i) { + if ((actual[i] == fill) != (reference[i] == fill)) return false; + } + return true; +} + +} // namespace + +TEST_CASE("TestCalculationKernel::test_window_lengths", "[CalculationKernel]") { + using namespace diffpy::srreal::kernels; + KernelGuard guard; + setCalculationKernel("scalar"); + const CalculationKernel& scalar = activeCalculationKernel(); + // around the vector and the reseed block sizes of all kernels, + // including the longest last block of 80, 160 or 320 points + const std::vector lengths = {0, 1, 2, 3, 4, 5, 7, 8, 9, + 15, 16, 17, 63, 64, 65, 79, 80, 81, + 127, 128, 129, 159, 160, 161, 255, 256, 257, + 300, 319, 320, 321, 511, 512, 513, 600}; + const double fill = 0.25; + const int padding = 20; + testAllKernels([&] { + const CalculationKernel& kernel = activeCalculationKernel(); + for (int n : lengths) { + for (int lo : {0, 1, 3}) { + const int hi = lo + n; + // one Gaussian peak across the window + for (bool cropped : {false, true}) { + INFO("lo " << lo << ", n " << n << ", cropped " << cropped); + QuantityType reference(hi + padding, fill); + QuantityType actual(reference); + const int rlosteps = 100; + const double rstep = 0.01; + GaussianWindowArgs args = { + .out = reference.data(), + .ilo = lo, + .ihi = hi, + .rlosteps = rlosteps, + .rstep = rstep, + .dist = (rlosteps + lo + 0.5 * n + 0.37) * rstep, + .fwhm = std::max(n, 8) * rstep / 4, + .scale = -1.5, + .peakheight = 1.0, + .cropped = cropped, + .halfboundrel = 1.5, + }; + scalar.gaussianWindow(args); + args.out = actual.data(); + kernel.gaussianWindow(args); + CHECK(maxRelativeError(actual, reference) <= KERNEL_TOLERANCE); + CHECK(sameUntouched(actual, reference, fill)); + } + // Debye sum over the whole range and stopped at kq = 298 or 115. + // The scale is negative when a calculator removes a pair. + std::vector sf(hi); + for (int kq = 0; kq < hi; ++kq) sf[kq] = std::exp(-0.01 * kq); + for (double scale : {2.0, -2.0}) { + for (double prec : {0.0, 1e-3, 0.1}) { + INFO("lo " << lo << ", n " << n << ", scale " << scale << ", prec " + << prec); + QuantityType reference(hi + padding, fill); + QuantityType actual(reference); + DebyeSumArgs args = { + .out = reference.data(), + .kqlo = lo, + .kqhi = hi, + .qstep = 0.05, + .dist = 1.7, + .dwsigma = 0.1, + .scale = scale, + .prec = prec, + .sf0 = sf.data(), + .sf1 = sf.data(), + }; + scalar.debyeSum(args); + args.out = actual.data(); + kernel.debyeSum(args); + CHECK(maxRelativeError(actual, reference) <= KERNEL_TOLERANCE); + CHECK(sameUntouched(actual, reference, fill)); + } + } + } + } + }); +} + +TEST_CASE("TestCalculationKernel::test_bond_candidates", + "[CalculationKernel]") { + const double nan = std::numeric_limits::quiet_NaN(); + const double inf = std::numeric_limits::infinity(); + const double d2lo = 0.25; + const double d2hi = 4.0; + // squares of 0.5 and 2.0 are the exact bounds + const std::vector pattern = {0.5, + std::nextafter(0.5, 0.0), + std::nextafter(0.5, 1.0), + 2.0, + std::nextafter(2.0, 1.0), + std::nextafter(2.0, 3.0), + nan, + -1.0, + 0.0, + 3.0, + inf, + -0.5, + nan}; + const std::vector lengths = {0, 1, 2, 3, 4, 5, 7, + 8, 9, 15, 16, 17, 31, 100}; + testAllKernels([&] { + const diffpy::srreal::kernels::CalculationKernel& kernel = + diffpy::srreal::kernels::activeCalculationKernel(); + for (int n : lengths) { + for (int shift = 0; shift < 3; ++shift) { + for (int axis = 0; axis < 3; ++axis) { + INFO("n = " << n << ", shift " << shift << ", axis " << axis); + std::vector xyz[3]; + for (std::vector& c : xyz) c.assign(n, 0.0); + std::vector expected; + for (int i = 0; i < n; ++i) { + const double v = pattern[(i + shift) % pattern.size()]; + xyz[axis][i] = v; + const double d2 = v * v; + if (std::isnan(d2) || (d2lo <= d2 && d2 <= d2hi)) { + expected.push_back(i); + } + } + std::vector idx(n + 1, -1); + diffpy::srreal::kernels::BondCandidatesArgs args; + args.x = xyz[0].data(); + args.y = xyz[1].data(); + args.z = xyz[2].data(); + args.n = n; + args.ox = 0.0; + args.oy = 0.0; + args.oz = 0.0; + args.d2lo = d2lo; + args.d2hi = d2hi; + args.idx = idx.data(); + const int m = kernel.bondCandidates(args); + REQUIRE(m == int(expected.size())); + CHECK(std::vector(idx.begin(), idx.begin() + m) == expected); + CHECK(-1 == idx[n]); + } + } + } + // All points are the same, so all or none are candidates. With the + // bound between the fused and the separately rounded x^2 + z^2, the + // count is partial when the last points round unlike the vector lanes. + int nrounded = 0; + for (int k = 1; k < 200; ++k) { + const double vx = 1.0 + k * 0x1p-30; + const double vz = 3.0 - k * 0x1p-31; + const volatile double xx = vx * vx; + const volatile double zz = vz * vz; + const double fused = std::fma(vz, vz, xx); + const double separate = zz + xx; + if (fused == separate) continue; + ++nrounded; + for (int n : {9, 15, 17}) { + INFO("k = " << k << ", n = " << n); + std::vector x(n, vx), y(n, 0.0), z(n, vz); + std::vector idx(n); + const diffpy::srreal::kernels::BondCandidatesArgs args = { + .x = x.data(), + .y = y.data(), + .z = z.data(), + .n = n, + .d2hi = std::min(fused, separate), + .idx = idx.data(), + }; + const int m = kernel.bondCandidates(args); + CHECK((m == 0 || m == n)); + } + } + CHECK(nrounded > 0); + std::vector x(9, 1.0), y(9, -2.0), z(9, 2.0); + std::vector idx(9); + diffpy::srreal::kernels::BondCandidatesArgs args = { + x.data(), y.data(), z.data(), 9, -1.0, 2.0, 0.0, d2lo, d2hi, idx.data()}; + CHECK(9 == kernel.bondCandidates(args)); + // z + oz is exactly the next double above 2 + args.oz = std::nextafter(2.0, 3.0) - 2.0; + CHECK(0 == kernel.bondCandidates(args)); + }); +} +#endif // LIBDIFFPY_TEST_KERNEL_INTERNALS diff --git a/src/tests/test_helpers.cpp b/src/tests/test_helpers.cpp index 33c37ba..ab9dfe9 100644 --- a/src/tests/test_helpers.cpp +++ b/src/tests/test_helpers.cpp @@ -20,7 +20,9 @@ #include #include #include +#include #include +#include #include #include @@ -135,4 +137,22 @@ diffpy::srreal::StructureAdapterPtr loadTestPeriodicStructure( return pstru; } +// Array comparison ---------------------------------------------------------- + +double maxRelativeError(const diffpy::srreal::QuantityType& actual, + const diffpy::srreal::QuantityType& reference, + double scale) { + const double inf = std::numeric_limits::infinity(); + if (actual.size() != reference.size()) return inf; + double maxdiff = 0.0; + double maxref = 0.0; + for (size_t i = 0; i < actual.size(); ++i) { + if (!std::isfinite(actual[i]) || !std::isfinite(reference[i])) return inf; + maxdiff = std::max(maxdiff, std::fabs(actual[i] - reference[i])); + maxref = std::max(maxref, std::fabs(reference[i])); + } + if (scale == 0.0) scale = maxref; + return (scale > 0.0) ? maxdiff / scale : inf; +} + // End of test_helpers.cpp diff --git a/src/tests/test_helpers.hpp b/src/tests/test_helpers.hpp index cb067a0..16ab01f 100644 --- a/src/tests/test_helpers.hpp +++ b/src/tests/test_helpers.hpp @@ -20,6 +20,10 @@ #define TEST_HELPERS_HPP_INCLUDED #include +#include +#include +#include +#include #include std::string prepend_tests_dir(const std::string& f); @@ -28,4 +32,36 @@ std::string prepend_testdata_dir(const std::string& f); diffpy::srreal::StructureAdapterPtr loadTestPeriodicStructure( const std::string& tailname); +/// max|actual - reference| / scale, where scale defaults to max|reference|. +/// Infinity for different sizes, a value that is not finite or zero scale. +double maxRelativeError(const diffpy::srreal::QuantityType& actual, + const diffpy::srreal::QuantityType& reference, + double scale = 0.0); + +// Calculation kernels ------------------------------------------------------- + +/// Restore the kernel that was active when the test started. +class KernelGuard { + public: + KernelGuard() : msaved(diffpy::srreal::getCalculationKernel()) {} + ~KernelGuard() { diffpy::srreal::setCalculationKernel(msaved); } + + private: + std::string msaved; +}; + +/// Call fn() with each available calculation kernel active, then restore the +/// kernel that was active before. Failed checks report their kernel. +/// Do not use SECTION in fn, it would run only with the first kernel. +template +void testAllKernels(Fn&& fn) { + KernelGuard guard; + for (const std::string& kname : + diffpy::srreal::availableCalculationKernels()) { + diffpy::srreal::setCalculationKernel(kname); + INFO("kernel " << kname); + fn(); + } +} + #endif // TEST_HELPERS_HPP_INCLUDED From e030b64522f73131c64d378f31dc45c264adbc69 Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Tue, 10 Nov 2026 21:48:56 -0800 Subject: [PATCH 3/6] perf: batch periodic bonds and reuse symmetry positions --- news/faster-pdf.rst | 13 +- src/diffpy/srreal/BaseBondGenerator.cpp | 25 +- src/diffpy/srreal/BaseBondGenerator.hpp | 5 + src/diffpy/srreal/CrystalStructureAdapter.cpp | 157 +++-- src/diffpy/srreal/CrystalStructureAdapter.hpp | 31 +- src/diffpy/srreal/PQEvaluator.cpp | 4 + .../srreal/PeriodicStructureAdapter.cpp | 453 ++++++++++++-- .../srreal/PeriodicStructureAdapter.hpp | 149 ++++- src/diffpy/srreal/StructureAdapter.cpp | 32 +- src/tests/TestBondGeneratorBatching.cpp | 575 ++++++++++++++++++ src/tests/TestPQEvaluator.cpp | 121 ++++ 11 files changed, 1407 insertions(+), 158 deletions(-) create mode 100644 src/tests/TestBondGeneratorBatching.cpp diff --git a/news/faster-pdf.rst b/news/faster-pdf.rst index 10b98bf..51964e7 100644 --- a/news/faster-pdf.rst +++ b/news/faster-pdf.rst @@ -25,12 +25,18 @@ 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. +* 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. **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:** @@ -41,3 +47,8 @@ * ``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. diff --git a/src/diffpy/srreal/BaseBondGenerator.cpp b/src/diffpy/srreal/BaseBondGenerator.cpp index 9e1abe4..919dc64 100644 --- a/src/diffpy/srreal/BaseBondGenerator.cpp +++ b/src/diffpy/srreal/BaseBondGenerator.cpp @@ -17,6 +17,8 @@ * *****************************************************************************/ +#include + #include #include #include @@ -38,7 +40,10 @@ BaseBondGenerator::BaseBondGenerator(StructureAdapterConstPtr stru) mr0(R3::zerovector), mr1(R3::zerovector), mr01(R3::zerovector), - mdistance(0.0) { + mdistance(0.0), + msite_consecutive(true), + mmultiplicity_site(-1), + mmultiplicity(0) { int cnt = stru->countSites(); msite_all.resize(cnt); for (int i = 0; i != cnt; ++i) msite_all[i] = i; @@ -72,6 +77,7 @@ void BaseBondGenerator::selectAnchorSite(int anchor) { assert(0 <= anchor && anchor < mstructure->countSites()); msite_anchor = anchor; mr0 = mstructure->siteCartesianPosition(msite_anchor); + mmultiplicity_site = -1; this->setFinishedFlag(); } @@ -80,20 +86,21 @@ void BaseBondGenerator::selectSiteRange(int first, int last) { assert(last <= mstructure->countSites()); msite_first = msite_all.begin() + first; msite_last = msite_all.begin() + last; + msite_consecutive = true; this->setFinishedFlag(); } void BaseBondGenerator::selectSites(const SiteIndices& selection) { msite_selection = selection; - msite_first = msite_selection.begin(); - msite_last = msite_selection.end(); - this->setFinishedFlag(); + this->selectSites(msite_selection.begin(), msite_selection.end()); } void BaseBondGenerator::selectSites(SiteIndices::const_iterator first, SiteIndices::const_iterator last) { msite_first = first; msite_last = last; + const auto gap = [](int i, int j) { return j != i + 1; }; + msite_consecutive = (std::adjacent_find(first, last, gap) == last); this->setFinishedFlag(); } @@ -118,7 +125,12 @@ int BaseBondGenerator::site0() const { return msite_anchor; } int BaseBondGenerator::site1() const { return *msite_current; } int BaseBondGenerator::multiplicity() const { - return mstructure->siteMultiplicity(this->site0()); + // Calculators ask for every bond, look it up once per anchor. + if (mmultiplicity_site != this->site0()) { + mmultiplicity = mstructure->siteMultiplicity(this->site0()); + mmultiplicity_site = this->site0(); + } + return mmultiplicity; } const R3::Vector& BaseBondGenerator::r0() const { return mr0; } @@ -138,8 +150,7 @@ const R3::Matrix& BaseBondGenerator::Ucartesian1() const { } double BaseBondGenerator::msd() const { - static R3::Vector s; - s = this->r01(); + const R3::Vector& s = this->r01(); double msd0 = meanSquareDisplacement( this->Ucartesian0(), s, mstructure->siteAnisotropy(this->site0())); double msd1 = meanSquareDisplacement( diff --git a/src/diffpy/srreal/BaseBondGenerator.hpp b/src/diffpy/srreal/BaseBondGenerator.hpp index 64133f3..62080af 100644 --- a/src/diffpy/srreal/BaseBondGenerator.hpp +++ b/src/diffpy/srreal/BaseBondGenerator.hpp @@ -81,6 +81,11 @@ class DLL_EXPORT BaseBondGenerator { double mdistance; SiteIndices msite_all; SiteIndices msite_selection; + /// True when the selected sites are consecutive indices. + bool msite_consecutive; + /// Multiplicity of the anchor site mmultiplicity_site, -1 when unknown. + mutable int mmultiplicity_site; + mutable int mmultiplicity; // methods virtual bool iterateSymmetry(); diff --git a/src/diffpy/srreal/CrystalStructureAdapter.cpp b/src/diffpy/srreal/CrystalStructureAdapter.cpp index 90da49c..38a5a79 100644 --- a/src/diffpy/srreal/CrystalStructureAdapter.cpp +++ b/src/diffpy/srreal/CrystalStructureAdapter.cpp @@ -20,7 +20,12 @@ * *****************************************************************************/ +#include #include +#include +#include +#include +#include #include #include @@ -37,6 +42,31 @@ namespace srreal { const double DEFAULT_SYMMETRY_PRECISION = 5e-5; +// Local Helpers ------------------------------------------------------------- + +namespace { + +// Bitwise comparisons detect stale symmetry positions. Unlike == they find +// NaN equal to itself, so NaN does not force an expansion on every check. + +bool sameBits(std::span x, std::span y) { + return std::ranges::equal(std::as_bytes(x), std::as_bytes(y)); +} + +bool sameAtomBits(const Atom& a0, const Atom& a1) { + return a0.atomtype == a1.atomtype && a0.anisotropy == a1.anisotropy && + sameBits({&a0.occupancy, 1}, {&a1.occupancy, 1}) && + sameBits(a0.xyz_cartn.data(), a1.xyz_cartn.data()) && + sameBits(a0.uij_cartn.data(), a1.uij_cartn.data()); +} + +/// Lattice equality is defined by its base vectors. +bool sameLatticeBits(const Lattice& L0, const Lattice& L1) { + return sameBits(L0.base().data(), L1.base().data()); +} + +} // namespace + ////////////////////////////////////////////////////////////////////////////// // class CrystalStructureAdapter ////////////////////////////////////////////////////////////////////////////// @@ -56,15 +86,14 @@ StructureAdapterPtr CrystalStructureAdapter::clone() const { } BaseBondGeneratorPtr CrystalStructureAdapter::createBondGenerator() const { - this->updateSymmetryPositions(); + // the generator brings the symmetry positions up to date BaseBondGeneratorPtr bnds( new CrystalStructureBondGenerator(shared_from_this())); return bnds; } int CrystalStructureAdapter::siteMultiplicity(int idx) const { - if (!this->isSymmetryCached()) this->updateSymmetryPositions(); - int rv = msymatoms[idx].size(); + int rv = this->getEquivalentAtoms(idx).size(); return rv; } @@ -122,7 +151,7 @@ const SymOpRotTrans& CrystalStructureAdapter::getSymOp(int i) const { const CrystalStructureAdapter::AtomVector& CrystalStructureAdapter::getEquivalentAtoms(int idx) const { assert(0 <= idx && idx < this->countSites()); - if (!this->isSymmetryCached()) this->updateSymmetryPositions(); + if (!this->isSiteSymmetryUpToDate(idx)) this->updateSymmetryPositions(); return msymatoms[idx]; } @@ -167,11 +196,9 @@ CrystalStructureAdapter::AtomVector CrystalStructureAdapter::expandLatticeAtom( // calculate mean values from equivalent sites and adjust any roundoffs assert(eqsites.size() == eqduplicity.size()); assert(eqsites.size() == eqsumpos.size()); - iterator ai = eqsites.begin(); - vector::const_iterator sii = eqsumpos.begin(); - vector::const_iterator dpi = eqduplicity.begin(); - for (; ai != eqsites.end(); ++ai, ++sii, ++dpi) { - ai->xyz_cartn = (*sii) / (*dpi); + for (auto&& [atom, sumpos, duplicity] : + std::views::zip(eqsites, eqsumpos, eqduplicity)) { + atom.xyz_cartn = sumpos / duplicity; } return eqsites; } @@ -191,6 +218,8 @@ void CrystalStructureAdapter::updateSymmetryPositions() const { iterator ai = saii->begin(); for (; ai != saii->end(); ++ai) this->toCartesian(*ai); } + msymsource.assign(this->begin(), this->end()); + msymlattice = this->getLattice(); msymmetry_cached = true; } @@ -200,12 +229,20 @@ int CrystalStructureAdapter::findEqualPosition(const AtomVector& eqsites, const Atom& a0) const { const double symeps = this->getSymmetryPrecision(); const Lattice& L = this->getLattice(); + // Over a distance d a fractional coordinate changes by at most d times its + // reciprocal vector length. Twice that bound skips most sites early. + const double xtol = 2 * symeps * L.ar(); + const double ytol = 2 * symeps * L.br(); R3::Vector dxyz; const_iterator ai = eqsites.begin(); for (; ai != eqsites.end(); ++ai) { - dxyz = ai->xyz_cartn - a0.xyz_cartn; + dxyz[0] = ai->xyz_cartn[0] - a0.xyz_cartn[0]; dxyz[0] -= round(dxyz[0]); + if (fabs(dxyz[0]) > xtol) continue; + dxyz[1] = ai->xyz_cartn[1] - a0.xyz_cartn[1]; dxyz[1] -= round(dxyz[1]); + if (fabs(dxyz[1]) > ytol) continue; + dxyz[2] = ai->xyz_cartn[2] - a0.xyz_cartn[2]; dxyz[2] -= round(dxyz[2]); if (L.norm(dxyz) <= symeps) return (ai - eqsites.begin()); } @@ -218,6 +255,27 @@ bool CrystalStructureAdapter::isSymmetryCached() const { return msymmetry_cached; } +bool CrystalStructureAdapter::isSiteSymmetryUpToDate(int idx) const { + // Changes of the symmetry operations or precision reset msymmetry_cached, + // so only the atom and the lattice need a check. + return this->isSymmetryCached() && msymsource.size() == msymatoms.size() && + sameLatticeBits(msymlattice, this->getLattice()) && + sameAtomBits(msymsource[idx], (*this)[idx]); +} + +bool CrystalStructureAdapter::isSymmetryUpToDate() const { + return this->isSymmetryCached() && msymsource.size() == msymatoms.size() && + sameLatticeBits(msymlattice, this->getLattice()) && + std::equal(this->begin(), this->end(), msymsource.begin(), + sameAtomBits); +} + +const std::vector& +CrystalStructureAdapter::symmetryAtoms() const { + if (!this->isSymmetryUpToDate()) this->updateSymmetryPositions(); + return msymatoms; +} + // Comparison functions ------------------------------------------------------ bool operator==(const CrystalStructureAdapter& stru0, @@ -240,70 +298,35 @@ bool operator!=(const CrystalStructureAdapter& stru0, // Constructor --------------------------------------------------------------- +// The positions of site1 are its symmetry equivalent atoms. They are the +// outer loop of the bond order, see PeriodicStructureBondGenerator. + CrystalStructureBondGenerator::CrystalStructureBondGenerator( StructureAdapterConstPtr adpt) - : PeriodicStructureBondGenerator(adpt) { - mcstructure = dynamic_cast(adpt.get()); - assert(mcstructure); - msymidx = 0; - mpuc1 = &(R3::zeromatrix()); + : PeriodicStructureBondGenerator(adpt, BondOrder::PositionFirst) { + const CrystalStructureAdapter* cstru = + dynamic_cast(adpt.get()); + if (!cstru) { + throw invalid_argument( + "CrystalStructureBondGenerator requires a CrystalStructureAdapter."); + } + // Copy the symmetry positions and displacements, the adapter replaces + // them whenever it expands the symmetry again. + const vector& symatoms = + cstru->symmetryAtoms(); + for (int i = 0; i < int(symatoms.size()); ++i) { + for (const Atom& a : symatoms[i]) { + this->addPosition(i, a.xyz_cartn); + mcelluij.push_back(a.uij_cartn); + } + } } // Public Methods ------------------------------------------------------------ -void CrystalStructureBondGenerator::selectAnchorSite(int anchor) { - this->BaseBondGenerator::selectAnchorSite(anchor); - const Atom& a0 = this->symatoms(anchor)[0]; - mr0 = a0.xyz_cartn; -} - const R3::Matrix& CrystalStructureBondGenerator::Ucartesian1() const { - return *mpuc1; -} - -// Protected Methods --------------------------------------------------------- - -bool CrystalStructureBondGenerator::iterateSymmetry() { - // Iterate the sphere at a fixed symmetry position. - if (this->PeriodicStructureBondGenerator::iterateSymmetry()) { - this->updater1(); - return true; - } - // Advance to the next symmetry position. We are done if they - // were all already used. - const AtomVector& sa = this->symatoms(this->site1()); - ++msymidx; - if (msymidx >= sa.size()) return false; - // rewind the sphere for the new symmetry position. - // this calls CrystalStructureBondGenerator::updater1() - this->PeriodicStructureBondGenerator::rewindSymmetry(); - return true; -} - -void CrystalStructureBondGenerator::rewindSymmetry() { - msymidx = 0; - // this calls CrystalStructureBondGenerator::updater1() - this->PeriodicStructureBondGenerator::rewindSymmetry(); -} - -void CrystalStructureBondGenerator::getNextBond() { - this->BaseBondGenerator::getNextBond(); -} - -void CrystalStructureBondGenerator::updater1() { - const AtomVector& sa = this->symatoms(this->site1()); - assert(msymidx < sa.size()); - mr1 = mrcsphere + sa[msymidx].xyz_cartn; - mpuc1 = &(sa[msymidx].uij_cartn); - this->updateDistance(); -} - -// Private Methods ----------------------------------------------------------- - -const CrystalStructureAdapter::AtomVector& -CrystalStructureBondGenerator::symatoms(int idx) { - assert(0 <= idx && idx < int(mcstructure->msymatoms.size())); - return mcstructure->msymatoms[idx]; + if (mcelluij.empty()) return R3::zeromatrix(); + return mcelluij[mcellindex]; } } // namespace srreal diff --git a/src/diffpy/srreal/CrystalStructureAdapter.hpp b/src/diffpy/srreal/CrystalStructureAdapter.hpp index 0a502c3..34f9246 100644 --- a/src/diffpy/srreal/CrystalStructureAdapter.hpp +++ b/src/diffpy/srreal/CrystalStructureAdapter.hpp @@ -90,6 +90,9 @@ class DLL_EXPORT CrystalStructureAdapter : public PeriodicStructureAdapter { double msymmetry_precision; mutable std::vector msymatoms; mutable bool msymmetry_cached; + /// Atoms and lattice that msymatoms were expanded from, not serialized. + mutable AtomVector msymsource; + mutable Lattice msymlattice; // symmetry helpers /// return index of AtomVector atom at an equal position or -1 @@ -98,6 +101,12 @@ class DLL_EXPORT CrystalStructureAdapter : public PeriodicStructureAdapter { /// this only detects addition or removal of atom in the asymmetric /// unit, but does not check for changes in atom positions. bool isSymmetryCached() const; + /// exact check of site idx, also for changes in place + bool isSiteSymmetryUpToDate(int idx) const; + /// exact check of all sites + bool isSymmetryUpToDate() const; + /// symmetry positions of all sites brought up to date + const std::vector& symmetryAtoms() const; // comparison friend DLL_EXPORT bool operator==(const CrystalStructureAdapter&, @@ -130,29 +139,13 @@ class DLL_EXPORT CrystalStructureBondGenerator // constructors CrystalStructureBondGenerator(StructureAdapterConstPtr); - // configuration - virtual void selectAnchorSite(int); - // data access virtual const R3::Matrix& Ucartesian1() const; - protected: - // methods - virtual bool iterateSymmetry(); - virtual void rewindSymmetry(); - virtual void getNextBond(); - virtual void updater1(); - - // data - const CrystalStructureAdapter* mcstructure; - size_t msymidx; - const R3::Matrix* mpuc1; - private: - typedef CrystalStructureAdapter::AtomVector AtomVector; - - // methods - const AtomVector& symatoms(int idx); + // data + /// Displacement tensor at each of the site1 positions. + std::vector mcelluij; }; } // namespace srreal diff --git a/src/diffpy/srreal/PQEvaluator.cpp b/src/diffpy/srreal/PQEvaluator.cpp index 3e67e75..aaa567b 100644 --- a/src/diffpy/srreal/PQEvaluator.cpp +++ b/src/diffpy/srreal/PQEvaluator.cpp @@ -73,6 +73,8 @@ void PQEvaluatorBasic::validate(PairQuantity& pq) const { return; } void PQEvaluatorBasic::updateValue(PairQuantity& pq, StructureAdapterPtr stru) { mtypeused = BASIC; + // The value is not valid until this update completes. + mvalue_ticker = eventticker::EventTicker(); pq.setStructure(stru); BaseBondGeneratorPtr bnds = pq.mstructure->createBondGenerator(); pq.configureBondGenerator(*bnds); @@ -155,6 +157,8 @@ void PQEvaluatorOptimized::updateValue(PairQuantity& pq, } // Remove contributions from the extra sites in the old structure assert(sd.stru0 == mlast_structure); + // The value belongs to no structure until this update completes. + mlast_structure.reset(); int cntsites0 = sd.stru0->countSites(); BaseBondGeneratorPtr bnds0 = sd.stru0->createBondGenerator(); pq.configureBondGenerator(*bnds0); diff --git a/src/diffpy/srreal/PeriodicStructureAdapter.cpp b/src/diffpy/srreal/PeriodicStructureAdapter.cpp index 280b41d..710a3b9 100644 --- a/src/diffpy/srreal/PeriodicStructureAdapter.cpp +++ b/src/diffpy/srreal/PeriodicStructureAdapter.cpp @@ -19,12 +19,16 @@ * *****************************************************************************/ +#include +#include #include +#include #include #include #include #include +#include using namespace std; @@ -106,90 +110,455 @@ bool operator!=(const PeriodicStructureAdapter& stru0, // class PeriodicStructureBondGenerator ////////////////////////////////////////////////////////////////////////////// +// Local Constants ----------------------------------------------------------- + +namespace { + +/// Fill candidate buffers with at least this many bonds when possible. +const int CANDIDATE_BATCH = 256; + +/// Candidate margin relative to the coordinate scale, well above round-off. +const double CANDIDATE_MARGIN = 1e-9; + +/// Major points per kernel call in searchAlongMajor. +const int MAJOR_BLOCK = 512; + +/// The kernel runs along the major points when there are fewer minor points +/// than this. Their candidate flags must fit in mblockbits. +const int MINOR_FEW = 64; + +/// Lattice vectors in one chunk of mlattice, which takes 0.8 MB. +const int LATTICE_CHUNK = 1 << 15; + +/// Squared limits for BondCandidatesArgs, negative means no bond or limit. +double lowerLimitSquared(double rlo) { + return (rlo > 0.0) ? (rlo * rlo) : -1.0; +} +double upperLimitSquared(double rhi) { + return (rhi < 0.0) ? -1.0 : (rhi * rhi); +} + +/// Larger of d and dmax, where NaN wins. +double maxOrNaN(double d, double dmax) { + return (std::isnan(d) || d > dmax) ? d : dmax; +} + +/// Narrow the integer range [lo, hi) to [dlo, dhi). NaN keeps the bound. +void narrowRange(int& lo, int& hi, double dlo, double dhi) { + if (dlo > lo) lo = (dlo < hi) ? int(dlo) : hi; + if (dhi < hi) hi = (dhi > lo) ? int(dhi) : lo; +} + +} // namespace + // Constructor --------------------------------------------------------------- PeriodicStructureBondGenerator::PeriodicStructureBondGenerator( StructureAdapterConstPtr adpt) - : BaseBondGenerator(adpt) { - mpstructure = dynamic_cast(adpt.get()); - assert(mpstructure); - int cntsites = mpstructure->countSites(); - mcartesian_positions_uc.reserve(cntsites); + : PeriodicStructureBondGenerator(adpt, BondOrder::LatticeFirst) { const Lattice& L = mpstructure->getLattice(); - PeriodicStructureAdapter::const_iterator ai = mpstructure->begin(); - R3::Vector xyzc; - for (; ai != mpstructure->end(); ++ai) { - xyzc = L.ucvCartesian(ai->xyz_cartn); - mcartesian_positions_uc.push_back(xyzc); + for (int i = 0; i < mpstructure->countSites(); ++i) { + this->addPosition(i, L.ucvCartesian((*mpstructure)[i].xyz_cartn)); } } +PeriodicStructureBondGenerator::PeriodicStructureBondGenerator( + StructureAdapterConstPtr adpt, BondOrder order) + : BaseBondGenerator(adpt), + mpstructure(dynamic_cast(adpt.get())), + mbondorder(order) { + assert(mpstructure); +} + // Public Methods ------------------------------------------------------------ void PeriodicStructureBondGenerator::rewind() { // Delay msphere instantiation to here instead of in constructor, // so it is possible to use setRmin, setRmax. - if (!msphere.get()) { - const Lattice& L = mpstructure->getLattice(); - double buffzone = L.ucMaxDiagonalLength(); - double rsphmin = this->getRmin() - buffzone; - double rsphmax = this->getRmax() + buffzone; - msphere.reset(new PointsInSphere(rsphmin, rsphmax, L)); - } - // BaseBondGenerator::rewind calls this->rewindSymmetry, - // which takes care of msphere configuration + if (!msphere_ready) this->buildSphere(); + // BaseBondGenerator::rewind calls this->rewindSymmetry this->BaseBondGenerator::rewind(); } void PeriodicStructureBondGenerator::selectAnchorSite(int anchor) { this->BaseBondGenerator::selectAnchorSite(anchor); - mr0 = mcartesian_positions_uc[anchor]; + const int i = msitebegin[anchor]; + mr0 = R3::Vector(mcell.x[i], mcell.y[i], mcell.z[i]); } void PeriodicStructureBondGenerator::setRmin(double rmin) { // destroy msphere so it will be created on rewind with new rmin - if (this->getRmin() != rmin) msphere.reset(); + if (this->getRmin() != rmin) this->releaseSphere(); this->BaseBondGenerator::setRmin(rmin); } void PeriodicStructureBondGenerator::setRmax(double rmax) { // destroy msphere so it will be created on rewind with new rmax - if (this->getRmax() != rmax) msphere.reset(); + if (this->getRmax() != rmax) this->releaseSphere(); this->BaseBondGenerator::setRmax(rmax); } // Protected Methods --------------------------------------------------------- -bool PeriodicStructureBondGenerator::iterateSymmetry() { - msphere->next(); - bool done = msphere->finished(); - mrcsphere = - done ? R3::zerovector : mpstructure->getLattice().cartesian(msphere->mno()); - return !done; -} +// site1 is a lattice vector plus a position of a selected site. Bonds are +// ordered by major then minor points as in the original nested loops, where +// major are the lattice vectors here and the positions in +// CrystalStructureBondGenerator. The kernel selects candidates and +// BaseBondGenerator checks their exact distance. +// +// Lattice vectors near the anchor come in chunks of LATTICE_CHUNK points. +// When positions are major, each position goes through all chunks. void PeriodicStructureBondGenerator::rewindSymmetry() { - msphere->rewind(); - mrcsphere = msphere->finished() - ? R3::zerovector - : mpstructure->getLattice().cartesian(msphere->mno()); - this->updater1(); + this->selectPositions(); + this->selectAnchorRows(); + this->loadLatticeChunk(true); + mnextmajor = 0; + mchunksearched = false; + mcandcount = 0; + mcandidx = 0; + this->seekCandidate(); } void PeriodicStructureBondGenerator::getNextBond() { - ++msite_current; - // go back to the first site if there is next symmetry element - if (msite_current >= msite_last && this->iterateSymmetry()) { - msite_current = msite_first; + ++mcandidx; + this->seekCandidate(); +} + +void PeriodicStructureBondGenerator::addPosition(int site, + const R3::Vector& xyz) { + assert(site + 1 == int(msitebegin.size()) || + site + 2 == int(msitebegin.size())); + msitebegin.resize(site + 2, mcell.size()); + mcell.push_back(xyz); + mcellsite.push_back(site); + msitebegin.back() = mcell.size(); + msphere_ready = false; +} + +double PeriodicStructureBondGenerator::PointArray::maxDistance( + const R3::Vector& c) const { + double d2max = 0.0; + for (int i = 0; i < this->size(); ++i) { + const double dx = x[i] - c[0], dy = y[i] - c[1], dz = z[i] - c[2]; + d2max = maxOrNaN(dx * dx + dy * dy + dz * dz, d2max); } - // update values only if not finished - if (!this->finished()) this->updater1(); + return sqrt(d2max); } // Private Methods ----------------------------------------------------------- -void PeriodicStructureBondGenerator::updater1() { - mr1 = mrcsphere + mcartesian_positions_uc[this->site1()]; +void PeriodicStructureBondGenerator::buildSphere() { + const Lattice& L = mpstructure->getLattice(); + const double buffzone = L.ucMaxDiagonalLength(); + PointsInSphere sph(this->getRmin() - buffzone, this->getRmax() + buffzone, L); + // Join consecutive points along o into rows and store the points of + // small spheres. L.cartesian(sph.mno()) gives the bonds exactly the + // coordinates of the original PointsInSphere loops. + msphere.clear(); + mspherepoints.clear(); + int npoints = 0; + for (sph.rewind(); !sph.finished(); sph.next(), ++npoints) { + if (npoints < LATTICE_CHUNK) { + mspherepoints.push_back(L.cartesian(sph.mno())); + } + const int m = sph.m(), n = sph.n(), o = sph.o(); + if (!msphere.empty()) { + LatticeRow& last = msphere.back(); + if (last.m == m && last.n == n && last.ohi == o) { + ++last.ohi; + continue; + } + } + msphere.push_back({m, n, o, o + 1, npoints}); + } + if (npoints > LATTICE_CHUNK) { + mspherepoints = {}; + } + // The norm along a row is largest at one of its ends. + double maxlatnorm = 0.0; + for (const LatticeRow& rw : msphere) { + for (int o : {rw.olo, rw.ohi - 1}) { + const double r = R3::norm(L.cartesian(R3::Vector(rw.m, rw.n, o))); + maxlatnorm = maxOrNaN(r, maxlatnorm); + } + } + // sphere around the unit cell positions + mcellcenter = R3::zerovector; + const int n = mcell.size(); + for (int i = 0; i < n; ++i) { + mcellcenter += R3::Vector(mcell.x[i], mcell.y[i], mcell.z[i]); + } + if (n) mcellcenter /= n; + mcellradius = mcell.maxDistance(mcellcenter); + // Widen the r-range by the rounding error of the kernel distances. + const double maxnorm = mcell.maxDistance() + maxlatnorm; + const double scale = + 2 * maxnorm + fabs(this->getRmin()) + fabs(this->getRmax()) + 1.0; + mroundmargin = CANDIDATE_MARGIN * scale; + mcandd2lo = lowerLimitSquared(this->getRmin() - mroundmargin); + mcandd2hi = upperLimitSquared(this->getRmax() + mroundmargin); + msphere_ready = true; +} + +void PeriodicStructureBondGenerator::releaseSphere() { + msphere = {}; + mspherepoints = {}; + msphere_ready = false; +} + +void PeriodicStructureBondGenerator::selectPositions() { + const int nsel = msite_last - msite_first; + assert(nsel > 0); + const int site0 = *msite_first; + // consecutive sites have their positions together in mcell + if (msite_consecutive) { + mpositionsbegin = msitebegin[site0]; + mpositions = PointSpan(mcell, mpositionsbegin, msitebegin[site0 + nsel]); + return; + } + mgathered.clear(); + mgatheredcell.clear(); + mgatheredsite.clear(); + for (int k = 0; k < nsel; ++k) { + const int site = msite_first[k]; + for (int i = msitebegin[site]; i < msitebegin[site + 1]; ++i) { + mgathered.push_back(mcell.x[i], mcell.y[i], mcell.z[i]); + mgatheredcell.push_back(i); + mgatheredsite.push_back(k); + } + } + mpositionsbegin = -1; + mpositions = PointSpan(mgathered); +} + +void PeriodicStructureBondGenerator::selectAnchorRows() { + manchorrows.clear(); + // The sphere around the unit cell positions shifted by the lattice vector + // must reach the r-range of the anchor. + const double margin = mcellradius + mroundmargin; + const double rlo = this->getRmin() - margin; + const double rhi = this->getRmax() + margin; + if (rhi < 0.0) return; + // Row points are v + o * c relative to the anchor, so their squared + // distance is dmin2 + cc * (o - os)^2. The o-limits get one more point + // on each side to cover their rounding. NaN keeps the whole row. + const Lattice& L = mpstructure->getLattice(); + const R3::Vector c = L.cartesian(R3::Vector(0, 0, 1)); + const double cc = R3::dot(c, c); + const R3::Vector center = mcellcenter - mr0; + for (const LatticeRow& rw : msphere) { + const R3::Vector v = L.cartesian(R3::Vector(rw.m, rw.n, 0)) + center; + const double vc = R3::dot(v, c); + const double os = -vc / cc; + double dmin2 = R3::dot(v, v) - vc * vc / cc; + if (dmin2 < 0.0) dmin2 = 0.0; + if (dmin2 > rhi * rhi) continue; + // points within rhi + const double w = sqrt((rhi * rhi - dmin2) / cc); + int olo = rw.olo; + int ohi = rw.ohi; + narrowRange(olo, ohi, ceil(os - w) - 1, floor(os + w) + 2); + if (olo >= ohi) continue; + // points well within rlo, they split the row + if (rlo > 0.0 && dmin2 < rlo * rlo) { + const double u = sqrt((rlo * rlo - dmin2) / cc); + int holelo = olo; + int holehi = ohi; + narrowRange(holelo, holehi, floor(os - u + 1) + 1, ceil(os + u - 1)); + if (holelo < holehi) { + if (olo < holelo) addAnchorRow(rw, olo, holelo); + if (holehi < ohi) addAnchorRow(rw, holehi, ohi); + continue; + } + } + addAnchorRow(rw, olo, ohi); + } +} + +void PeriodicStructureBondGenerator::addAnchorRow(const LatticeRow& rw, int olo, + int ohi) { + manchorrows.push_back({rw.m, rw.n, olo, ohi, rw.point + (olo - rw.olo)}); +} + +void PeriodicStructureBondGenerator::loadLatticeChunk(bool first) { + LatticeChunk& ch = mchunk; + ch.points.clear(); + ch.first = first; + if (first) { + ch.nextrow = 0; + ch.nextpoint = 0; + } + const int nrows = manchorrows.size(); + if (!mspherepoints.x.empty()) { + // copy the stored sphere points, which are all in one chunk + for (; ch.nextrow < nrows; ++ch.nextrow) { + const LatticeRow& rw = manchorrows[ch.nextrow]; + ch.points.append(mspherepoints, rw.point, rw.point + (rw.ohi - rw.olo)); + } + } else { + // The retained rows already contain the exact integer coordinates. + const Lattice& L = mpstructure->getLattice(); + while (ch.nextrow < nrows && ch.points.size() < LATTICE_CHUNK) { + const LatticeRow& rw = manchorrows[ch.nextrow]; + const int count = + min(rw.ohi - rw.olo - ch.nextpoint, LATTICE_CHUNK - ch.points.size()); + for (int i = 0; i < count; ++i) { + ch.points.push_back( + L.cartesian(R3::Vector(rw.m, rw.n, rw.olo + ch.nextpoint + i))); + } + ch.nextpoint += count; + if (ch.nextpoint == rw.ohi - rw.olo) { + ++ch.nextrow; + ch.nextpoint = 0; + } + } + } + ch.last = (ch.nextrow >= nrows); +} + +bool PeriodicStructureBondGenerator::refillCandidates() { + mcandcount = 0; + mcandidx = 0; + const PointSpan lattice(mchunk.points); + if (mbondorder == BondOrder::LatticeFirst) { + if (mnextmajor >= lattice.n) { + if (mchunk.last) return false; + this->loadLatticeChunk(false); + mnextmajor = 0; + return true; + } + this->searchMajorPoints(lattice, mpositions); + return true; + } + if (mnextmajor >= mpositions.n) return false; + // All positions see the same lattice vectors when they fit in one chunk. + if (mchunk.first && mchunk.last) { + this->searchMajorPoints(mpositions, lattice); + return true; + } + // Otherwise every position goes through all chunks. + if (mchunksearched) { + if (!mchunk.last) { + this->loadLatticeChunk(false); + } else { + if (++mnextmajor >= mpositions.n) return false; + this->loadLatticeChunk(true); + } + mchunksearched = false; + return true; + } + this->growCandidates(lattice.n); + mcandcount = + this->searchMajorPoint(mnextmajor, mpositions, lattice, mcandminor.data()); + fill_n(mcandmajor.begin(), mcandcount, mnextmajor); + mchunksearched = true; + return true; +} + +void PeriodicStructureBondGenerator::growCandidates(size_t capacity) { + if (mcandmajor.size() < capacity) { + mcandmajor.resize(capacity); + mcandminor.resize(capacity); + } +} + +int PeriodicStructureBondGenerator::searchMajorPoint(int imajor, + const PointSpan& major, + const PointSpan& minor, + int* idx) const { + return kernels::activeCalculationKernel().bondCandidates({ + .x = minor.x, + .y = minor.y, + .z = minor.z, + .n = minor.n, + .ox = major.x[imajor] - mr0[0], + .oy = major.y[imajor] - mr0[1], + .oz = major.z[imajor] - mr0[2], + .d2lo = mcandd2lo, + .d2hi = mcandd2hi, + .idx = idx, + }); +} + +void PeriodicStructureBondGenerator::searchMajorPoints(const PointSpan& major, + const PointSpan& minor) { + if (minor.n == 0) { + mnextmajor = major.n; + return; + } + // The kernel is faster along the major points when the minor ones are few. + if (minor.n < MINOR_FEW && minor.n < major.n) { + this->searchAlongMajor(major, minor); + return; + } + this->growCandidates(CANDIDATE_BATCH + size_t(minor.n)); + for (; mnextmajor < major.n && mcandcount < CANDIDATE_BATCH; ++mnextmajor) { + int* idx = mcandminor.data() + mcandcount; + const int m = this->searchMajorPoint(mnextmajor, major, minor, idx); + fill_n(mcandmajor.begin() + mcandcount, m, mnextmajor); + mcandcount += m; + } +} + +void PeriodicStructureBondGenerator::searchAlongMajor(const PointSpan& major, + const PointSpan& minor) { + this->growCandidates(CANDIDATE_BATCH + size_t(minor.n) * MAJOR_BLOCK); + // Bit j of mblockbits[i] flags the candidate bond of the major point i + // in the block and the minor point j. + mblockidx.resize(MAJOR_BLOCK); + mblockbits.resize(MAJOR_BLOCK); + const auto bondCandidates = kernels::activeCalculationKernel().bondCandidates; + kernels::BondCandidatesArgs args{.d2lo = mcandd2lo, .d2hi = mcandd2hi}; + args.idx = mblockidx.data(); + while (mnextmajor < major.n && mcandcount < CANDIDATE_BATCH) { + const int nblock = min(MAJOR_BLOCK, major.n - mnextmajor); + fill_n(mblockbits.begin(), nblock, 0); + args.x = major.x + mnextmajor; + args.y = major.y + mnextmajor; + args.z = major.z + mnextmajor; + args.n = nblock; + for (int j = 0; j < minor.n; ++j) { + args.ox = minor.x[j] - mr0[0]; + args.oy = minor.y[j] - mr0[1]; + args.oz = minor.z[j] - mr0[2]; + const int m = bondCandidates(args); + for (int k = 0; k < m; ++k) mblockbits[mblockidx[k]] |= uint64_t(1) << j; + } + for (int i = 0; i < nblock; ++i) { + for (uint64_t bits = mblockbits[i]; bits; bits &= bits - 1) { + mcandmajor[mcandcount] = mnextmajor + i; + mcandminor[mcandcount] = countr_zero(bits); + ++mcandcount; + } + } + mnextmajor += nblock; + } +} + +void PeriodicStructureBondGenerator::seekCandidate() { + while (mcandidx >= mcandcount) { + if (!this->refillCandidates()) { + msite_current = msite_last; + return; + } + } + const int imajor = mcandmajor[mcandidx]; + const int iminor = mcandminor[mcandidx]; + const bool latticefirst = (mbondorder == BondOrder::LatticeFirst); + const int ipos = latticefirst ? iminor : imajor; + const int ilat = latticefirst ? imajor : iminor; + if (mpositionsbegin >= 0) { + mcellindex = mpositionsbegin + ipos; + msite_current = msite_first + (mcellsite[mcellindex] - *msite_first); + } else { + mcellindex = mgatheredcell[ipos]; + msite_current = msite_first + mgatheredsite[ipos]; + } + const PointArray& lattice = mchunk.points; + mr1[0] = lattice.x[ilat] + mpositions.x[ipos]; + mr1[1] = lattice.y[ilat] + mpositions.y[ipos]; + mr1[2] = lattice.z[ilat] + mpositions.z[ipos]; this->updateDistance(); } diff --git a/src/diffpy/srreal/PeriodicStructureAdapter.hpp b/src/diffpy/srreal/PeriodicStructureAdapter.hpp index 3b2cb43..6e3bcff 100644 --- a/src/diffpy/srreal/PeriodicStructureAdapter.hpp +++ b/src/diffpy/srreal/PeriodicStructureAdapter.hpp @@ -22,6 +22,11 @@ #ifndef PERIODICSTRUCTUREADAPTER_HPP_INCLUDED #define PERIODICSTRUCTUREADAPTER_HPP_INCLUDED +#include +#include +#include +#include + #include #include #include @@ -74,6 +79,11 @@ class DLL_EXPORT PeriodicStructureBondGenerator : public BaseBondGenerator { public: // constructors PeriodicStructureBondGenerator(StructureAdapterConstPtr); + // mpositions points into the generator's own arrays + PeriodicStructureBondGenerator(const PeriodicStructureBondGenerator&) = + delete; + PeriodicStructureBondGenerator& operator=( + const PeriodicStructureBondGenerator&) = delete; // methods // loop control @@ -85,20 +95,149 @@ class DLL_EXPORT PeriodicStructureBondGenerator : public BaseBondGenerator { virtual void setRmax(double); protected: + /// Outer loop of the bond order, the lattice vectors or the positions of + /// the selected sites. + enum class BondOrder { LatticeFirst, PositionFirst }; + + /// Cartesian points stored by coordinate for the calculation kernels. + struct PointArray { + std::vector x; + std::vector y; + std::vector z; + + int size() const { return int(x.size()); } + void clear() { + x.clear(); + y.clear(); + z.clear(); + } + void push_back(double px, double py, double pz) { + x.push_back(px); + y.push_back(py); + z.push_back(pz); + } + void push_back(const R3::Vector& v) { this->push_back(v[0], v[1], v[2]); } + /// Append the points [first, last) of a. + void append(const PointArray& a, int first, int last) { + x.insert(x.end(), a.x.begin() + first, a.x.begin() + last); + y.insert(y.end(), a.y.begin() + first, a.y.begin() + last); + z.insert(z.end(), a.z.begin() + first, a.z.begin() + last); + } + /// Largest distance of the points from c. + double maxDistance(const R3::Vector& c = R3::zerovector) const; + }; + + /// Range of points in a PointArray. + struct PointSpan { + const double* x = nullptr; + const double* y = nullptr; + const double* z = nullptr; + int n = 0; + + PointSpan() = default; + PointSpan(const PointArray& a, int first, int last) + : x(a.x.data() + first), + y(a.y.data() + first), + z(a.z.data() + first), + n(last - first) {} + PointSpan(const PointArray& a) : PointSpan(a, 0, a.size()) {} + }; + + /// Lattice vectors [m, n, o] for o in [olo, ohi). point is the index + /// of [m, n, olo] in the sequence of all sphere points. + struct LatticeRow { + int m; + int n; + int olo; + int ohi; + int point; + }; + + /// One chunk of the lattice vectors near the anchor. + struct LatticeChunk { + PointArray points; + /// The chunk is the first or the last one. + bool first = true; + bool last = true; + /// The next chunk starts at this anchor row and offset within it. + int nextrow = 0; + int nextpoint = 0; + }; + + /// Create generator without site1 positions, add them with addPosition. + PeriodicStructureBondGenerator(StructureAdapterConstPtr, BondOrder); + // data const PeriodicStructureAdapter* mpstructure; - std::unique_ptr msphere; - R3::Vector mrcsphere; + /// Index in the site1 positions of the current bond. + int mcellindex = 0; // methods - virtual bool iterateSymmetry(); virtual void rewindSymmetry(); virtual void getNextBond(); - virtual void updater1(); + /// Add unit cell position of site1. Call in the order of the sites. + void addPosition(int site, const R3::Vector& xyz); private: // data - std::vector mcartesian_positions_uc; + BondOrder mbondorder; + /// Unit cell positions of site1 for all sites, ordered by site. + PointArray mcell; + /// Site of each position in mcell. + std::vector mcellsite; + /// Positions of site i are mcell points [msitebegin[i], msitebegin[i + 1]). + std::vector msitebegin = {0}; + /// Rows of lattice vectors within the r-range plus a unit cell buffer. + std::vector msphere; + /// All msphere lattice vectors when they fit in one chunk, else empty. + PointArray mspherepoints; + bool msphere_ready = false; + /// Center and radius of a sphere around the mcell points. + R3::Vector mcellcenter; + double mcellradius = 0.0; + /// Widening of the r-range for the rounding error of kernel distances. + double mroundmargin = 0.0; + /// Positions of the selected sites. + PointSpan mpositions; + /// Index of mpositions[0] in mcell, or -1 when they are in mgathered. + int mpositionsbegin = 0; + /// Selected positions gathered when the sites are not consecutive. + PointArray mgathered; + std::vector mgatheredcell; + std::vector mgatheredsite; + /// Parts of the msphere rows that can make a bond with the anchor. + std::vector manchorrows; + LatticeChunk mchunk; + /// Squared distance limits of the candidate bonds. + double mcandd2lo = -1.0; + double mcandd2hi = -1.0; + /// Candidate bonds as indices of the major (outer loop) and minor points. + std::vector mcandmajor; + std::vector mcandminor; + int mcandcount = 0; + int mcandidx = 0; + int mnextmajor = 0; + /// True when the lattice chunk was searched for the major point + /// mnextmajor, used when positions are major and there are more chunks. + bool mchunksearched = false; + /// Kernel output and candidate bits for one block of major points. + std::vector mblockidx; + std::vector mblockbits; + + // methods + void buildSphere(); + void releaseSphere(); + void selectPositions(); + void selectAnchorRows(); + void addAnchorRow(const LatticeRow& rw, int olo, int ohi); + void loadLatticeChunk(bool first); + bool refillCandidates(); + void growCandidates(size_t capacity); + int searchMajorPoint(int imajor, const PointSpan& major, + const PointSpan& minor, int* idx) const; + void searchMajorPoints(const PointSpan& major, const PointSpan& minor); + void searchAlongMajor(const PointSpan& major, const PointSpan& minor); + void seekCandidate(); }; } // namespace srreal diff --git a/src/diffpy/srreal/StructureAdapter.cpp b/src/diffpy/srreal/StructureAdapter.cpp index ad1e541..295ca39 100644 --- a/src/diffpy/srreal/StructureAdapter.cpp +++ b/src/diffpy/srreal/StructureAdapter.cpp @@ -63,24 +63,22 @@ StructureDifference StructureAdapter::diff( double meanSquareDisplacement(const R3::Matrix& Uijcartn, const R3::Vector& s, bool anisotropy) { - double rv; - if (anisotropy) { - assert(R3::norm(s) > 0); - assert(eps_eq(Uijcartn(0, 1), Uijcartn(1, 0))); - assert(eps_eq(Uijcartn(0, 2), Uijcartn(2, 0))); - assert(eps_eq(Uijcartn(1, 2), Uijcartn(2, 1))); - static R3::Vector sn; - sn = s / R3::norm(s); - rv = Uijcartn(0, 0) * sn(0) * sn(0) + Uijcartn(1, 1) * sn(1) * sn(1) + - Uijcartn(2, 2) * sn(2) * sn(2) + 2 * Uijcartn(0, 1) * sn(0) * sn(1) + - 2 * Uijcartn(0, 2) * sn(0) * sn(2) + - 2 * Uijcartn(1, 2) * sn(1) * sn(2); - } else { - assert(eps_eq(Uijcartn(0, 0), Uijcartn(1, 1))); - assert(eps_eq(Uijcartn(0, 0), Uijcartn(2, 2))); - rv = Uijcartn(0, 0); + const R3::Matrix& U = Uijcartn; + if (!anisotropy) { + assert(eps_eq(U(0, 0), U(1, 1))); + assert(eps_eq(U(0, 0), U(2, 2))); + return U(0, 0); } - return rv; + assert(R3::norm(s) > 0); + assert(eps_eq(U(0, 1), U(1, 0))); + assert(eps_eq(U(0, 2), U(2, 0))); + assert(eps_eq(U(1, 2), U(2, 1))); + // s'Us / s's is the MSD along s without normalizing s first. + const double s0 = s[0], s1 = s[1], s2 = s[2]; + const double sUs = + U(0, 0) * s0 * s0 + U(1, 1) * s1 * s1 + U(2, 2) * s2 * s2 + + 2 * (U(0, 1) * s0 * s1 + U(0, 2) * s0 * s2 + U(1, 2) * s1 * s2); + return sUs / (s0 * s0 + s1 * s1 + s2 * s2); } double maxUii(StructureAdapterPtr stru) { diff --git a/src/tests/TestBondGeneratorBatching.cpp b/src/tests/TestBondGeneratorBatching.cpp new file mode 100644 index 0000000..33bb24f --- /dev/null +++ b/src/tests/TestBondGeneratorBatching.cpp @@ -0,0 +1,575 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * Compare the periodic and crystal bond generators with a reference loop. + * + *****************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "test_helpers.hpp" + +using namespace diffpy::srreal; + +namespace { + +struct Bond { + int site0; + int site1; + R3::Vector r1; + R3::Vector r01; + double distance; + R3::Matrix u1; +}; + +struct Config { + double rmin; + double rmax; +}; + +constexpr std::array CONFIGS = { + {{0.0, 9.7}, {5.0, 4.0}, {0.0, 0.0}, {6.1, 6.9}}}; + +/// Hexagonal, orthorhombic and triclinic cells. +const char* const STRUCTURES[] = {"ZnS_wurtzite.stru", "CaTiO3.stru", + "alpha_K2Bi8Se13.stru"}; + +/// Site selections exercised for every anchor. +std::vector siteSelections(int cntsites, int anchor) { + std::vector rv(4); + for (int i = 0; i <= anchor; ++i) rv[0].push_back(i); + for (int i = cntsites - 1; i >= 0; i -= 2) rv[1].push_back(i); + for (int i = 1; i < cntsites - 1; ++i) rv[2].push_back(i); + // rv[3] stays empty + return rv; +} + +bool inRange(const Bond& b, const Config& cfg) { + using diffpy::mathutils::eps_eq; + const bool outofrange = (b.distance < cfg.rmin) || (b.distance > cfg.rmax); + return !outofrange && !eps_eq(b.distance, 0.0); +} + +/// Reference: msphere points outer loop, selected sites inner loop. +std::vector referencePeriodic(const PeriodicStructureAdapter& stru, + int anchor, const SiteIndices& sel, + const Config& cfg) { + const Lattice& L = stru.getLattice(); + const double buffzone = L.ucMaxDiagonalLength(); + PointsInSphere sph(cfg.rmin - buffzone, cfg.rmax + buffzone, L); + const R3::Vector r0 = L.ucvCartesian(stru[anchor].xyz_cartn); + std::vector rv; + for (sph.rewind(); !sph.finished(); sph.next()) { + const R3::Vector rc = L.cartesian(sph.mno()); + for (int j : sel) { + Bond b; + b.site0 = anchor; + b.site1 = j; + const R3::Vector xyz = L.ucvCartesian(stru[j].xyz_cartn); + for (int k = 0; k < 3; ++k) b.r1[k] = rc[k] + xyz[k]; + for (int k = 0; k < 3; ++k) b.r01[k] = b.r1[k] - r0[k]; + b.distance = R3::norm(b.r01); + b.u1 = stru.siteCartesianUij(j); + if (inRange(b, cfg)) rv.push_back(b); + } + } + return rv; +} + +/// Reference: sites, then symmetry positions, then msphere points. +std::vector referenceCrystal(const CrystalStructureAdapter& stru, + int anchor, const SiteIndices& sel, + const Config& cfg) { + const Lattice& L = stru.getLattice(); + const double buffzone = L.ucMaxDiagonalLength(); + PointsInSphere sph(cfg.rmin - buffzone, cfg.rmax + buffzone, L); + const R3::Vector r0 = stru.getEquivalentAtoms(anchor)[0].xyz_cartn; + std::vector rv; + for (int j : sel) { + for (const Atom& a : stru.getEquivalentAtoms(j)) { + for (sph.rewind(); !sph.finished(); sph.next()) { + const R3::Vector rc = L.cartesian(sph.mno()); + Bond b; + b.site0 = anchor; + b.site1 = j; + for (int k = 0; k < 3; ++k) b.r1[k] = rc[k] + a.xyz_cartn[k]; + for (int k = 0; k < 3; ++k) b.r01[k] = b.r1[k] - r0[k]; + b.distance = R3::norm(b.r01); + b.u1 = a.uij_cartn; + if (inRange(b, cfg)) rv.push_back(b); + } + } + } + return rv; +} + +std::vector generated(BaseBondGenerator& bnds) { + std::vector rv; + for (bnds.rewind(); !bnds.finished(); bnds.next()) { + Bond b; + b.site0 = bnds.site0(); + b.site1 = bnds.site1(); + b.r1 = bnds.r1(); + b.r01 = bnds.r01(); + b.distance = bnds.distance(); + b.u1 = bnds.Ucartesian1(); + rv.push_back(b); + } + return rv; +} + +bool sameBond(const Bond& a, const Bond& e) { + bool rv = (a.site0 == e.site0) && (a.site1 == e.site1) && + (std::fabs(a.distance - e.distance) <= 1e-14 * e.distance); + for (int k = 0; k < 3; ++k) { + rv = rv && (a.r1[k] == e.r1[k]) && (a.r01[k] == e.r01[k]); + for (int l = 0; l < 3; ++l) rv = rv && (a.u1(k, l) == e.u1(k, l)); + } + return rv; +} + +void checkSameBonds(const std::vector& actual, + const std::vector& expected) { + REQUIRE(actual.size() == expected.size()); + size_t i = 0; + while (i < actual.size() && sameBond(actual[i], expected[i])) ++i; + if (i == actual.size()) return; + INFO("first different bond " << i); + const Bond& a = actual[i]; + const Bond& e = expected[i]; + CHECK(a.site0 == e.site0); + CHECK(a.site1 == e.site1); + for (int k = 0; k < 3; ++k) { + CHECK(a.r1[k] == e.r1[k]); + CHECK(a.r01[k] == e.r01[k]); + } + CHECK(std::fabs(a.distance - e.distance) <= 1e-14 * e.distance); + for (int k = 0; k < 3; ++k) { + for (int l = 0; l < 3; ++l) CHECK(a.u1(k, l) == e.u1(k, l)); + } + FAIL("bond sequences differ"); +} + +/// Crystal structure from a PDFfit file with a few symmetry operations. +CrystalStructureAdapterPtr loadTestCrystalStructure(const std::string& f) { + StructureAdapterPtr p = loadTestPeriodicStructure(f); + const PeriodicStructureAdapter& ps = + dynamic_cast(*p); + const Lattice& L = ps.getLattice(); + CrystalStructureAdapterPtr rv(new CrystalStructureAdapter); + rv->setLatPar(L.a(), L.b(), L.c(), L.alpha(), L.beta(), L.gamma()); + for (int i = 0; i < ps.countSites(); ++i) rv->append(ps[i]); + R3::Vector t = R3::zerovector; + rv->addSymOp(R3::identity(), t); + rv->addSymOp(-1.0 * R3::identity(), t); + R3::Matrix R = R3::identity(); + R(0, 0) = -1.0; + R(1, 1) = -1.0; + t[0] = 0.5; + t[1] = 0.5; + rv->addSymOp(R, t); + return rv; +} + +/// Select bonds up to 6 from the site 1, used with CaTiO3. +void nearSite1(BaseBondGenerator& bnds) { + bnds.setRmax(6.0); + bnds.selectAnchorSite(1); +} + +/// New bond generator for the bonds up to 6 from the site 1. +BaseBondGeneratorPtr nearSite1(const StructureAdapter& stru) { + BaseBondGeneratorPtr rv = stru.createBondGenerator(); + nearSite1(*rv); + return rv; +} + +template +void compareAllConfigurations(const Tstru& stru, Treference reference) { + BaseBondGeneratorPtr bnds = stru.createBondGenerator(); + const int n = stru.countSites(); + for (const Config& cfg : CONFIGS) { + bnds->setRmin(cfg.rmin); + bnds->setRmax(cfg.rmax); + for (int anchor = 0; anchor < n; ++anchor) { + bnds->selectAnchorSite(anchor); + const std::vector sels = siteSelections(n, anchor); + for (size_t isel = 0; isel < sels.size(); ++isel) { + INFO("rmin " << cfg.rmin << ", rmax " << cfg.rmax << ", anchor " + << anchor << ", selection " << isel); + if (isel == 0) { + bnds->selectSiteRange(0, anchor + 1); + } else { + bnds->selectSites(sels[isel]); + } + checkSameBonds(generated(*bnds), + reference(stru, anchor, sels[isel], cfg)); + } + } + } +} + +} // namespace + +TEST_CASE("TestBondGeneratorBatching::test_periodic", + "[BondGeneratorBatching]") { + testAllKernels([] { + for (const char* f : STRUCTURES) { + INFO("structure " << f); + StructureAdapterPtr p = loadTestPeriodicStructure(f); + const PeriodicStructureAdapter& stru = + dynamic_cast(*p); + compareAllConfigurations(stru, referencePeriodic); + } + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_crystal", + "[BondGeneratorBatching]") { + testAllKernels([] { + for (const char* f : STRUCTURES) { + INFO("structure " << f); + CrystalStructureAdapterPtr stru = loadTestCrystalStructure(f); + // the symmetry operations must give extra equivalent positions + int cnteq = 0; + for (int i = 0; i < stru->countSites(); ++i) { + cnteq += stru->getEquivalentAtoms(i).size(); + } + CHECK(cnteq > stru->countSites()); + compareAllConfigurations(*stru, referenceCrystal); + } + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_crystal_changes", + "[BondGeneratorBatching]") { + // A new bond generator must follow every change, also in place. + testAllKernels([] { + CrystalStructureAdapterPtr stru = loadTestCrystalStructure("CaTiO3.stru"); + auto bonds = [](const CrystalStructureAdapter& s) { + BaseBondGeneratorPtr bnds = s.createBondGenerator(); + bnds->setRmax(6.0); + bnds->selectAnchorSite(0); + return generated(*bnds); + }; + // the same structure rebuilt without cached symmetry positions + auto rebuilt = [](const CrystalStructureAdapter& s) { + const Lattice& L = s.getLattice(); + CrystalStructureAdapterPtr rv(new CrystalStructureAdapter); + rv->setLatPar(L.a(), L.b(), L.c(), L.alpha(), L.beta(), L.gamma()); + rv->assign(s.begin(), s.end()); + rv->setSymmetryPrecision(s.getSymmetryPrecision()); + for (int i = 0; i < s.countSymOps(); ++i) rv->addSymOp(s.getSymOp(i)); + return rv; + }; + auto same = [](const std::vector& a, const std::vector& b) { + return a.size() == b.size() && + std::equal(a.begin(), a.end(), b.begin(), sameBond); + }; + std::vector last = bonds(*stru); + REQUIRE(!last.empty()); + checkSameBonds(last, bonds(*rebuilt(*stru))); + R3::Matrix R = R3::identity(); + R(2, 2) = -1.0; + const std::pair> changes[] = { + {"atom position", [&] { stru->at(1).xyz_cartn[0] += 0.3; }}, + {"atom displacement", [&] { stru->at(2).uij_cartn(0, 0) += 0.01; }}, + {"lattice", [&] { stru->setLatPar(5.5, 7.7, 5.4, 90.0, 95.0, 90.0); }}, + {"new atom", [&] { stru->append(stru->at(0)); }}, + {"removed atom", [&] { stru->erase(0); }}, + {"symmetry operation", [&] { stru->addSymOp(R, R3::zerovector); }}, + {"symmetry precision", [&] { stru->setSymmetryPrecision(1.0); }}, + }; + for (const auto& [name, change] : changes) { + INFO("changed " << name); + change(); + const std::vector now = bonds(*stru); + CHECK(!same(last, now)); + checkSameBonds(now, bonds(*rebuilt(*stru))); + last = now; + } + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_crystal_refresh", + "[BondGeneratorBatching]") { + // A generator keeps its bonds when the structure expands symmetry again. + testAllKernels([] { + CrystalStructureAdapterPtr stru = loadTestCrystalStructure("CaTiO3.stru"); + BaseBondGeneratorPtr bnds = nearSite1(*stru); + const std::vector expected = generated(*bnds); + REQUIRE(!expected.empty()); + // An unchanged structure and its copy keep their symmetry positions when + // they create another generator. + const Atom* symatoms = stru->getEquivalentAtoms(0).data(); + BaseBondGeneratorPtr bnds1 = stru->createBondGenerator(); + CHECK(symatoms == stru->getEquivalentAtoms(0).data()); + auto copy = std::make_shared(*stru); + const Atom* copyatoms = copy->getEquivalentAtoms(0).data(); + BaseBondGeneratorPtr bnds2 = copy->createBondGenerator(); + CHECK(copyatoms == copy->getEquivalentAtoms(0).data()); + // changes of the structure do not reach the first generator + stru->updateSymmetryPositions(); + checkSameBonds(generated(*bnds), expected); + stru->at(2).uij_cartn(0, 0) += 0.01; + stru->at(3).xyz_cartn[1] += 0.2; + bnds1 = stru->createBondGenerator(); + checkSameBonds(generated(*bnds), expected); + stru->addSymOp(R3::identity(), R3::Vector(0.5, 0.5, 0.5)); + CHECK(stru->siteMultiplicity(0) > 0); + checkSameBonds(generated(*bnds), expected); + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_crystal_constructor", + "[BondGeneratorBatching]") { + // The public constructor finds the symmetry positions when needed. + testAllKernels([] { + CrystalStructureAdapterPtr stru = loadTestCrystalStructure("CaTiO3.stru"); + CrystalStructureBondGenerator direct(stru); + nearSite1(direct); + const std::vector expected = generated(*nearSite1(*stru)); + REQUIRE(!expected.empty()); + checkSameBonds(generated(direct), expected); + }); + // a structure without symmetry operations is rejected + StructureAdapterPtr periodic = loadTestPeriodicStructure("CaTiO3.stru"); + CHECK_THROWS_AS(CrystalStructureBondGenerator(periodic), + std::invalid_argument); + // The generators point into their own arrays and cannot be copied. + static_assert(!std::is_copy_constructible_v); + static_assert(!std::is_copy_assignable_v); + static_assert(!std::is_copy_constructible_v); +} + +TEST_CASE("TestBondGeneratorBatching::test_crystal_multiplicity", + "[BondGeneratorBatching]") { + // Site multiplicities follow atoms moved in place to and from a special + // position without any bond generator. + auto stru = std::make_shared(); + stru->setLatPar(4.0, 5.0, 6.0, 90.0, 90.0, 90.0); + stru->addSymOp(R3::identity(), R3::zerovector); + stru->addSymOp(-1.0 * R3::identity(), R3::zerovector); + Atom a; + a.atomtype = "C"; + stru->append(a); + stru->append(a); + stru->at(1).xyz_cartn = R3::Vector(1.0, 1.0, 1.0); + CHECK(1 == stru->siteMultiplicity(0)); + CHECK(2 == stru->siteMultiplicity(1)); + CHECK(3.0 == stru->totalOccupancy()); + stru->at(0).xyz_cartn = R3::Vector(0.4, 0.5, 0.6); + CHECK(2 == stru->siteMultiplicity(0)); + CHECK(2u == stru->getEquivalentAtoms(0).size()); + CHECK(4.0 == stru->totalOccupancy()); + // fractional (1/2, 1/2, 1/2) is a special position as well + stru->at(1).xyz_cartn = R3::Vector(2.0, 2.5, 3.0); + CHECK(1 == stru->siteMultiplicity(1)); + CHECK(3.0 == stru->totalOccupancy()); + BaseBondGeneratorPtr bnds = stru->createBondGenerator(); + for (int anchor : {0, 1}) { + bnds->selectAnchorSite(anchor); + CHECK(stru->siteMultiplicity(anchor) == bnds->multiplicity()); + } +} + +TEST_CASE("TestBondGeneratorBatching::test_empty_first_rewind", + "[BondGeneratorBatching]") { + // Empty selection at the first rewind. + testAllKernels([] { + StructureAdapterPtr p = loadTestPeriodicStructure("CaTiO3.stru"); + const PeriodicStructureAdapter& stru = + dynamic_cast(*p); + CrystalStructureAdapterPtr cstru = loadTestCrystalStructure("CaTiO3.stru"); + BaseBondGeneratorPtr bnds = stru.createBondGenerator(); + BaseBondGeneratorPtr cbnds = cstru->createBondGenerator(); + for (BaseBondGeneratorPtr b : {bnds, cbnds}) { + b->setRmax(300.0); + b->selectAnchorSite(1); + b->selectSiteRange(2, 2); + b->rewind(); + CHECK(b->finished()); + b->setRmax(6.0); + } + const Config cfg = {0.0, 6.0}; + SiteIndices sel = {0, 3, 7}; + bnds->selectSites(sel); + cbnds->selectSites(sel); + checkSameBonds(generated(*bnds), referencePeriodic(stru, 1, sel, cfg)); + checkSameBonds(generated(*cbnds), referenceCrystal(*cstru, 1, sel, cfg)); + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_range_bounds", + "[BondGeneratorBatching]") { + // Bonds exactly at rmin or rmax are included, one ulp beyond are not. + StructureAdapterPtr p = loadTestPeriodicStructure("ZnS_wurtzite.stru"); + CrystalStructureAdapterPtr c = loadTestCrystalStructure("ZnS_wurtzite.stru"); + testAllKernels([&] { + for (StructureAdapterPtr stru : {p, StructureAdapterPtr(c)}) { + BaseBondGeneratorPtr bnds = stru->createBondGenerator(); + bnds->setRmax(5.0); + bnds->selectAnchorSite(0); + bnds->selectSiteRange(0, stru->countSites()); + const std::vector all = generated(*bnds); + REQUIRE(all.size() > 10); + for (size_t ib : {size_t(0), all.size() / 2, all.size() - 1}) { + const double d = all[ib].distance; + INFO("distance " << d); + auto count = [&](double rmin, double rmax) { + bnds->setRmin(rmin); + bnds->setRmax(rmax); + int rv = 0; + for (const Bond& b : generated(*bnds)) { + rv += (b.site1 == all[ib].site1 && b.r01[0] == all[ib].r01[0] && + b.r01[1] == all[ib].r01[1] && b.r01[2] == all[ib].r01[2]); + } + return rv; + }; + const double below = std::nextafter(d, 0.0); + const double above = std::nextafter(d, 2 * d); + CHECK(1 == count(0.0, d)); + CHECK(0 == count(0.0, below)); + CHECK(1 == count(d, 5.0)); + CHECK(0 == count(above, 5.0)); + CHECK(1 == count(d, d)); + CHECK(1 == count(below, above)); + } + } + }); +} + +TEST_CASE("TestBondGeneratorBatching::test_large_sphere", + "[BondGeneratorBatching]") { + // Spheres with many points and few candidates per point. + StructureAdapterPtr p = loadTestPeriodicStructure("Ni.stru"); + const PeriodicStructureAdapter& stru = + dynamic_cast(*p); + CrystalStructureAdapterPtr cstru = loadTestCrystalStructure("Ni.stru"); + const Config cfg = {100.0, 110.0}; + const Lattice& L = stru.getLattice(); + const double buffzone = L.ucMaxDiagonalLength(); + PointsInSphere sph(cfg.rmin - buffzone, cfg.rmax + buffzone, L); + int npts = 0; + for (sph.rewind(); !sph.finished(); sph.next()) ++npts; + REQUIRE(npts > (1 << 16)); + SiteIndices sel = {0, 1, 2, 3}; + testAllKernels([&] { + BaseBondGeneratorPtr bnds = stru.createBondGenerator(); + BaseBondGeneratorPtr cbnds = cstru->createBondGenerator(); + for (BaseBondGeneratorPtr b : {bnds, cbnds}) { + b->setRmin(cfg.rmin); + b->setRmax(cfg.rmax); + b->selectSites(sel); + } + for (int anchor : {0, 3}) { + INFO("anchor " << anchor); + bnds->selectAnchorSite(anchor); + cbnds->selectAnchorSite(anchor); + const std::vector expected = + referencePeriodic(stru, anchor, sel, cfg); + REQUIRE(!expected.empty()); + checkSameBonds(generated(*bnds), expected); + checkSameBonds(generated(*bnds), expected); + checkSameBonds(generated(*cbnds), + referenceCrystal(*cstru, anchor, sel, cfg)); + } + }); +} + +namespace { + +/// Set isotropic displacement parameters of an atom. +void setUiso(Atom& a, double uiso) { + a.anisotropy = false; + a.uij_cartn = R3::zeromatrix(); + for (int k = 0; k < 3; ++k) a.uij_cartn(k, k) = uiso; +} + +/// Add two atoms with Uiso = 0.01, which are 1 apart along the x-axis. +template +void appendDimer(Tstru& stru) { + Atom a; + setUiso(a, 0.01); + stru.append(a); + a.xyz_cartn[0] = 1.0; + stru.append(a); +} + +/// msd() at the only bond of the dimer after a rewind. +double dimerMSD(BaseBondGenerator& bnds) { + bnds.rewind(); + REQUIRE(!bnds.finished()); + REQUIRE(1 == bnds.site1()); + return bnds.msd(); +} + +} // namespace + +TEST_CASE("TestBondGeneratorBatching::test_msd", "[BondGeneratorBatching]") { + // msd() of a bond along (1, 2, 2) for isotropic, mixed and anisotropic + // sites in every kind of bond generator. The anisotropic tensor gives + // s'Us / s's = 0.18 / 9 = 0.02 along the bond. + testAllKernels([] { + R3::Matrix U = R3::zeromatrix(); + U(0, 0) = 0.01; + U(1, 1) = 0.03; + U(2, 2) = 0.01; + U(0, 1) = U(1, 0) = 0.0025; + auto atomic = std::make_shared(); + auto periodic = std::make_shared(); + periodic->setLatPar(10.0, 10.0, 10.0, 90.0, 90.0, 90.0); + auto crystal = std::make_shared(); + crystal->setLatPar(10.0, 10.0, 10.0, 90.0, 90.0, 90.0); + crystal->addSymOp(R3::identity(), R3::zerovector); + const std::pair> + adapters[] = { + {"atomic", atomic}, {"periodic", periodic}, {"crystal", crystal}}; + const double eps = 1e-15; + for (const auto& [name, stru] : adapters) { + INFO("adapter " << name); + appendDimer(*stru); + stru->at(1).xyz_cartn = R3::Vector(1.0, 2.0, 2.0); + // a new generator, the crystal one keeps earlier symmetry positions + auto msd = [&stru] { + BaseBondGeneratorPtr bnds = stru->createBondGenerator(); + bnds->setRmax(3.5); + bnds->selectAnchorSite(0); + return dimerMSD(*bnds); + }; + CHECK(std::fabs(msd() - 0.02) < eps); + stru->at(1).anisotropy = true; + stru->at(1).uij_cartn = U; + CHECK(std::fabs(msd() - 0.03) < eps); + stru->at(0).anisotropy = true; + stru->at(0).uij_cartn = U; + CHECK(std::fabs(msd() - 0.04) < eps); + } + }); +} + +// End of file diff --git a/src/tests/TestPQEvaluator.cpp b/src/tests/TestPQEvaluator.cpp index c2884f1..08cb8c5 100644 --- a/src/tests/TestPQEvaluator.cpp +++ b/src/tests/TestPQEvaluator.cpp @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -45,6 +46,48 @@ class BadPairCounter : public PairCounter { virtual void restorePartialValue() {} }; +// calculator that fails after a given number of pair contributions + +class FailingPDFCalculator : public PDFCalculator { + public: + /// Pair contributions left before the failure, negative for no failure. + int countdown = -1; + /// Pair contributions added so far. + int calls = 0; + + protected: + virtual void addPairContribution(const BaseBondGenerator& bnds, + int summationscale) { + if (countdown == 0) throw runtime_error("failed pair contribution"); + if (countdown > 0) --countdown; + ++calls; + this->PDFCalculator::addPairContribution(bnds, summationscale); + } +}; + +// Crystal with an inversion center at the origin. The first atom sits +// on it at x = 0 and has two symmetry positions elsewhere. + +CrystalStructureAdapterPtr inversionCrystal(double x) { + CrystalStructureAdapterPtr rv(new CrystalStructureAdapter); + rv->setLatPar(4.0, 4.5, 5.0, 90.0, 90.0, 90.0); + rv->addSymOp(R3::identity(), R3::zerovector); + rv->addSymOp(-1.0 * R3::identity(), R3::zerovector); + const double positions[][3] = {{x, 0.0, 0.0}, + {0.9, 0.6, 1.3}, + {1.7, 1.2, 0.5}, + {0.4, 1.9, 2.1}, + {1.3, 0.3, 1.9}}; + Atom a; + a.atomtype = "Ni"; + a.uij_cartn = 0.005 * R3::identity(); + for (const auto& xyz : positions) { + a.xyz_cartn = R3::Vector(xyz[0], xyz[1], xyz[2]); + rv->append(a); + } + return rv; +} + ////////////////////////////////////////////////////////////////////////////// // TestPQEvaluator ////////////////////////////////////////////////////////////////////////////// @@ -272,6 +315,84 @@ TEST_CASE_METHOD(TestPQEvaluatorFixture, CHECK(PQEvaluatorType::CHECK == badcounter.getEvaluatorTypeUsed()); } +TEST_CASE_METHOD(TestPQEvaluatorFixture, "TestPQEvaluator::test_failed_update", + "[PQEvaluator]") { + // A value left by a failed update must not pass for a valid one. + mpdfcb.eval(mstru10); + const QuantityType v10 = mpdfcb.value(); + mpdfcb.eval(mstru10d1); + const QuantityType v10d1 = mpdfcb.value(); + // count the pair contributions in a fast update + FailingPDFCalculator counter; + counter.eval(mstru10); + counter.calls = 0; + counter.eval(mstru10d1); + REQUIRE(OPTIMIZED == counter.getEvaluatorTypeUsed()); + const int nfast = counter.calls; + REQUIRE(nfast > 4); + FailingPDFCalculator pdfc; + SECTION("full calculation") { + // all atoms are different, which needs a full calculation + AtomicStructureAdapterPtr shifted = + std::make_shared(*mstru10); + for (Atom& a : *shifted) a.xyz_cartn[1] += 0.5; + pdfc.eval(mstru10); + pdfc.countdown = 3; + CHECK_THROWS_AS(pdfc.eval(shifted), runtime_error); + pdfc.countdown = -1; + pdfc.eval(mstru10); + CHECK(BASIC == pdfc.getEvaluatorTypeUsed()); + CHECK(allclose(v10, pdfc.value())); + } + SECTION("fast update") { + pdfc.eval(mstru10); + pdfc.countdown = nfast / 2; + CHECK_THROWS_AS(pdfc.eval(mstru10d1), runtime_error); + pdfc.countdown = -1; + pdfc.eval(mstru10); + CHECK(BASIC == pdfc.getEvaluatorTypeUsed()); + CHECK(allclose(v10, pdfc.value())); + } + SECTION("check of a fast update") { + // fail in the full calculation that follows the fast update + pdfc.setEvaluatorType(PQEvaluatorType::CHECK); + pdfc.eval(mstru10); + pdfc.countdown = nfast + 2; + CHECK_THROWS_AS(pdfc.eval(mstru10d1), runtime_error); + pdfc.countdown = -1; + pdfc.eval(mstru10d1); + CHECK(allclose(v10d1, pdfc.value())); + } +} + +TEST_CASE_METHOD(TestPQEvaluatorFixture, + "TestPQEvaluator::test_structure_changed_in_place", + "[PQEvaluator]") { + // The first evaluation after an atom moves off a special position must + // include its new symmetry positions. + mpdfcb.setRmax(7.0); + mpdfcb.eval(inversionCrystal(0.3)); + const QuantityType expected = mpdfcb.getPDF(); + for (PQEvaluatorType tp : {BASIC, OPTIMIZED, PQEvaluatorType::CHECK}) { + INFO("evaluator " << tp); + CrystalStructureAdapterPtr stru = inversionCrystal(0.0); + PDFCalculator pdfc; + pdfc.setEvaluatorType(tp); + pdfc.setRmax(7.0); + pdfc.eval(stru); + REQUIRE(1 == stru->siteMultiplicity(0)); + REQUIRE(!allclose(expected, pdfc.getPDF())); + stru->at(0).xyz_cartn[0] = 0.3; + for (int i = 0; i < 2; ++i) { + pdfc.eval(stru); + // the first evaluation is a fast update unless the type is BASIC + if (i == 0) CHECK(tp == pdfc.getEvaluatorTypeUsed()); + CHECK(2 == stru->siteMultiplicity(0)); + CHECK(allclose(expected, pdfc.getPDF())); + } + } +} + } // namespace srreal } // namespace diffpy From e227189888d1d814c98e1cf60c4cb26c76122bbd Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Tue, 10 Nov 2026 21:49:09 -0800 Subject: [PATCH 4/6] build: optimize the shared library --- CMakeLists.txt | 24 +++++++++++++++++++++++- docs/source/development.rst | 3 ++- news/faster-pdf.rst | 4 ++++ 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0564a61..a78169d 100644 --- a/CMakeLists.txt +++ b/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,7 +342,6 @@ target_link_libraries(diffpy if(LIBDIFFPY_HAS_OBJCRYST) target_link_libraries(diffpy PUBLIC ObjCryst::ObjCryst) endif() - include(CalculationKernels) # Installation ---------------------------------------------------------------- diff --git a/docs/source/development.rst b/docs/source/development.rst index 41db078..465fe30 100644 --- a/docs/source/development.rst +++ b/docs/source/development.rst @@ -45,7 +45,8 @@ 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 -G ""``. 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 ------------------- diff --git a/news/faster-pdf.rst b/news/faster-pdf.rst index 51964e7..f2b1ec7 100644 --- a/news/faster-pdf.rst +++ b/news/faster-pdf.rst @@ -28,6 +28,10 @@ * 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. +* 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:** From 07a2a0403467b250e456b023fa5e8f6fa90795d6 Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Wed, 11 Nov 2026 12:12:23 -0800 Subject: [PATCH 5/6] perf: faster fft with cache --- news/faster-pdf.rst | 1 + src/diffpy/srreal/PDFUtils.cpp | 138 +++++++++++++++++-------- src/tests/CMakeLists.txt | 4 +- src/tests/TestPDFUtils.cpp | 184 +++++++++++++++++++++++++++++++++ 4 files changed, 282 insertions(+), 45 deletions(-) create mode 100644 src/tests/TestPDFUtils.cpp diff --git a/news/faster-pdf.rst b/news/faster-pdf.rst index f2b1ec7..918b974 100644 --- a/news/faster-pdf.rst +++ b/news/faster-pdf.rst @@ -28,6 +28,7 @@ * 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 diff --git a/src/diffpy/srreal/PDFUtils.cpp b/src/diffpy/srreal/PDFUtils.cpp index 70dbdf7..038de3d 100644 --- a/src/diffpy/srreal/PDFUtils.cpp +++ b/src/diffpy/srreal/PDFUtils.cpp @@ -16,74 +16,124 @@ * *****************************************************************************/ +#include +#include +#include +#include +#include +#include +#include +#include #include +#include #include #include #include -#include -#include using namespace std; -using namespace diffpy::validators; -using namespace diffpy::mathutils; namespace diffpy { namespace srreal { -// Local Constants ----------------------------------------------------------- - -namespace { - -const char* EMSGFFT = "Fourier Transformation failed."; - -} // namespace - // PDFUtils functions -------------------------------------------------------- QuantityType fftgtof(const QuantityType& g, double rstep, double rmin) { if (g.empty()) return g; - int padrmin = int(round(rmin / rstep)); - int Npad1 = padrmin + g.size(); - // pad to the next power of 2 for fast Fourier transformation - int Npad2 = (1 << int(ceil(log2(Npad1)))); - // sine transformations needs an odd extension - // gpadc array has to be doubled for complex coefficients - int Npad4 = 4 * Npad2; - valarray gpadc(0.0, Npad4); - QuantityType::const_iterator gi; - // copy the original g signal - int ilo = padrmin; - for (gi = g.begin(); gi != g.end(); ++gi, ++ilo) { - gpadc[2 * ilo] = *gi; + // std::isfinite is always true with -ffast-math. + const auto isfinitebits = [](double v) { + const uint64_t expmask = 0x7ff0000000000000; + return (bit_cast(v) & expmask) != expmask; + }; + if (!(isfinitebits(rstep) && rstep > 0.0)) { + throw invalid_argument("rstep must be finite and positive."); + } + if (!isfinitebits(rmin)) { + throw invalid_argument("rmin must be finite."); + } + // rmin within half a step below zero rounds to no padding. + const double pad = round(rmin / rstep); + if (pad < 0.0) { + throw invalid_argument("rmin rounds to a negative offset."); } - // copy the odd part of g skipping the first point, - // because it is periodic image of gpadc[0] - int ihi = 2 * Npad2 - 1; - for (ilo = 1; ilo < Npad2; ++ilo, --ihi) { - gpadc[2 * ihi] = -1 * gpadc[2 * ilo]; + // Keep size_t(pad), bit_ceil and the FFT buffer of 2*N doubles in range. + const size_t maxpoints = bit_floor(vector().max_size() / 2); + if (!isfinitebits(pad) || pad > maxpoints || + g.size() > maxpoints - size_t(pad)) { + throw length_error("Too many points for the Fourier transformation."); } - int status; - status = gsl_fft_complex_radix2_inverse(&(gpadc[0]), 1, 2 * Npad2); - if (status != GSL_SUCCESS) throw invalid_argument(EMSGFFT); - QuantityType f(Npad2); - for (int i = 0; i < Npad2; ++i) { - f[i] = gpadc[2 * i + 1] * Npad2 * rstep; + const size_t padding = size_t(pad); + const size_t N = bit_ceil(padding + g.size()); + QuantityType f(N, 0.0); + copy(g.begin(), g.end(), f.begin() + padding); + // sin(0) removes f[0] from the sum; zero it to keep it out of round-off. + f[0] = 0.0; + if (N < 4) { + if (N == 2) f[1] *= rstep; + return f; + } + + // Per-thread FFT buffers and tables for each power-of-two length. + struct FFTCache { + vector data, wr, wi; + unique_ptr + table; + unique_ptr + work; + + explicit FFTCache(size_t npts) + : data(2 * npts), + wr(npts / 2 + 1), + wi(npts / 2 + 1), + table(gsl_fft_complex_wavetable_alloc(npts), + gsl_fft_complex_wavetable_free), + work(gsl_fft_complex_workspace_alloc(npts), + gsl_fft_complex_workspace_free) { + if (!table || !work) throw bad_alloc(); + for (size_t j = 0; j <= npts / 2; ++j) { + wr[j] = cos(M_PI * j / npts); + wi[j] = sin(M_PI * j / npts); + } + } + }; + thread_local array, numeric_limits::digits> + transforms; + auto& plan = transforms[countr_zero(N)]; + if (!plan) plan.emplace(N); + FFTCache& cache = *plan; + + // Need Im X[m] for X[m] = sum_k f[k] exp(i*pi*m*k/N). The N-point FFT Z + // of z[j] = f[2*j] + i*f[2*j+1], padded with zeros, is enough. + copy(f.begin(), f.end(), cache.data.begin()); + fill(cache.data.begin() + N, cache.data.end(), 0.0); + const int status = gsl_fft_complex_backward( + cache.data.data(), 1, N, cache.table.get(), cache.work.get()); + if (status != GSL_SUCCESS) + throw invalid_argument("Fourier Transformation failed."); + + // Z[m] = E[m] + i*O[m] and conj(Z[N-m]) = E[m] - i*O[m] yield the even + // and odd sample transforms E, O; X[m] = E[m] + exp(i*pi*m/N) * O[m]. + const double c = 0.5 * rstep; + const auto& z = cache.data; + for (size_t m = 1; m <= N / 2; ++m) { + const size_t mc = N - m; + const double ei = c * (z[2 * m + 1] - z[2 * mc + 1]); + const double orr = c * (z[2 * m + 1] + z[2 * mc + 1]); + const double oi = c * (z[2 * mc] - z[2 * m]); + const double t = cache.wr[m] * oi + cache.wi[m] * orr; + f[m] = ei + t; + f[mc] = t - ei; } - // real components should be all close to zero - assert(fabs(valarray(gpadc[slice(0, 2 * Npad2, 2)]).min()) <= - SQRT_DOUBLE_EPS * gpadc.max()); - assert(fabs(valarray(gpadc[slice(0, 2 * Npad2, 2)]).max()) <= - SQRT_DOUBLE_EPS * gpadc.max()); return f; } QuantityType fftftog(const QuantityType& f, double qstep, double qmin) { - // ftog is the same as gtof, just normalized by 2 / PI + // ftog is the same as gtof, just normalized by 2 / PI. QuantityType g = fftgtof(f, qstep, qmin); - QuantityType::iterator gi; - for (gi = g.begin(); gi != g.end(); ++gi) *gi *= 2.0 / M_PI; + for (double& gi : g) gi *= 2.0 / M_PI; return g; } diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 5f879eb..08632a3 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -4,6 +4,7 @@ # suite separately, preserving filtering without executing test discovery. find_package(Catch2 3.5 CONFIG REQUIRED) +find_package(Threads REQUIRED) message(STATUS "Using Catch2 ${Catch2_VERSION} from ${Catch2_DIR}") # Select the test suites ----------------------------------------------------- @@ -58,7 +59,8 @@ endif() set_property(SOURCE TestCalculationKernel.cpp test_helpers.cpp APPEND PROPERTY COMPILE_OPTIONS "${_test_math_options}") target_include_directories(alltests PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}") -target_link_libraries(alltests PRIVATE libdiffpy::diffpy Catch2::Catch2WithMain) +target_link_libraries(alltests PRIVATE + libdiffpy::diffpy Catch2::Catch2WithMain Threads::Threads) # test_helpers.cpp finds the testdata directory from this macro. set_property(SOURCE test_helpers.cpp APPEND PROPERTY COMPILE_DEFINITIONS "DIFFPYTESTSDIRPATH=${CMAKE_CURRENT_SOURCE_DIR}") diff --git a/src/tests/TestPDFUtils.cpp b/src/tests/TestPDFUtils.cpp new file mode 100644 index 0000000..f4fbd87 --- /dev/null +++ b/src/tests/TestPDFUtils.cpp @@ -0,0 +1,184 @@ +/***************************************************************************** + * + * libdiffpy Complex Modeling Initiative + * (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.rst for license information. + * + ****************************************************************************** + * + * TestPDFUtils -- unit test suite for the Fourier transformations + * fftgtof and fftftog + * + *****************************************************************************/ + +#include +#include +#include +#include + +#include + +#include + +#include "test_helpers.hpp" + +using namespace std; +using namespace diffpy::srreal; + +namespace { + +/// Tolerance relative to the largest value, round-off is below 1e-13 here. +const double RELTOL = 1e-12; + +/// oscillating signal without any symmetry +QuantityType makeSignal(int npts) { + QuantityType g(npts); + for (int i = 0; i < npts; ++i) { + g[i] = sin(0.37 * i + 0.2) * exp(-0.002 * i) + 0.3 * cos(1.9 * i); + } + return g; +} + +/// direct summation of the sine transformation evaluated by fftgtof +QuantityType sumgtof(const QuantityType& g, double rstep, int padrmin = 0) { + int Npad2 = 1; + while (Npad2 < padrmin + int(g.size())) Npad2 *= 2; + QuantityType f(Npad2); + for (int j = 0; j < Npad2; ++j) { + double fj = 0.0; + for (int i = 0; i < int(g.size()); ++i) { + // reduce the sine argument for an accurate result at large j * k + const int jk = j * (padrmin + i) % (2 * Npad2); + fj += g[i] * sin(M_PI * jk / Npad2); + } + f[j] = fj * rstep; + } + return f; +} + +} // namespace + +TEST_CASE("TestPDFUtils::test_fftgtof_empty", "[PDFUtils]") { + QuantityType g; + CHECK(fftgtof(g, 0.01).empty()); + CHECK(fftftog(g, 0.01).empty()); +} + +TEST_CASE("TestPDFUtils::test_fftgtof_small", "[PDFUtils]") { + const double rstep = 0.25; + QuantityType g(1, 3.0); + QuantityType f = fftgtof(g, rstep); + REQUIRE(1u == f.size()); + CHECK(0.0 == f[0]); + g.push_back(5.0); + f = fftgtof(g, rstep); + REQUIRE(2u == f.size()); + CHECK(0.0 == f[0]); + CHECK(5.0 * rstep == f[1]); + g.push_back(-2.0); + f = fftgtof(g, rstep); + REQUIRE(4u == f.size()); + CHECK(maxRelativeError(f, sumgtof(g, rstep)) < RELTOL); +} + +TEST_CASE("TestPDFUtils::test_fftgtof_sizes", "[PDFUtils]") { + const double rstep = 0.01; + // one length for each padded size from 4 to 1024 + const int sizes[] = {4, 7, 9, 32, 33, 100, 256, 500, 999}; + for (int npts : sizes) { + INFO("npts = " << npts); + QuantityType g = makeSignal(npts); + QuantityType f = fftgtof(g, rstep); + QuantityType fsum = sumgtof(g, rstep); + REQUIRE(fsum.size() == f.size()); + CHECK(0.0 == f[0]); + CHECK(maxRelativeError(f, fsum) < RELTOL); + } +} + +TEST_CASE("TestPDFUtils::test_fftgtof_rmin", "[PDFUtils]") { + const double rstep = 0.02; + QuantityType g = makeSignal(200); + // rmin is rounded to the nearest multiple of rstep + const double rmins[] = {0.0, 0.02, 0.029, 0.031, 1.0, 1.12, 5.6}; + const int padrmins[] = {0, 1, 1, 2, 50, 56, 280}; + for (int i = 0; i < 7; ++i) { + INFO("rmin = " << rmins[i]); + QuantityType f = fftgtof(g, rstep, rmins[i]); + QuantityType fsum = sumgtof(g, rstep, padrmins[i]); + REQUIRE(fsum.size() == f.size()); + CHECK(maxRelativeError(f, fsum) < RELTOL); + } + // shifted signal is equivalent to explicit zero padding + QuantityType gpad(56, 0.0); + gpad.insert(gpad.end(), g.begin(), g.end()); + CHECK(fftgtof(gpad, rstep) == fftgtof(g, rstep, 1.12)); +} + +TEST_CASE("TestPDFUtils::test_fftgtof_edge_points", "[PDFUtils]") { + // Edge-point signals have the largest round-off; it must not grow with N. + const double rstep = 0.01; + for (int npts : {2048, 1 << 16}) { + for (int k : {1, 2, npts - 2, npts - 1}) { + INFO("npts = " << npts << ", signal at point " << k); + QuantityType g(npts, 0.0); + g[k] = 1.0; + // the exact transformation of a single point + QuantityType fexact(npts); + for (int j = 0; j < npts; ++j) { + const long long jk = (static_cast(j) * k) % (2 * npts); + fexact[j] = rstep * sin(M_PI * double(jk) / npts); + } + CHECK(maxRelativeError(fftgtof(g, rstep), fexact) < 1e-13); + // The zero-coordinate sample must not contaminate the sine sum. + g[0] = 1e100; + CHECK(maxRelativeError(fftgtof(g, rstep), fexact) < 1e-13); + } + } +} + +TEST_CASE("TestPDFUtils::test_fftftog", "[PDFUtils]") { + const double qstep = 0.05; + QuantityType f = makeSignal(300); + QuantityType g = fftftog(f, qstep, 0.5); + QuantityType gsum = sumgtof(f, qstep, 10); + for (double& gi : gsum) gi *= 2.0 / M_PI; + REQUIRE(gsum.size() == g.size()); + CHECK(maxRelativeError(g, gsum) < RELTOL); +} + +TEST_CASE("TestPDFUtils::test_fft_threads", "[PDFUtils]") { + // Same length, different data and steps: threads must not share FFT state. + constexpr int nthreads = 4; + constexpr int npts = 1024; + array signals, expected; + array steps; + for (int i = 0; i < nthreads; ++i) { + signals[i] = makeSignal(npts); + for (double& v : signals[i]) v = (i + 1) * v + 0.1 * i; + signals[i][0] = 0.0; + steps[i] = 0.01 * (i + 1); + expected[i] = sumgtof(signals[i], steps[i]); + } + barrier ready(nthreads); + array, nthreads> workers; + for (int i = 0; i < nthreads; ++i) { + workers[i] = async(launch::async, [&, i] { + ready.arrive_and_wait(); + double error = 0.0; + for (int repeat = 0; repeat < 16; ++repeat) { + const auto f = fftgtof(signals[i], steps[i]); + error = max(error, maxRelativeError(f, expected[i])); + const auto g = fftftog(f, M_PI / (npts * steps[i])); + error = max(error, maxRelativeError(g, signals[i])); + } + return error; + }); + } + for (auto& worker : workers) CHECK(worker.get() < RELTOL); +} + +// End of file From 980c6627b326d36dda2fe83757618638e1adeb3b Mon Sep 17 00:00:00 2001 From: SYMSCAE Date: Thu, 12 Nov 2026 23:04:34 -0800 Subject: [PATCH 6/6] fix: improve gaussian peak numerical stability --- .../srreal/kernels/CalculationKernelImpl.hpp | 3 +- .../srreal/kernels/CalculationKernelSIMD.ipp | 45 ++++++------------- .../kernels/CalculationKernelScalar.cpp | 3 +- src/tests/TestCalculationKernel.cpp | 30 ------------- 4 files changed, 18 insertions(+), 63 deletions(-) diff --git a/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp b/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp index 5a32f94..8e1a3d5 100644 --- a/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp +++ b/src/diffpy/srreal/kernels/CalculationKernelImpl.hpp @@ -26,7 +26,8 @@ namespace kernels { /// Add one Gaussian peak to the PDFCalculator RDF array. /// -/// For each i in [ilo, ihi) with x = (rlosteps + i) * rstep - dist: +/// For each i in [ilo, ihi) with x = xlo + (i - ilo) * rstep, where +/// xlo = (rlosteps + ilo) * rstep - dist: /// out[i] += scale * G(x, fwhm) * (x / dist + 1) /// where G(x, fwhm) = peakheight * exp(-4 * ln(2) * (x / fwhm)^2), or zero /// for |x / fwhm| >= halfboundrel when cropped. peakheight is the profile diff --git a/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp b/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp index bf2e5c4..9548423 100644 --- a/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp +++ b/src/diffpy/srreal/kernels/CalculationKernelSIMD.ipp @@ -40,7 +40,6 @@ #error "Define LIBDIFFPY_XSIMD_NAMESPACE, see the comment above." #endif -#include #include #include @@ -66,13 +65,6 @@ constexpr int RESEED_SLACK = 8; /// normal and rho finite. exp(-SAFE_EXPONENT) is about 1e-261. constexpr double SAFE_EXPONENT = 600.0; -/// True when xsimd::fma of Arch is a fused multiply-add. SSE2 has none and -/// xsimd evaluates u * v + w there. -template -constexpr bool FUSED_FMA = true; -template <> -constexpr bool FUSED_FMA = false; - template struct SIMDKernel { using B = xsimd::batch; @@ -111,11 +103,6 @@ struct SIMDKernel { for (int j = 0; j < n; ++j) out[j] += uv[j]; } - /// u * v + w for one value, rounded as in the lanes of xsimd::fma. - static double mulAdd(double u, double v, double w) { - return FUSED_FMA ? fma(u, v, w) : u * v + w; - } - static void gaussianWindow(const GaussianWindowArgs& args) { // Vector stores may alias anything, work on a copy of the arguments. const GaussianWindowArgs a = args; @@ -123,16 +110,13 @@ struct SIMDKernel { // exp(-b x^2) and its ratio over one vector step. const double b = 4 * M_LN2 / (a.fwhm * a.fwhm); const double wh = W * a.rstep; - // Use the scalar kernel for NaN, for windows too far from the peak, where - // b * x^2 or 2 * b * wh * |x| reach SAFE_EXPONENT, and for narrow peaks - // far out, where the grid rounding (~eps * dist) is too coarse. + // Use the scalar kernel for NaN and for windows too far from the peak, + // where b * x^2 or 2 * b * wh * |x| reach SAFE_EXPONENT. const double xlo = (a.rlosteps + a.ilo) * a.rstep - a.dist; const double xhi = (a.rlosteps + a.ihi) * a.rstep - a.dist; const double xfar = (-xlo > xhi) ? -xlo : xhi; const double reach = (xfar > 2 * wh) ? xfar : 2 * wh; - const bool recurrence = (b * xfar * reach < SAFE_EXPONENT) && - (DBL_EPSILON * a.dist <= 1e-13 * a.fwhm); - if (!recurrence) { + if (!(b * xfar * reach < SAFE_EXPONENT)) { scalarCalculationKernel().gaussianWindow(args); return; } @@ -142,27 +126,26 @@ struct SIMDKernel { const B ratiostep(exp(-2 * b * wh * wh)); const B vw(1.0 * W); const B vrstep(a.rstep); - const B vdist(a.dist); + const B vxlo(xlo); const B vinvdist(1.0 / a.dist); const B vnorm(a.scale * a.peakheight); const B vfwhm(a.fwhm); const B vhalfbound(a.halfboundrel); const B one(1.0); const B zero(0.0); - // Grid steps of the lanes, whole numbers are exact in double. - B steps = B(double(a.rlosteps + a.ilo)) + laneIndex(); - // x / dist + 1 is r / dist, which grows by a constant per vector step. - const B vrscale(a.rstep / a.dist); + // Grid steps of the lanes from ilo, whole numbers are exact in double. + B steps = laneIndex(); + // x / dist + 1 grows by a constant per vector step. const B vrscalestep(wh / a.dist); for (int i = a.ilo; i < a.ihi;) { const int iend = blockEnd(i, a.ihi); - const B xseed = steps * vrstep - vdist; + // Same x as the scalar kernel, so cropping agrees. + const B xseed = steps * vrstep + vxlo; B g = vnorm * xsimd::exp(negb * xseed * xseed); B ratio = xsimd::exp(negb * (twowh * xseed + whwh)); if (a.cropped) { for (; i < iend; i += W) { - // Same x as the scalar kernel, so cropping agrees. - const B x = steps * vrstep - vdist; + const B x = steps * vrstep + vxlo; const B y = xsimd::select(xsimd::abs(x / vfwhm) >= vhalfbound, zero, g); addProduct(a.out + i, y, xsimd::fma(x, vinvdist, one), iend - i); @@ -171,7 +154,7 @@ struct SIMDKernel { steps += vw; } } else { - B rscaled = steps * vrscale; + B rscaled = xsimd::fma(xseed, vinvdist, one); const int istart = i; for (; i < iend; i += W) { addProduct(a.out + i, g, rscaled, iend - i); @@ -255,13 +238,13 @@ struct SIMDKernel { m += int((bits >> j) & 1); } } - // The last points go one by one with the rounding of the vector lanes, - // a partial vector load would be much slower. + // The last points go one by one, a partial vector load would be much + // slower. for (; i < a.n; ++i) { const double dx = a.x[i] + a.ox; const double dy = a.y[i] + a.oy; const double dz = a.z[i] + a.oz; - const double d2 = mulAdd(dz, dz, mulAdd(dy, dy, dx * dx)); + const double d2 = dx * dx + dy * dy + dz * dz; a.idx[m] = i; m += int(!(d2 < a.d2lo || d2 > a.d2hi)); } diff --git a/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp b/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp index 5849f12..b0a1f80 100644 --- a/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp +++ b/src/diffpy/srreal/kernels/CalculationKernelScalar.cpp @@ -26,8 +26,9 @@ namespace { void gaussianWindowScalar(const GaussianWindowArgs& a) { if (a.fwhm <= 0.0) return; + const double xlo = (a.rlosteps + a.ilo) * a.rstep - a.dist; for (int i = a.ilo; i < a.ihi; ++i) { - const double x = (a.rlosteps + i) * a.rstep - a.dist; + const double x = xlo + (i - a.ilo) * a.rstep; const double xrel = x / a.fwhm; if (a.cropped && fabs(xrel) >= a.halfboundrel) continue; const double y = a.peakheight * exp(-4 * std::numbers::ln2 * xrel * xrel); diff --git a/src/tests/TestCalculationKernel.cpp b/src/tests/TestCalculationKernel.cpp index faa7fc8..7f3c7a0 100644 --- a/src/tests/TestCalculationKernel.cpp +++ b/src/tests/TestCalculationKernel.cpp @@ -444,36 +444,6 @@ TEST_CASE("TestCalculationKernel::test_bond_candidates", } } } - // All points are the same, so all or none are candidates. With the - // bound between the fused and the separately rounded x^2 + z^2, the - // count is partial when the last points round unlike the vector lanes. - int nrounded = 0; - for (int k = 1; k < 200; ++k) { - const double vx = 1.0 + k * 0x1p-30; - const double vz = 3.0 - k * 0x1p-31; - const volatile double xx = vx * vx; - const volatile double zz = vz * vz; - const double fused = std::fma(vz, vz, xx); - const double separate = zz + xx; - if (fused == separate) continue; - ++nrounded; - for (int n : {9, 15, 17}) { - INFO("k = " << k << ", n = " << n); - std::vector x(n, vx), y(n, 0.0), z(n, vz); - std::vector idx(n); - const diffpy::srreal::kernels::BondCandidatesArgs args = { - .x = x.data(), - .y = y.data(), - .z = z.data(), - .n = n, - .d2hi = std::min(fused, separate), - .idx = idx.data(), - }; - const int m = kernel.bondCandidates(args); - CHECK((m == 0 || m == n)); - } - } - CHECK(nrounded > 0); std::vector x(9, 1.0), y(9, -2.0), z(9, 2.0); std::vector idx(9); diffpy::srreal::kernels::BondCandidatesArgs args = {