Stop paying for workers and threads that do no work

Four changes, all from measuring why this 48-core machine was SLOWER than a
quarter of itself on a single job.

The image loop now takes a worker count of its own. It is GPU-bound, every
worker builds a private analysis engine of tens of megabytes of device and
pinned memory - and every one of those allocations implicitly synchronises the
device - so past a handful per card another worker adds setup and contention and
no throughput. Measured at about 23 ms of pure setup per extra worker, which is
why the penalty is WORSE on short runs: 200 images cost 0.86 s of loop at 12
workers and 2.25 s at 48. Capped at four per GPU, floor of eight. Every other
phase still gets the full thread count, because each one starts its own workers.

Ingest built its array with a serial push_back over every observation of every
frame - 63 million of them on the largest crystal here. Each frame's block offset
is known before anything is written, so the frames convert together, each still
written by one thread in its own order.

The pass that buckets observations by h was the single most expensive thing in a
large run - 24% of all cycles, in five instructions. It strided an array of
80-byte observations to read one 4-byte field, and its store address depended on
the loaded value, so the store buffer could not retire and the misses stopped
overlapping. The sweep that already reads every observation now copies h out as
it goes, and the bucketing walks that instead.

And the post-refine passes took the raw thread count. One of them runs 134 times
inside the rotation-scale fit, starting 48 threads each time to divide 390k terms
among them; it is gated on the work now, like everything else.

Measured on one crystal, N=48: 12.80 s -> 11.45 s, which is what 12 threads used
to cost, and on the best-matched pair the two are now level. On the heaviest
crystal ingest goes 13.5 s -> 9.9 s and the run 77.8 s -> 69.4 s. Battery 9m01s
-> 8m24s, space group 21/24, no failures.

Also restores get_gpu_numa_node() - the sysfs lookup deleted with NUMAHWPolicy -
and an opt-in CPU pin to that node behind JFJOCH_PIN_CPU_TO_GPU_NODE. It is off
because it measured neutral here: all four GPUs hang off two of the four nodes,
so pinning to them costs a worker the other half of the machine. It is kept for
boxes whose GPUs are spread over every socket. The old lookup had a latent bug -
CUDA reports the PCI id with upper-case hex and sysfs paths are lower case, so on
three of this machine's four GPUs it would have silently returned "unknown".

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

common/CUDAWrapper.cpp

@@ -11,6 +11,8 @@ int32_t get_gpu_count() {
 
 void set_gpu(int32_t dev_id) {}
 
+int get_gpu_numa_node(int32_t dev_id) { return -1; }
+
 void pin_gpu() {}
 
 #endif

common/CUDAWrapper.cu

@@ -1,8 +1,16 @@
 // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
-#include <fstream>
+#include <algorithm>
 #include <atomic>
+#include <cctype>
+#include <cstdlib>
+#include <fstream>
+#include <string>
+#ifdef __linux__
+#include <pthread.h>
+#include <sched.h>
+#endif
 
 #include "CUDAWrapper.h"
 #include "JFJochException.h"
@@ -39,9 +47,92 @@ void set_gpu(int32_t dev_id) {
     }
 }
 
+// The NUMA node a GPU is attached to, from sysfs; -1 when unknown or not applicable. No libnuma:
+// this is a file read. (This lookup existed before cc925b26 and is restored unchanged apart from
+// lower-casing the bus id - CUDA reports hex digits in upper case on some drivers and sysfs paths
+// are lower case.)
+int get_gpu_numa_node(int32_t dev_id) {
+#ifdef __linux__
+    if (dev_id < 0 || dev_id >= get_gpu_count())
+        return -1;
+    char buf[64] = {};
+    if (cudaDeviceGetPCIBusId(buf, static_cast<int>(sizeof(buf)), dev_id) != cudaSuccess)
+        return -1;
+    std::string bus_id(buf);
+    std::transform(bus_id.begin(), bus_id.end(), bus_id.begin(),
+                   [](unsigned char c) { return static_cast<char>(std::tolower(c)); });
+    std::ifstream f("/sys/bus/pci/devices/" + bus_id + "/numa_node");
+    int node = -1;
+    if (!(f >> node))
+        return -1;
+    return node;
+#else
+    (void) dev_id;
+    return -1;
+#endif
+}
+
+#ifdef __linux__
+namespace {
+    // Confine the calling thread to the cores of one NUMA node, read from sysfs. Its memory then
+    // follows by first touch, which is what actually matters - the placement of what the thread
+    // allocates, not the affinity itself.
+    void run_on_numa_node(int node) {
+        std::ifstream f("/sys/devices/system/node/node" + std::to_string(node) + "/cpulist");
+        std::string list;
+        if (!std::getline(f, list) || list.empty())
+            return;
+        cpu_set_t set;
+        CPU_ZERO(&set);
+        // cpulist is comma-separated singles and a-b ranges, e.g. "24-35" or "0,2,4-7".
+        size_t pos = 0;
+        while (pos < list.size()) {
+            size_t comma = list.find(',', pos);
+            const std::string item = list.substr(pos, comma == std::string::npos ? comma : comma - pos);
+            const size_t dash = item.find('-');
+            const int lo = std::atoi(item.c_str());
+            const int hi = dash == std::string::npos ? lo : std::atoi(item.c_str() + dash + 1);
+            for (int c = lo; c <= hi && c < CPU_SETSIZE; ++c)
+                CPU_SET(c, &set);
+            if (comma == std::string::npos) break;
+            pos = comma + 1;
+        }
+        if (CPU_COUNT(&set) > 0)
+            pthread_setaffinity_np(pthread_self(), sizeof(set), &set);
+    }
+
+    int numa_node_count() {
+        int n = 0;
+        while (std::ifstream("/sys/devices/system/node/node" + std::to_string(n) + "/cpulist").good())
+            ++n;
+        return n;
+    }
+}
+#endif
+
 void pin_gpu() {
     static std::atomic<uint32_t> counter{0};
     auto dev_count = get_gpu_count();
-    if (dev_count > 0)
-        set_gpu(counter.fetch_add(1) % dev_count);
+    if (dev_count <= 0)
+        return;
+    const int32_t dev = static_cast<int32_t>(counter.fetch_add(1) % dev_count);
+    set_gpu(dev);
+
+#ifdef __linux__
+    // Optionally also confine the thread to the cores local to that GPU. On a machine whose GPUs all
+    // hang off some of the sockets, a worker can otherwise sit on a socket with no GPU of its own and
+    // cross the interconnect for every transfer and every page it allocates. Off by default: it is a
+    // real trade - it also denies the thread the cores of the other sockets - and only measurement on
+    // a given box can say which way it goes. No-op unless there is more than one node and the GPU's
+    // node is known.
+    static const bool pin_cpu = [] {
+        const char *e = std::getenv("JFJOCH_PIN_CPU_TO_GPU_NODE");
+        return e && *e && *e != '0';
+    }();
+    if (!pin_cpu || numa_node_count() < 2)
+        return;
+    const int node = get_gpu_numa_node(dev);
+    if (node >= 0)
+        run_on_numa_node(node);
+#endif
 }

common/CUDAWrapper.h

@@ -8,7 +8,16 @@
 int32_t get_gpu_count();
 void set_gpu(int32_t dev_id);
 
+// NUMA node the given GPU is attached to, read from sysfs; -1 when unknown, when there is no such
+// device, or on a platform where the question does not apply.
+int get_gpu_numa_node(int32_t dev_id);
+
 // Pin the calling thread to the next GPU in round-robin order, using a process-wide counter
 // (counter++ % get_gpu_count()). Call once per thread; no thread id needed. No-op when no GPU
 // is visible. Honours CUDA_VISIBLE_DEVICES via get_gpu_count().
+//
+// With JFJOCH_PIN_CPU_TO_GPU_NODE set in the environment it also confines the thread to the cores of
+// that GPU's NUMA node, so the memory it allocates lands beside the GPU it feeds. Off by default -
+// it also takes the other sockets' cores away from the thread, and which way that trade goes depends
+// on the machine.
 void pin_gpu();

image_analysis/geom_refinement/PostRefine.cpp

@@ -116,7 +116,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
         const size_t nthreads = std::max(1, settings.num_threads);
         const int n_out = static_cast<int>(outcomes.size());
         std::vector<size_t> pts_offset(n_out + 1, 0);
-        ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
+        ParallelChunks(n_out, ThreadsForWork(static_cast<size_t>(n_out), nthreads), [&](int lo, int hi) {
             for (int o = lo; o < hi; o++) {
                 size_t keep = 0;
                 for (const auto &r : outcomes[o].reflections)
@@ -129,7 +129,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             pts_offset[o + 1] += pts_offset[o];
 
         std::vector<Partial> pts(pts_offset[n_out]);
-        ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
+        ParallelChunks(n_out, ThreadsForWork(static_cast<size_t>(n_out), nthreads), [&](int lo, int hi) {
             for (int o = lo; o < hi; o++) {
                 size_t at = pts_offset[o];
                 for (const auto &r : outcomes[o].reflections) {
@@ -364,7 +364,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             std::vector<ScaleTerm> terms(scale_events.size());
             std::vector<int> fifth_of(scale_events.size());
             const double aa_c[3] = {-phi_c * u[0], -phi_c * u[1], -phi_c * u[2]};
-            ParallelChunks(static_cast<int>(scale_events.size()), nthreads, [&](int lo, int hi) {
+            ParallelChunks(static_cast<int>(scale_events.size()), ThreadsForWork(scale_events.size(), nthreads), [&](int lo, int hi) {
                 for (int e = lo; e < hi; ++e) {
                     const double p[3] = {scale_events[e].e_ref[0] / s, scale_events[e].e_ref[1] / s,
                                          scale_events[e].e_ref[2] / s};
@@ -405,7 +405,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             const auto cost_by_fifth = [&](double k) {
                 std::array<double, 5> total{};
                 std::mutex mx;
-                ParallelChunks(static_cast<int>(terms.size()), nthreads, [&](int lo, int hi) {
+                ParallelChunks(static_cast<int>(terms.size()), ThreadsForWork(terms.size(), nthreads), [&](int lo, int hi) {
                     std::array<double, 5> acc{};
                     for (int e = lo; e < hi; ++e) {
                         const double r = residual_at(terms[e], k);

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -256,14 +256,28 @@ void RotationScaleMerge::Ingest() {
     const float dist_tol = x.GetIndexingSettings().GetUnitCellDistTolerance();
     const float ang_tol = x.GetIndexingSettings().GetUnitCellAngleTolerance_deg();
 
-    for (int o = 0; o < n_frames; ++o) {
-        frame_start[o] = static_cast<int32_t>(partials.size());
+    // A frame's observations occupy one contiguous block and every block's offset is known before a
+    // byte is written, so the frames convert in parallel. This was a serial push_back over every
+    // observation of every frame - 63 million of them on the largest crystal in the test set, and the
+    // single largest serial stretch of its ingest. Each frame is still written by one thread in its
+    // own order, so the result is the same array.
+    {
+        int32_t acc = 0;
+        for (int o = 0; o < n_frames; ++o) {
+            frame_start[o] = acc;
+            frame_count[o] = static_cast<int32_t>(partials_out[o].reflections.size());
+            acc += frame_count[o];
+        }
+    }
+    partials.resize(total);
+    ParallelFor(n_frames, ThreadsForWork(total, nthreads), [&](int o) {
         if (reference_cell) {
             const auto cell = partials_out[o].latt.GetUnitCell();
             frame_cell_ok[o] = cell.is_close(*reference_cell, dist_tol, ang_tol) ? 1 : 0;
         }
+        int32_t at = frame_start[o];
         for (const auto &r : partials_out[o].reflections) {
-            Obs obs{};
+            Obs &obs = partials[at++];
             obs.h = r.h; obs.k = r.k; obs.l = r.l;
             obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; obs.rlp = r.rlp;
             obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg;
@@ -273,10 +287,8 @@ void RotationScaleMerge::Ingest() {
             obs.on_ice = r.on_ice_ring ? 1 : 0;
             obs.corr = r.image_scale_corr;
             obs.group = -1;
-            partials.push_back(obs);
         }
-        frame_count[o] = static_cast<int32_t>(partials.size()) - frame_start[o];
-    }
+    });
 
     DivideOutIncidentFlux();
 
@@ -284,7 +296,13 @@ void RotationScaleMerge::Ingest() {
     // obs from a flat 1-byte array instead of re-reading the fat Obs struct for every space group.
     finite_ok.resize(partials.size());
     // The h range comes out of the same sweep - the sort below buckets by h and needs to know how
-    // many buckets that is, and this pass already reads every observation.
+    // many buckets that is, and this pass already reads every observation. It also copies out h on
+    // its way past, because the bucketing below reads nothing else from the observation: walking an
+    // array of 4-byte h instead of striding 80-byte Obs turns a pass that was the single most
+    // expensive thing in a large run into one that barely registers. Striding hurts there and not
+    // here because the histogram's STORE address depends on the loaded value, so the store buffer
+    // cannot retire and the misses stop overlapping.
+    std::vector<int32_t> h_of(partials.size());
     int h_min = 0, h_max = 0;
     {
         const int n_obs = static_cast<int>(partials.size());
@@ -301,6 +319,7 @@ void RotationScaleMerge::Ingest() {
                     const auto &o = partials[i];
                     finite_ok[i] = (std::isfinite(o.I) && std::isfinite(o.rlp) && o.rlp != 0.0f
                                     && std::isfinite(o.sigma) && o.sigma > 0.0f) ? 1 : 0;
+                    h_of[i] = o.h;
                     lmin = std::min(lmin, o.h);
                     lmax = std::max(lmax, o.h);
                 }
@@ -354,7 +373,7 @@ void RotationScaleMerge::Ingest() {
         ParallelChunks(nt, nthreads, [&](int tlo, int thi) {
             for (int t = tlo; t < thi; ++t) {
                 const int lo = t * chunk, hi = std::min(n, lo + chunk);
-                for (int i = lo; i < hi; ++i) hist[t][partials[i].h - h_min]++;
+                for (int i = lo; i < hi; ++i) hist[t][h_of[i] - h_min]++;
             }
         });
 
@@ -379,7 +398,7 @@ void RotationScaleMerge::Ingest() {
                 const int lo = t * chunk, hi = std::min(n, lo + chunk);
                 for (int i = lo; i < hi; ++i) {
                     const auto &o = partials[i];
-                    keys[fill[o.h - h_min]++] =
+                    keys[fill[h_of[i] - h_min]++] =
                             SortKey{o.h, o.k, o.l, o.image_number, static_cast<int32_t>(i)};
                 }
             }
@@ -3051,8 +3070,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
     // Put back the observations the previous pass filtered out of its own merge (see where this is
     // filled). This pass decides for itself which ones to drop, and the final merge drops none.
     if (!corr_before_pass_filters.empty()) {
-        for (size_t i = 0; i < partials.size(); ++i)
-            partials[i].corr = corr_before_pass_filters[i];
+        ParallelChunks(static_cast<int>(partials.size()),
+                       ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+            for (int i = lo; i < hi; ++i) partials[i].corr = corr_before_pass_filters[i];
+        });
 #ifdef JFJOCH_USE_CUDA
         if (gpu_active_)
             gpu_->SetCorr(corr_before_pass_filters.data());
@@ -3149,7 +3170,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
     if (gpu_active_) {
         std::vector<float> corr(partials.size());
         gpu_->GetCorr(corr.data());
-        for (size_t i = 0; i < partials.size(); ++i) partials[i].corr = corr[i];
+        ParallelChunks(static_cast<int>(partials.size()),
+                       ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+            for (int i = lo; i < hi; ++i) partials[i].corr = corr[i];
+        });
     }
 #endif
 
@@ -3356,7 +3380,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
     // the corrected corr back to the device so the merge reads it.
     if ((corrections || rejected_full_scales) && combined_on_gpu && scaled_fulls_on_gpu) {
         std::vector<float> fcorr(fulls.size());
-        for (size_t i = 0; i < fulls.size(); ++i) fcorr[i] = fulls[i].corr;
+        ParallelChunks(static_cast<int>(fulls.size()),
+                       ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) {
+            for (int i = lo; i < hi; ++i) fcorr[i] = fulls[i].corr;
+        });
         gpu_->SetFullsCorr(fcorr.data());
     }
 #endif

rugnux/Rugnux.cpp

@@ -1662,10 +1662,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
 
     std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
                                                       : std::function<void()>(azint_worker);
+    // Workers for this loop only. It is GPU-bound - every worker builds its own analysis engine, tens
+    // of megabytes of device and pinned memory whose allocation implicitly synchronises the device,
+    // and they all queue on the same few cards - so past a handful per GPU another worker adds setup
+    // and contention and no throughput. Measured on a four-GPU node: the loop is flat from 8 workers
+    // up and slower at 48 than at 12. Every other phase still gets the full thread count, because
+    // each one spawns its own workers.
+    const int loop_workers = get_gpu_count() > 0
+            ? std::min(config_.nthreads, std::max(8, 4 * get_gpu_count()))
+            : config_.nthreads;
     std::vector<std::future<void> > futures;
-    futures.reserve(config_.nthreads);
+    futures.reserve(loop_workers);
     const auto image_loop_start = std::chrono::steady_clock::now();
-    for (int i = 0; i < config_.nthreads; ++i)
+    for (int i = 0; i < loop_workers; ++i)
         futures.push_back(std::async(std::launch::async, worker));
     for (auto &f: futures)
         f.get();