mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 80 Wiki Activity

RotationScaleMergeGPU, HotPixelFinderGPU: allocate synchronously, like BeamCenterFFTGPU

Browse Source 
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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_f
2026-09-27 17:11:38 +02:00
co-authored by Claude Opus 5.5
parent 77710e1cec
commit ca1baffef4
2 changed files with 70 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMergeGPU.cu
+59 -50
View File
@@ -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 <typename T>
void Upload(CudaDevicePtr<T> &dst, const T *src, int n) {
dst = CudaDevicePtr<T>(std::max(1, n));
dst = CudaDevicePtr<T>(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<float>(n); d.sigma = CudaDevicePtr<float>(n);
d.prescaling_corr = CudaDevicePtr<float>(n); d.partiality = CudaDevicePtr<float>(n);
d.zeta = CudaDevicePtr<float>(n); d.corr = CudaDevicePtr<float>(n);
d.bkg = CudaDevicePtr<float>(n); d.var_bkg = CudaDevicePtr<float>(n);
d.image_number = CudaDevicePtr<float>(n); d.d_obs = CudaDevicePtr<float>(n);
d.px_obs = CudaDevicePtr<float>(n); d.py_obs = CudaDevicePtr<float>(n);
d.frame = CudaDevicePtr<int32_t>(n);
d.on_ice = CudaDevicePtr<uint8_t>(n);
d.clipped = CudaDevicePtr<uint8_t>(n);
d.g = CudaDevicePtr<double>(n_frames);
d.scaled = CudaDevicePtr<uint8_t>(n_frames);
d.inv_sigma = CudaDevicePtr<double>(n_obs);
d.sco_coeff = CudaDevicePtr<float>(n_obs);
d.sco_ok = CudaDevicePtr<uint8_t>(n_obs);
d.cc = CudaDevicePtr<double>(std::max(1, n_frames));
d.cc_n = CudaDevicePtr<int64_t>(std::max(1, n_frames));
d.I = CudaDevicePtr<float>(n, ALLOC); d.sigma = CudaDevicePtr<float>(n, ALLOC);
d.prescaling_corr = CudaDevicePtr<float>(n, ALLOC); d.partiality = CudaDevicePtr<float>(n, ALLOC);
d.zeta = CudaDevicePtr<float>(n, ALLOC); d.corr = CudaDevicePtr<float>(n, ALLOC);
d.bkg = CudaDevicePtr<float>(n, ALLOC); d.var_bkg = CudaDevicePtr<float>(n, ALLOC);
d.image_number = CudaDevicePtr<float>(n, ALLOC); d.d_obs = CudaDevicePtr<float>(n, ALLOC);
d.px_obs = CudaDevicePtr<float>(n, ALLOC); d.py_obs = CudaDevicePtr<float>(n, ALLOC);
d.frame = CudaDevicePtr<int32_t>(n, ALLOC);
d.on_ice = CudaDevicePtr<uint8_t>(n, ALLOC);
d.clipped = CudaDevicePtr<uint8_t>(n, ALLOC);
d.g = CudaDevicePtr<double>(n_frames, ALLOC);
d.scaled = CudaDevicePtr<uint8_t>(n_frames, ALLOC);
d.inv_sigma = CudaDevicePtr<double>(n_obs, ALLOC);
d.sco_coeff = CudaDevicePtr<float>(n_obs, ALLOC);
d.sco_ok = CudaDevicePtr<uint8_t>(n_obs, ALLOC);
d.cc = CudaDevicePtr<double>(std::max(1, n_frames), ALLOC);
d.cc_n = CudaDevicePtr<int64_t>(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<double>(std::max(1, n_groups));
d.group_mean = CudaDevicePtr<double>(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<double>(std::max(1, ng)); d.m_swI = CudaDevicePtr<double>(std::max(1, ng));
d.m_em_mean = CudaDevicePtr<double>(std::max(1, ng)); d.m_cnt = CudaDevicePtr<int32_t>(std::max(1, ng));
d.m_sw = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.m_swI = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.m_em_mean = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.m_cnt = CudaDevicePtr<int32_t>(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<double>(std::max(1, 2 * ng)); d.m_swIh = CudaDevicePtr<double>(std::max(1, 2 * ng));
d.m_cnth = CudaDevicePtr<int32_t>(std::max(1, 2 * ng));
d.m_swh = CudaDevicePtr<double>(std::max(1, 2 * ng), ALLOC); d.m_swIh = CudaDevicePtr<double>(std::max(1, 2 * ng), ALLOC);
d.m_cnth = CudaDevicePtr<int32_t>(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<int32_t>();
d.m_hand = CudaDevicePtr<uint8_t>(); d.m_has_hands = CudaDevicePtr<uint8_t>();
}
d.m_s2 = CudaDevicePtr<double>(std::max(1, nf)); d.m_I2 = CudaDevicePtr<double>(std::max(1, nf));
d.m_dev2 = CudaDevicePtr<double>(std::max(1, nf)); d.m_valid = CudaDevicePtr<uint8_t>(std::max(1, nf));
d.m_s2 = CudaDevicePtr<double>(std::max(1, nf), ALLOC); d.m_I2 = CudaDevicePtr<double>(std::max(1, nf), ALLOC);
d.m_dev2 = CudaDevicePtr<double>(std::max(1, nf), ALLOC); d.m_valid = CudaDevicePtr<uint8_t>(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<double>(std::max(1, ng)); d.a_sw = CudaDevicePtr<double>(std::max(1, ng));
d.a_swIh0 = CudaDevicePtr<double>(std::max(1, ng)); d.a_swIh1 = CudaDevicePtr<double>(std::max(1, ng));
d.a_swh0 = CudaDevicePtr<double>(std::max(1, ng)); d.a_swh1 = CudaDevicePtr<double>(std::max(1, ng));
d.a_swht0 = CudaDevicePtr<double>(std::max(1, ng)); d.a_swht1 = CudaDevicePtr<double>(std::max(1, ng));
d.a_nh0 = CudaDevicePtr<int32_t>(std::max(1, ng)); d.a_nh1 = CudaDevicePtr<int32_t>(std::max(1, ng));
d.a_d = CudaDevicePtr<double>(std::max(1, ng)); d.a_rejected = CudaDevicePtr<int32_t>(std::max(1, ng));
d.a_on_ice = CudaDevicePtr<uint8_t>(std::max(1, ng));
d.a_swI = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.a_sw = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.a_swIh0 = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.a_swIh1 = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.a_swh0 = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.a_swh1 = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.a_swht0 = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.a_swht1 = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.a_nh0 = CudaDevicePtr<int32_t>(std::max(1, ng), ALLOC); d.a_nh1 = CudaDevicePtr<int32_t>(std::max(1, ng), ALLOC);
d.a_d = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.a_rejected = CudaDevicePtr<int32_t>(std::max(1, ng), ALLOC);
d.a_on_ice = CudaDevicePtr<uint8_t>(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<double>(std::max(1, ng)); d.r_sumI = CudaDevicePtr<double>(std::max(1, ng));
d.r_wabsdev = CudaDevicePtr<double>(std::max(1, ng)); d.r_wsumI = CudaDevicePtr<double>(std::max(1, ng));
d.r_sumv = CudaDevicePtr<double>(std::max(1, ng)); d.r_sumv2 = CudaDevicePtr<double>(std::max(1, ng));
d.r_n = CudaDevicePtr<int32_t>(std::max(1, ng)); d.r_nusable = CudaDevicePtr<int32_t>(std::max(1, ng));
d.r_absdev = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.r_sumI = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.r_wabsdev = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.r_wsumI = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.r_sumv = CudaDevicePtr<double>(std::max(1, ng), ALLOC); d.r_sumv2 = CudaDevicePtr<double>(std::max(1, ng), ALLOC);
d.r_n = CudaDevicePtr<int32_t>(std::max(1, ng), ALLOC); d.r_nusable = CudaDevicePtr<int32_t>(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<unsigned long long> dropped(1);
CudaDevicePtr<unsigned long long> 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<<<blocks, BLK>>>(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<int32_t>(std::max(1, n_runs));
d.rr_nevents = CudaDevicePtr<int32_t>(std::max(1, n_runs));
d.rr_offset = CudaDevicePtr<int32_t>(std::max(1, n_runs));
d.rr_group = CudaDevicePtr<int32_t>(std::max(1, n_runs), ALLOC);
d.rr_nevents = CudaDevicePtr<int32_t>(std::max(1, n_runs), ALLOC);
d.rr_offset = CudaDevicePtr<int32_t>(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<int32_t>(nf); d.f_k = CudaDevicePtr<int32_t>(nf); d.f_l = CudaDevicePtr<int32_t>(nf);
d.f_frame = CudaDevicePtr<int32_t>(nf); d.f_group = CudaDevicePtr<int32_t>(nf);
d.f_I = CudaDevicePtr<float>(nf); d.f_sigma = CudaDevicePtr<float>(nf);
d.f_d = CudaDevicePtr<float>(nf); d.f_img = CudaDevicePtr<float>(nf);
d.f_px = CudaDevicePtr<float>(nf); d.f_py = CudaDevicePtr<float>(nf);
d.f_var_bkg = CudaDevicePtr<float>(nf); d.f_var_per_I = CudaDevicePtr<float>(nf);
d.f_on_ice = CudaDevicePtr<uint8_t>(nf); d.f_clipped = CudaDevicePtr<uint8_t>(nf);
d.f_corr = CudaDevicePtr<float>(nf); d.f_partiality = CudaDevicePtr<float>(nf);
d.f_rlp = CudaDevicePtr<float>(nf); d.f_zeta = CudaDevicePtr<float>(nf);
d.f_inv_sigma = CudaDevicePtr<double>(nf);
d.f_sco_coeff = CudaDevicePtr<float>(nf); d.f_sco_ok = CudaDevicePtr<uint8_t>(nf);
d.f_h = CudaDevicePtr<int32_t>(nf, ALLOC); d.f_k = CudaDevicePtr<int32_t>(nf, ALLOC); d.f_l = CudaDevicePtr<int32_t>(nf, ALLOC);
d.f_frame = CudaDevicePtr<int32_t>(nf, ALLOC); d.f_group = CudaDevicePtr<int32_t>(nf, ALLOC);
d.f_I = CudaDevicePtr<float>(nf, ALLOC); d.f_sigma = CudaDevicePtr<float>(nf, ALLOC);
d.f_d = CudaDevicePtr<float>(nf, ALLOC); d.f_img = CudaDevicePtr<float>(nf, ALLOC);
d.f_px = CudaDevicePtr<float>(nf, ALLOC); d.f_py = CudaDevicePtr<float>(nf, ALLOC);
d.f_var_bkg = CudaDevicePtr<float>(nf, ALLOC); d.f_var_per_I = CudaDevicePtr<float>(nf, ALLOC);
d.f_on_ice = CudaDevicePtr<uint8_t>(nf, ALLOC); d.f_clipped = CudaDevicePtr<uint8_t>(nf, ALLOC);
d.f_corr = CudaDevicePtr<float>(nf, ALLOC); d.f_partiality = CudaDevicePtr<float>(nf, ALLOC);
d.f_rlp = CudaDevicePtr<float>(nf, ALLOC); d.f_zeta = CudaDevicePtr<float>(nf, ALLOC);
d.f_inv_sigma = CudaDevicePtr<double>(nf, ALLOC);
d.f_sco_coeff = CudaDevicePtr<float>(nf, ALLOC); d.f_sco_ok = CudaDevicePtr<uint8_t>(nf, ALLOC);
CudaCheck(cudaMemcpy(d.rr_offset.get(), offset.data(), size_t(d.n_runs) * sizeof(int32_t),
cudaMemcpyHostToDevice), "upload offset");
rugnux/HotPixelsGPU.cu
+11 -3
View File
@@ -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<uint32_t> &host_key_begin, int nrings, int sectors)
: npixels(npixels), nkeys(static_cast<size_t>(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<uint32_t> pixels(host_key_begin.back());
std::vector<uint32_t> filled(host_key_begin.begin(), host_key_begin.end() - 1);
for (size_t i = 0; i < npixels; i++)