From ca1baffef458e2b4462aab802ab86c22f2a64ef4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 27 Sep 2026 17:11:38 +0200 Subject: [PATCH] RotationScaleMergeGPU, HotPixelFinderGPU: allocate synchronously, like BeamCenterFFTGPU Both queue their work (kernels, cudaMemcpy) on the legacy NULL stream - HotPixelFinderGPU for its download, its sums otherwise on the workers' streams - while their buffers came from the stream-ordered pool, whose cudaFreeAsync is ordered on the thread's non-blocking allocation stream and so after none of that work. Unlike BeamCenterFFTGPU no free has overtaken a read: every entry point waits for its work on the host (cudaDeviceSynchronize, or a blocking device-to-host copy) before it returns. But that holds only by convention, not for a free while unwinding from a failed call, and compute-sanitizer --track-stream-ordered-races cannot see host synchronisation, so it reported every reassigned merge buffer as a use-after-free - noise that buries a real race like the BeamCenterFFTGPU one. compute-sanitizer --tool memcheck --track-stream-ordered-races all: - rugnux --mode scale on a myob _process.h5: 22 use-after-free reports before, 0 after (MergeAccum/MergeAccumRange/MergeRmeas buffers freed by MergeAccum's reassignment or ~Impl). - rugnux -e 150 -N 8 on myob: 92 before, 0 after (together with the next commit). myob, cytc, lyso: p.hkl, p.mtz, p_P1.mtz, p_unmerged.mtz byte-identical to r4-integration. Full-run wall time unchanged within noise (17.4/25.9/18.6 s before, 17.6/25.8/18.7 s after); the scale/merge phase alone (--mode scale, myob) 2.23 -> 2.31 s. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- .../scale_merge/RotationScaleMergeGPU.cu | 109 ++++++++++-------- rugnux/HotPixelsGPU.cu | 14 ++- 2 files changed, 70 insertions(+), 53 deletions(-) diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index b4fcc4352..d5da6e73f 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -16,6 +16,15 @@ namespace { constexpr int BLK = 256; constexpr int MIN_REFLECTIONS = 20; + // Every kernel and copy here is queued on the legacy NULL stream, so - as in BeamCenterFFTGPU - + // no buffer comes from the pool. A pooled buffer is freed with cudaFreeAsync on the thread's + // non-blocking allocation stream, which is not ordered after the NULL stream. Each entry point + // below waits for its own work before it returns, so no free here has yet overtaken a read, but + // that holds only by that convention, and not at all for a free while unwinding from a failed + // call; nor can compute-sanitizer --track-stream-ordered-races see it, and it reported every + // reassigned merge buffer as a use-after-free. cudaFree synchronises the device first. + constexpr CudaAlloc ALLOC = CudaAlloc::Synchronous; + __device__ __forceinline__ double SafeInvD(double x, double fallback) { return (isfinite(x) && x != 0.0) ? 1.0 / x : fallback; } @@ -633,7 +642,7 @@ namespace { template void Upload(CudaDevicePtr &dst, const T *src, int n) { - dst = CudaDevicePtr(std::max(1, n)); + dst = CudaDevicePtr(std::max(1, n), ALLOC); if (n > 0) CudaCheck(cudaMemcpy(dst.get(), src, size_t(n) * sizeof(T), cudaMemcpyHostToDevice), "upload"); } @@ -757,22 +766,22 @@ void RotationScaleMergeGPU::SetPartialsLayout(int n_obs, int n_frames, d.n_frames = n_frames; Upload(d.frame_start, frame_start, n_frames); Upload(d.frame_count, frame_count, n_frames); const int n = std::max(1, n_obs); - d.I = CudaDevicePtr(n); d.sigma = CudaDevicePtr(n); - d.prescaling_corr = CudaDevicePtr(n); d.partiality = CudaDevicePtr(n); - d.zeta = CudaDevicePtr(n); d.corr = CudaDevicePtr(n); - d.bkg = CudaDevicePtr(n); d.var_bkg = CudaDevicePtr(n); - d.image_number = CudaDevicePtr(n); d.d_obs = CudaDevicePtr(n); - d.px_obs = CudaDevicePtr(n); d.py_obs = CudaDevicePtr(n); - d.frame = CudaDevicePtr(n); - d.on_ice = CudaDevicePtr(n); - d.clipped = CudaDevicePtr(n); - d.g = CudaDevicePtr(n_frames); - d.scaled = CudaDevicePtr(n_frames); - d.inv_sigma = CudaDevicePtr(n_obs); - d.sco_coeff = CudaDevicePtr(n_obs); - d.sco_ok = CudaDevicePtr(n_obs); - d.cc = CudaDevicePtr(std::max(1, n_frames)); - d.cc_n = CudaDevicePtr(std::max(1, n_frames)); + d.I = CudaDevicePtr(n, ALLOC); d.sigma = CudaDevicePtr(n, ALLOC); + d.prescaling_corr = CudaDevicePtr(n, ALLOC); d.partiality = CudaDevicePtr(n, ALLOC); + d.zeta = CudaDevicePtr(n, ALLOC); d.corr = CudaDevicePtr(n, ALLOC); + d.bkg = CudaDevicePtr(n, ALLOC); d.var_bkg = CudaDevicePtr(n, ALLOC); + d.image_number = CudaDevicePtr(n, ALLOC); d.d_obs = CudaDevicePtr(n, ALLOC); + d.px_obs = CudaDevicePtr(n, ALLOC); d.py_obs = CudaDevicePtr(n, ALLOC); + d.frame = CudaDevicePtr(n, ALLOC); + d.on_ice = CudaDevicePtr(n, ALLOC); + d.clipped = CudaDevicePtr(n, ALLOC); + d.g = CudaDevicePtr(n_frames, ALLOC); + d.scaled = CudaDevicePtr(n_frames, ALLOC); + d.inv_sigma = CudaDevicePtr(n_obs, ALLOC); + d.sco_coeff = CudaDevicePtr(n_obs, ALLOC); + d.sco_ok = CudaDevicePtr(n_obs, ALLOC); + d.cc = CudaDevicePtr(std::max(1, n_frames), ALLOC); + d.cc_n = CudaDevicePtr(std::max(1, n_frames), ALLOC); } namespace { @@ -828,7 +837,7 @@ void RotationScaleMergeGPU::SetGroups(int n_groups, const int32_t *group, const Upload(d.group_perm, group_perm, n_group_perm); // obs with group >= 0, in group order Upload(d.group_start, group_start, n_groups); Upload(d.group_count, group_count, n_groups); - d.group_mean = CudaDevicePtr(std::max(1, n_groups)); + d.group_mean = CudaDevicePtr(std::max(1, n_groups), ALLOC); } void RotationScaleMergeGPU::SetCorr(const float *corr) { @@ -895,14 +904,14 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit auto &d = *impl_; const int ng = d.n_groups, nf = d.n_fulls; d.merge_for_search = for_search ? 1 : 0; d.merge_min_part = min_partiality; - d.m_sw = CudaDevicePtr(std::max(1, ng)); d.m_swI = CudaDevicePtr(std::max(1, ng)); - d.m_em_mean = CudaDevicePtr(std::max(1, ng)); d.m_cnt = CudaDevicePtr(std::max(1, ng)); + d.m_sw = CudaDevicePtr(std::max(1, ng), ALLOC); d.m_swI = CudaDevicePtr(std::max(1, ng), ALLOC); + d.m_em_mean = CudaDevicePtr(std::max(1, ng), ALLOC); d.m_cnt = CudaDevicePtr(std::max(1, ng), ALLOC); // The error model is fitted on the Bijvoet hands (see the host obs_hand block). Absent - a merge // that already separates them - the hand sums are not allocated and the kernels take the pooled // group, which is what they did before. if (hand && has_hands) { - d.m_swh = CudaDevicePtr(std::max(1, 2 * ng)); d.m_swIh = CudaDevicePtr(std::max(1, 2 * ng)); - d.m_cnth = CudaDevicePtr(std::max(1, 2 * ng)); + d.m_swh = CudaDevicePtr(std::max(1, 2 * ng), ALLOC); d.m_swIh = CudaDevicePtr(std::max(1, 2 * ng), ALLOC); + d.m_cnth = CudaDevicePtr(std::max(1, 2 * ng), ALLOC); if (nf > 0) Upload(d.m_hand, hand, nf); if (ng > 0) Upload(d.m_has_hands, has_hands, ng); } else { @@ -910,8 +919,8 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit d.m_cnth = CudaDevicePtr(); d.m_hand = CudaDevicePtr(); d.m_has_hands = CudaDevicePtr(); } - d.m_s2 = CudaDevicePtr(std::max(1, nf)); d.m_I2 = CudaDevicePtr(std::max(1, nf)); - d.m_dev2 = CudaDevicePtr(std::max(1, nf)); d.m_valid = CudaDevicePtr(std::max(1, nf)); + d.m_s2 = CudaDevicePtr(std::max(1, nf), ALLOC); d.m_I2 = CudaDevicePtr(std::max(1, nf), ALLOC); + d.m_dev2 = CudaDevicePtr(std::max(1, nf), ALLOC); d.m_valid = CudaDevicePtr(std::max(1, nf), ALLOC); MergeParams p{}; p.n_groups = ng; p.min_partiality = min_partiality; @@ -953,13 +962,13 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int ng = d.n_groups; - d.a_swI = CudaDevicePtr(std::max(1, ng)); d.a_sw = CudaDevicePtr(std::max(1, ng)); - d.a_swIh0 = CudaDevicePtr(std::max(1, ng)); d.a_swIh1 = CudaDevicePtr(std::max(1, ng)); - d.a_swh0 = CudaDevicePtr(std::max(1, ng)); d.a_swh1 = CudaDevicePtr(std::max(1, ng)); - d.a_swht0 = CudaDevicePtr(std::max(1, ng)); d.a_swht1 = CudaDevicePtr(std::max(1, ng)); - d.a_nh0 = CudaDevicePtr(std::max(1, ng)); d.a_nh1 = CudaDevicePtr(std::max(1, ng)); - d.a_d = CudaDevicePtr(std::max(1, ng)); d.a_rejected = CudaDevicePtr(std::max(1, ng)); - d.a_on_ice = CudaDevicePtr(std::max(1, ng)); + d.a_swI = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_sw = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_swIh0 = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_swIh1 = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_swh0 = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_swh1 = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_swht0 = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_swht1 = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_nh0 = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_nh1 = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_d = CudaDevicePtr(std::max(1, ng), ALLOC); d.a_rejected = CudaDevicePtr(std::max(1, ng), ALLOC); + d.a_on_ice = CudaDevicePtr(std::max(1, ng), ALLOC); Upload(d.m_rejected, rejected_obs, d.n_fulls); Upload(d.reject_median, reject_median, ng); if (reject_var_add && ng > 0) Upload(d.reject_var_add, reject_var_add, ng); @@ -1019,10 +1028,10 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d auto &d = *impl_; const int ng = d.n_groups; Upload(d.merged_I, merged_I, ng); - d.r_absdev = CudaDevicePtr(std::max(1, ng)); d.r_sumI = CudaDevicePtr(std::max(1, ng)); - d.r_wabsdev = CudaDevicePtr(std::max(1, ng)); d.r_wsumI = CudaDevicePtr(std::max(1, ng)); - d.r_sumv = CudaDevicePtr(std::max(1, ng)); d.r_sumv2 = CudaDevicePtr(std::max(1, ng)); - d.r_n = CudaDevicePtr(std::max(1, ng)); d.r_nusable = CudaDevicePtr(std::max(1, ng)); + d.r_absdev = CudaDevicePtr(std::max(1, ng), ALLOC); d.r_sumI = CudaDevicePtr(std::max(1, ng), ALLOC); + d.r_wabsdev = CudaDevicePtr(std::max(1, ng), ALLOC); d.r_wsumI = CudaDevicePtr(std::max(1, ng), ALLOC); + d.r_sumv = CudaDevicePtr(std::max(1, ng), ALLOC); d.r_sumv2 = CudaDevicePtr(std::max(1, ng), ALLOC); + d.r_n = CudaDevicePtr(std::max(1, ng), ALLOC); d.r_nusable = CudaDevicePtr(std::max(1, ng), ALLOC); MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; @@ -1079,7 +1088,7 @@ void RotationScaleMergeGPU::SmoothFullsCorr(const uint8_t *apply, const double * int64_t RotationScaleMergeGPU::FilterCorrByZeta(double min_zeta) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; - CudaDevicePtr dropped(1); + CudaDevicePtr dropped(1, ALLOC); CudaCheck(cudaMemset(dropped.get(), 0, sizeof(unsigned long long)), "zero zeta drop count"); const int blocks = std::min(65535, (d.n_obs + BLK - 1) / BLK); FilterZetaKernel<<>>(d.n_obs, min_zeta, d.zeta.get(), d.corr.get(), dropped.get()); @@ -1136,9 +1145,9 @@ void RotationScaleMergeGPU::SetRawRuns(int n_runs, int n_perm, const int32_t *pe Upload(d.rr_k, rr_k, n_runs); Upload(d.rr_l, rr_l, n_runs); - d.rr_group = CudaDevicePtr(std::max(1, n_runs)); - d.rr_nevents = CudaDevicePtr(std::max(1, n_runs)); - d.rr_offset = CudaDevicePtr(std::max(1, n_runs)); + d.rr_group = CudaDevicePtr(std::max(1, n_runs), ALLOC); + d.rr_nevents = CudaDevicePtr(std::max(1, n_runs), ALLOC); + d.rr_offset = CudaDevicePtr(std::max(1, n_runs), ALLOC); } int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_partiality, @@ -1180,17 +1189,17 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti // Allocate the fulls SoA and emit. const int nf = std::max(1, d.n_fulls); - d.f_h = CudaDevicePtr(nf); d.f_k = CudaDevicePtr(nf); d.f_l = CudaDevicePtr(nf); - d.f_frame = CudaDevicePtr(nf); d.f_group = CudaDevicePtr(nf); - d.f_I = CudaDevicePtr(nf); d.f_sigma = CudaDevicePtr(nf); - d.f_d = CudaDevicePtr(nf); d.f_img = CudaDevicePtr(nf); - d.f_px = CudaDevicePtr(nf); d.f_py = CudaDevicePtr(nf); - d.f_var_bkg = CudaDevicePtr(nf); d.f_var_per_I = CudaDevicePtr(nf); - d.f_on_ice = CudaDevicePtr(nf); d.f_clipped = CudaDevicePtr(nf); - d.f_corr = CudaDevicePtr(nf); d.f_partiality = CudaDevicePtr(nf); - d.f_rlp = CudaDevicePtr(nf); d.f_zeta = CudaDevicePtr(nf); - d.f_inv_sigma = CudaDevicePtr(nf); - d.f_sco_coeff = CudaDevicePtr(nf); d.f_sco_ok = CudaDevicePtr(nf); + d.f_h = CudaDevicePtr(nf, ALLOC); d.f_k = CudaDevicePtr(nf, ALLOC); d.f_l = CudaDevicePtr(nf, ALLOC); + d.f_frame = CudaDevicePtr(nf, ALLOC); d.f_group = CudaDevicePtr(nf, ALLOC); + d.f_I = CudaDevicePtr(nf, ALLOC); d.f_sigma = CudaDevicePtr(nf, ALLOC); + d.f_d = CudaDevicePtr(nf, ALLOC); d.f_img = CudaDevicePtr(nf, ALLOC); + d.f_px = CudaDevicePtr(nf, ALLOC); d.f_py = CudaDevicePtr(nf, ALLOC); + d.f_var_bkg = CudaDevicePtr(nf, ALLOC); d.f_var_per_I = CudaDevicePtr(nf, ALLOC); + d.f_on_ice = CudaDevicePtr(nf, ALLOC); d.f_clipped = CudaDevicePtr(nf, ALLOC); + d.f_corr = CudaDevicePtr(nf, ALLOC); d.f_partiality = CudaDevicePtr(nf, ALLOC); + d.f_rlp = CudaDevicePtr(nf, ALLOC); d.f_zeta = CudaDevicePtr(nf, ALLOC); + d.f_inv_sigma = CudaDevicePtr(nf, ALLOC); + d.f_sco_coeff = CudaDevicePtr(nf, ALLOC); d.f_sco_ok = CudaDevicePtr(nf, ALLOC); CudaCheck(cudaMemcpy(d.rr_offset.get(), offset.data(), size_t(d.n_runs) * sizeof(int32_t), cudaMemcpyHostToDevice), "upload offset"); diff --git a/rugnux/HotPixelsGPU.cu b/rugnux/HotPixelsGPU.cu index b17ecb098..62a3d22b6 100644 --- a/rugnux/HotPixelsGPU.cu +++ b/rugnux/HotPixelsGPU.cu @@ -145,15 +145,23 @@ __global__ void accumulate_kernel(const int32_t *__restrict__ image, const int32 } } +// The shared tables and sums are filled on a stream of their own, added to on the workers' streams +// and downloaded on the NULL stream, so they are allocated synchronously rather than from the pool: +// a pooled buffer is freed on the thread's allocation stream, which none of those is ordered before. +// Each of those steps is waited for on the host, so this costs nothing but the device-wide sync of +// eight allocations per run, and it leaves compute-sanitizer --track-stream-ordered-races nothing +// to report here, where it flagged every one of them. +constexpr CudaAlloc ALLOC = CudaAlloc::Synchronous; + } // namespace HotPixelFinderGPU::HotPixelFinderGPU(const int32_t *host_key, size_t npixels, const std::vector &host_key_begin, int nrings, int sectors) : npixels(npixels), nkeys(static_cast(nrings) * sectors), nrings(nrings), sectors(sectors), - key(npixels), pixels_by_key(host_key_begin.back()), key_begin(host_key_begin.size()), - n_lit(npixels), n_error(npixels), n_error_ring_ok(npixels), - sum_value(npixels), error_level_sum(npixels) { + key(npixels, ALLOC), pixels_by_key(host_key_begin.back(), ALLOC), key_begin(host_key_begin.size(), ALLOC), + n_lit(npixels, ALLOC), n_error(npixels, ALLOC), n_error_ring_ok(npixels, ALLOC), + sum_value(npixels, ALLOC), error_level_sum(npixels, ALLOC) { std::vector pixels(host_key_begin.back()); std::vector filled(host_key_begin.begin(), host_key_begin.end() - 1); for (size_t i = 0; i < npixels; i++)