perf: Parallelize siteRDF with ProductIterator

benchmark/CMakeLists.txt

@@ -10,6 +10,7 @@ dissolve_benchmark(energy)
 dissolve_benchmark(forces)
 dissolve_benchmark(region)
 dissolve_benchmark(math)
+dissolve_benchmark(siterdf)
 
 # create a single benchmark executable
 add_executable(benchmarks benchmarks.cpp ${benchmark_files})

benchmark/benchmark_files/argon10k.txt

@@ -5,6 +5,10 @@ Species 'Ball'
   # Atoms
   Atom    1    Ar     0.0  0.0  0.0  'Ar'
 
+  Site 'Ar'
+    Origin  1
+  EndSite
+
   Isotopologue  'Natural'
 EndSpecies
 

benchmark/benchmark_files/argon1k.txt

@@ -5,6 +5,10 @@ Species 'Ball'
   # Atoms
   Atom    1    Ar     0.0  0.0  0.0  'Ar'
 
+  Site 'Ar'
+    Origin  1
+  EndSite
+
   Isotopologue  'Natural'
 EndSpecies
 

benchmark/benchmark_files/argon5k.txt

@@ -5,6 +5,10 @@ Species 'Ball'
   # Atoms
   Atom    1    Ar     0.0  0.0  0.0  'Ar'
 
+  Site 'Ar'
+    Origin  1
+  EndSite
+
   Isotopologue  'Natural'
 EndSpecies
 

benchmark/benchmark_files/water1k.txt

@@ -25,6 +25,10 @@ Species 'Water'
 
   # Angles
   Angle    3    1    2  @HW-OW-HW
+
+  Site  'O'
+    Origin  1
+  EndSite
 EndSpecies
 
 #------------------------------------------------------------------------------#

benchmark/siterdf/siterdf_benchmark.cpp

@@ -0,0 +1,28 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+
+#include "benchmark/benchmark.h"
+#include "common/problems.h"
+#include "module/context.h"
+#include "modules/siteRDF/siteRDF.h"
+
+template <ProblemType problem, Population population>
+static void BM_Module_SiteRDF(benchmark::State &state)
+{
+    Problem<problem, population> problemDef;
+    SiteRDFModule module;
+    std::vector<const SpeciesSite *> sites;
+    sites.push_back(problemDef.dissolve_.coreData().species().front()->sites().front().get());
+    module.keywords().set("Configuration", problemDef.cfg_);
+    module.keywords().set("SiteA", sites);
+    module.keywords().set("SiteB", sites);
+    ModuleContext context(problemDef.dissolve_.worldPool(), problemDef.dissolve_);
+    for (auto _ : state)
+    {
+        problemDef.cfg_->incrementContentsVersion();
+        module.executeProcessing(context);
+    }
+}
+
+BENCHMARK_TEMPLATE(BM_Module_SiteRDF, ProblemType::smallMolecule, Population::small)->Unit(benchmark::kMillisecond);;
+BENCHMARK_TEMPLATE(BM_Module_SiteRDF, ProblemType::atomic, Population::large)->Unit(benchmark::kMillisecond);;

src/classes/productIterator.h

@@ -0,0 +1,93 @@
+// SPDX-License-Identifier: GPL-3.0-or-later
+// Copyright (c) 2024 Team Dissolve and contributors
+#pragma once
+
+#include <iterator>
+#include <tuple>
+
+template <typename A, typename B> class ProductIterator
+{
+    private:
+    A x_, xBegin_, xEnd_;
+    B y_, yBegin_, yEnd_;
+    std::decay_t<decltype(yEnd_ - yBegin_)> ySize_;
+
+    public:
+    using valueA = std::decay_t<decltype(*x_)>;
+    using valueB = std::decay_t<decltype(*y_)>;
+    using value = std::tuple<valueA, valueB>;
+    using value_type = std::tuple<valueA, valueB>;
+    using difference_type = int;
+    using reference = std::tuple<const valueA &, const valueB &>;
+    using pointer = std::tuple<valueA, valueB> *;
+    using iterator_category = std::random_access_iterator_tag;
+
+    public:
+    int toIndex() const { return (x_ - xBegin_) * ySize_ + (y_ - yBegin_); }
+    void fromIndex(int index)
+    {
+        x_ = xBegin_ + (index / ySize_);
+        y_ = yBegin_ + (index % ySize_);
+    }
+
+    public:
+    ProductIterator(A xBegin, A xEnd, B yBegin, B yEnd, int index = 0)
+        : xBegin_(xBegin), xEnd_(xEnd), yBegin_(yBegin), yEnd_(yEnd), ySize_(yEnd_ - yBegin_), x_(xBegin), y_(yBegin)
+    {
+        if (index != 0)
+        {
+            fromIndex(index);
+        }
+    }
+
+    ProductIterator begin() const { return ProductIterator(xBegin_, xEnd_, yBegin_, yEnd_); }
+    ProductIterator end() const
+    {
+        ProductIterator it(xBegin_, xEnd_, yBegin_, yEnd_);
+        it.x_ = xEnd_;
+        it.y_ = yBegin_;
+        return it;
+    }
+    reference operator*() { return std::tie(*x_, *y_); }
+    reference operator[](difference_type i) const { return *(*this + i); }
+
+    difference_type operator-(const ProductIterator &it) const { return toIndex() - it.toIndex(); }
+    ProductIterator &operator+=(difference_type forward)
+    {
+        fromIndex(forward + toIndex());
+        return *this;
+    }
+
+    ProductIterator &operator-=(difference_type backward)
+    {
+        fromIndex(toIndex() - backward);
+        return *this;
+    }
+
+    // Operators : arithmetic
+    ProductIterator &operator++()
+    {
+        ++y_;
+        if (y_ >= yEnd_)
+        {
+            ++x_;
+            y_ = yBegin_;
+        }
+        return *this;
+    }
+
+    ProductIterator operator+(difference_type forward) const
+    {
+        return ProductIterator(xBegin_, xEnd_, yBegin_, yEnd_, forward + toIndex());
+    }
+
+    // Operators : comparison
+    bool operator==(const ProductIterator &other) const { return x_ == other.x_ && y_ == other.y_; }
+    bool operator!=(const ProductIterator &other) const { return x_ != other.x_ || y_ != other.y_; }
+    bool operator<(const ProductIterator &other) const
+    {
+        if (x_ != other.x_)
+            return x_ < other.x_;
+        return y_ < other.y_;
+    }
+};

src/keywords/store.cpp

@@ -16,6 +16,7 @@
 #include "keywords/nodeVector.h"
 #include "keywords/range.h"
 #include "keywords/species.h"
+#include "keywords/speciesSiteVector.h"
 #include "keywords/stdString.h"
 #include "keywords/vec3Double.h"
 #include "keywords/vec3Integer.h"
@@ -189,6 +190,11 @@ bool KeywordStore::set(std::string_view name, const std::shared_ptr<Collect1DPro
 {
     return getKeyword<NodeKeyword<Collect1DProcedureNode>>(name, find(name))->setData(value);
 }
+bool KeywordStore::set(std::string_view name, const std::vector<const SpeciesSite *> value)
+{
+    getKeyword<SpeciesSiteVectorKeyword>(name, find(name))->data() = value;
+    return true;
+}
 bool KeywordStore::set(std::string_view name, const std::vector<std::shared_ptr<const Collect1DProcedureNode>> value)
 {
     return getKeyword<NodeVectorKeyword<Collect1DProcedureNode>>(name, find(name))->setData(value);

src/keywords/store.h

@@ -172,6 +172,7 @@ class KeywordStore
     bool set(std::string_view name, const Vec3<NodeValue> value);
     bool set(std::string_view name, const Range value);
     bool set(std::string_view name, const std::shared_ptr<Collect1DProcedureNode> value);
+    bool set(std::string_view name, const std::vector<const SpeciesSite *> value);
     bool set(std::string_view name, const std::vector<std::shared_ptr<const Collect1DProcedureNode>> value);
     bool set(std::string_view name, const std::shared_ptr<RegionProcedureNodeBase> value);
     bool set(std::string_view name, const std::shared_ptr<SelectProcedureNode> value);

src/modules/siteRDF/process.cpp

@@ -3,6 +3,7 @@
 
 #include "analyser/dataExporter.h"
 #include "analyser/dataOperator1D.h"
+#include "classes/productIterator.h"
 #include "base/sysFunc.h"
 #include "io/export/data1D.h"
 #include "main/dissolve.h"
@@ -11,6 +12,8 @@
 #include "math/sampledDouble.h"
 #include "module/context.h"
 #include "modules/siteRDF/siteRDF.h"
+#include "templates/algorithms.h"
+#include "templates/combinable.h"
 
 // Run main processing
 Module::ExecutionResult SiteRDFModule::process(ModuleContext &moduleContext)
@@ -36,15 +39,25 @@ Module::ExecutionResult SiteRDFModule::process(ModuleContext &moduleContext)
         histAB.initialise(distanceRange_.x, distanceRange_.y, distanceRange_.z);
     histAB.zeroBins();
 
-    for (const auto &[siteA, indexA] : a.sites())
-    {
+    auto combinableHistograms = dissolve::CombinableValue<Histogram1D>(
+        [&histAB]()
+        {
+            return histAB;
+        });
+
+    dissolve::for_each(std::execution::seq, a.sites().begin(), a.sites().end(), [this,  &b, &combinableHistograms](const auto& pair) {
+        const auto &[siteA, indexA] = pair;
+
+        auto &hist = combinableHistograms.local();
         for (const auto &[siteB, indexB] : b.sites())
         {
             if (excludeSameMolecule_ && (siteB->molecule() == siteA->molecule()))
-                continue;
-            histAB.bin(targetConfiguration_->box()->minimumDistance(siteA->origin(), siteB->origin()));
+              continue;
+            hist.bin(targetConfiguration_->box()->minimumDistance(siteA->origin(), siteB->origin()));
         }
-    }
+      });
+
+    histAB = combinableHistograms.finalize();
 
     // Accumulate histogram
     histAB.accumulate();