diff --git a/common/CUDAWrapper.cpp b/common/CUDAWrapper.cpp index ab3e2e6c..ece8f4be 100644 --- a/common/CUDAWrapper.cpp +++ b/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 diff --git a/common/CUDAWrapper.cu b/common/CUDAWrapper.cu index 554c54c6..40ba2e37 100644 --- a/common/CUDAWrapper.cu +++ b/common/CUDAWrapper.cu @@ -1,8 +1,16 @@ // SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#include +#include #include +#include +#include +#include +#include +#ifdef __linux__ +#include +#include +#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(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(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 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(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 } diff --git a/common/CUDAWrapper.h b/common/CUDAWrapper.h index 4b288ee0..f29de57c 100644 --- a/common/CUDAWrapper.h +++ b/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(); diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 838a4a45..e00bbeac 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -116,7 +116,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(outcomes.size()); std::vector pts_offset(n_out + 1, 0); - ParallelChunks(n_out, nthreads, [&](int lo, int hi) { + ParallelChunks(n_out, ThreadsForWork(static_cast(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 pts(pts_offset[n_out]); - ParallelChunks(n_out, nthreads, [&](int lo, int hi) { + ParallelChunks(n_out, ThreadsForWork(static_cast(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 terms(scale_events.size()); std::vector 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(scale_events.size()), nthreads, [&](int lo, int hi) { + ParallelChunks(static_cast(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 total{}; std::mutex mx; - ParallelChunks(static_cast(terms.size()), nthreads, [&](int lo, int hi) { + ParallelChunks(static_cast(terms.size()), ThreadsForWork(terms.size(), nthreads), [&](int lo, int hi) { std::array acc{}; for (int e = lo; e < hi; ++e) { const double r = residual_at(terms[e], k); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7cbc6aae..2ce92ebc 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/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(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(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(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 h_of(partials.size()); int h_min = 0, h_max = 0; { const int n_obs = static_cast(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(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(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 corr(partials.size()); gpu_->GetCorr(corr.data()); - for (size_t i = 0; i < partials.size(); ++i) partials[i].corr = corr[i]; + ParallelChunks(static_cast(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 fcorr(fulls.size()); - for (size_t i = 0; i < fulls.size(); ++i) fcorr[i] = fulls[i].corr; + ParallelChunks(static_cast(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 diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 75d3e4d3..343d56ba 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1662,10 +1662,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b std::function worker = per_image_analysis ? std::function(full_worker) : std::function(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 > 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();