diff --git a/image_analysis/scale_merge/ErrorModelGPU.cu b/image_analysis/scale_merge/ErrorModelGPU.cu index ef1302b34..8d56d5631 100644 --- a/image_analysis/scale_merge/ErrorModelGPU.cu +++ b/image_analysis/scale_merge/ErrorModelGPU.cu @@ -89,8 +89,8 @@ ErrorModelFit FitErrorModelGPU(const std::vector &pool) { return fit; const int n = static_cast(pool.size()); - // The rank sort is the peak: the pool, two key and two index arrays and the sort scratch. - constexpr size_t BYTES_PER_SAMPLE = 32 + 2 * 8 + 2 * 4 + 16; + // The rank sort is the peak: the pool, two key and two index arrays. + constexpr size_t BYTES_PER_SAMPLE = 32 + 2 * 8 + 2 * 4; size_t free_bytes = 0; if (!cuda_memory_available(BYTES_PER_SAMPLE * pool.size(), free_bytes)) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, @@ -104,22 +104,25 @@ ErrorModelFit FitErrorModelGPU(const std::vector &pool) { stream)); // The samples in the host's rank order - (I2/s2, s2, I2, dev2), a total order - by four stable - // passes, least significant key first. The equal-count bins are then consecutive stretches. + // passes, least significant key first. The equal-count bins are then consecutive stretches. Each + // pass goes back and forth between the two key and the two index arrays, so that cub holds no copy + // of them of its own; the same stable sort, so the same order. CudaDevicePtr idx(n, ALLOC), idx_next(n, ALLOC); CudaDevicePtr key(n, ALLOC), key_sorted(n, ALLOC); IotaKernel<<>>(n, idx); cuda_err(cudaGetLastError()); + cub::DoubleBuffer keys(key.get(), key_sorted.get()); + cub::DoubleBuffer order(idx.get(), idx_next.get()); size_t bytes = 0; - cuda_err(cub::DeviceRadixSort::SortPairs(nullptr, bytes, key.get(), key_sorted.get(), idx.get(), idx_next.get(), - n, 0, 64, stream)); + cuda_err(cub::DeviceRadixSort::SortPairs(nullptr, bytes, keys, order, n, 0, 64, stream)); CudaDevicePtr scratch(std::max(bytes, 1), ALLOC); for (int field = 0; field < 4; ++field) { - KeyKernel<<>>(n, d_pool, idx, field, key); + KeyKernel<<>>(n, d_pool, order.Current(), field, keys.Current()); cuda_err(cudaGetLastError()); - cuda_err(cub::DeviceRadixSort::SortPairs(scratch.get(), bytes, key.get(), key_sorted.get(), idx.get(), - idx_next.get(), n, 0, 64, stream)); - std::swap(idx, idx_next); + cuda_err(cub::DeviceRadixSort::SortPairs(scratch.get(), bytes, keys, order, n, 0, 64, stream)); } + if (keys.Current() != key_sorted.get()) std::swap(key, key_sorted); + if (order.Current() != idx.get()) std::swap(idx, idx_next); scratch = CudaDevicePtr(); idx_next = CudaDevicePtr(); key = CudaDevicePtr(); // After the last pass key_sorted holds I2/s2 in rank order. CudaDevicePtr s2(n, ALLOC), I2(n, ALLOC), dev2(n, ALLOC); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 7595dcca8..7ab124f79 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -608,18 +608,23 @@ void RotationScaleMerge::IngestHost() { } #ifdef JFJOCH_USE_CUDA -// The resident path's ingest. The reflections go to the card once, frame by frame, and the card measures -// the per-frame background, sorts, cuts the raw-hkl runs and drops what the resolution range excludes -// (RotationScaleMergeGPU::IngestUpload .. IngestSort) - the observations, the runs and the order of each -// that BuildInRangeObservations lays down on the host, to the bit, and left where the device pipeline -// reads them. The geometry smoothing then runs here on what the card hands back. +// The resident path's ingest. The reflections go to the card frame by frame, and the card measures the +// per-frame background, sorts, cuts the raw-hkl runs and drops what the resolution range excludes +// (RotationScaleMergeGPU::IngestUpload .. IngestSort); once the scaling is known to fit, they go to the +// card a second time for the observations the cut kept (IngestObservations) - the observations, the runs +// and the order of each that BuildInRangeObservations lays down on the host, to the bit, and left where +// the device pipeline reads them. The geometry smoothing then runs here on what the card hands back. void RotationScaleMerge::IngestOnDevice() { std::vector refl(n_frames); std::vector count(n_frames); + size_t n_reflections = 0; for (int o = 0; o < n_frames; ++o) { refl[o] = partials_out[o].reflections.data(); count[o] = static_cast(partials_out[o].reflections.size()); + n_reflections += count[o]; } + if (ScaleOnCPUWhereSortTooLarge(n_reflections)) + return; TimingMark upload_mark(*logger, "rsm ingest upload"); gpu_->IngestUpload(n_frames, refl.data(), count.data(), nthreads); upload_mark.Stop(); @@ -627,7 +632,7 @@ void RotationScaleMerge::IngestOnDevice() { std::vector mean_bkg(n_frames); gpu_->IngestMeanBackground(mean_bkg.data()); MeasureIncidentFlux(mean_bkg); - const int n = gpu_->IngestSort(frame_flux.data(), d_min_limit, d_max_limit); + const int n = gpu_->IngestSort(d_min_limit, d_max_limit); n_partials_obs = n; RotationScaleMergeGPU::IngestLayout layout; gpu_->IngestGetLayout(layout); @@ -644,9 +649,13 @@ void RotationScaleMerge::IngestOnDevice() { rawrun_l = std::move(layout.run_l); rawrun_d = std::move(layout.run_d); // The scaling's footprint is known from here on: where the card cannot take it, the card's work - // ends before the smoothing. + // ends before a single observation is on it. if (ScaleOnCPUWhereTooLarge()) return; + auto pos = std::make_shared(); + pos->px.resize(n); pos->py.resize(n); + gpu_->IngestObservations(refl.data(), count.data(), nthreads, frame_flux.data(), pos->px.data(), pos->py.data()); + obs_position = std::move(pos); const SmoothedGeometry sg = SmoothedFrameGeometry(); int64_t n_resolved = 0; @@ -703,13 +712,28 @@ void RotationScaleMerge::IngestOnDevice() { // scales on the CPU - decided here, once, on the card's total memory, so the answer depends on the card // but not on the moment. bool RotationScaleMerge::ScaleOnCPUWhereTooLarge() { + return ScaleOnCPUWhereOverCard(GPUScalingPeakBytes(), fmt::format("scaling {} observations of {} distinct hkl", + n_partials_obs, rawrun_start.size())); +} + +// Before the sort no observation is on the card yet, nor is it known how many there are: what the sort +// and the cut hold (IngestUpload .. IngestSort), measured at 32 bytes per reflection, or 20 and 56 per +// raw-hkl run where the runs are many - 48 per reflection with half as many runs as reflections, the +// most measured being 0.47 - and the two staging buffers. A set that the card cannot even sort scales +// on the CPU without touching the card. +bool RotationScaleMerge::ScaleOnCPUWhereSortTooLarge(size_t n_reflections) { + constexpr size_t SORT_BYTES_PER_REFLECTION = 48, STAGING_BYTES = size_t(128) << 20; + return ScaleOnCPUWhereOverCard(SORT_BYTES_PER_REFLECTION * n_reflections + STAGING_BYTES, + fmt::format("sorting {} reflections", n_reflections)); +} + +bool RotationScaleMerge::ScaleOnCPUWhereOverCard(size_t bytes, const std::string &what) { const size_t card = get_gpu_total_memory(gpu_device); - if (GPUScalingPeakBytes() <= CPU_SCALING_ABOVE_CARD_FRACTION * card) + if (bytes <= CPU_SCALING_ABOVE_CARD_FRACTION * card) return false; - logger->Warning("RotationScaleMerge: scaling {} observations of {} distinct hkl on the GPU needs about " - "{:.1f} GB, over {:.0f}% of the card's {:.1f} GB - scaling on the CPU instead", - n_partials_obs, rawrun_start.size(), GPUScalingPeakBytes() / 1e9, - 100.0 * CPU_SCALING_ABOVE_CARD_FRACTION, card / 1e9); + logger->Warning("RotationScaleMerge: {} on the GPU needs about {:.1f} GB, over {:.0f}% of the card's {:.1f} GB " + "- scaling on the CPU instead", what, bytes / 1e9, 100.0 * CPU_SCALING_ABOVE_CARD_FRACTION, + card / 1e9); scale_on_cpu_ = true; gpu_active_ = false; return true; @@ -725,12 +749,13 @@ bool RotationScaleMerge::OnGPU() const { } size_t RotationScaleMerge::GPUScalingPeakBytes() const { - // The device scaling's peak from this decision on, fitted to the engine's device buffers measured at - // their peak on fourteen engines of five large data sets (23-123 M observations, 1.6-15 M distinct - // hkl): within 12% on each, the furthest under a P1 merge of a large cell. The ingest's own peak, 82 - // to 100 bytes per reflection, has passed by the time this decides. Unlike DeviceBytes, an estimate - // and not an upper bound - it decides where the GPU clearly cannot take the scaling. - constexpr size_t PEAK_BYTES_PER_OBS = 63, PEAK_BYTES_PER_HKL = 115; + // The device scaling's peak from this decision on - the observations taken from the reflections + // (IngestObservations) included, the sort before it not - fitted to the engine's device buffers + // measured at their peak on fourteen engines of five large data sets (23-123 M observations, + // 1.6-15 M distinct hkl): within 9% on each, the furthest under a P1 merge of a large cell. Unlike + // DeviceBytes, an estimate and not an upper bound - it decides where the GPU clearly cannot take the + // scaling. + constexpr size_t PEAK_BYTES_PER_OBS = 56, PEAK_BYTES_PER_HKL = 99; return PEAK_BYTES_PER_OBS * static_cast(n_partials_obs) + PEAK_BYTES_PER_HKL * rawrun_start.size(); } @@ -738,22 +763,21 @@ size_t RotationScaleMerge::DeviceBytes() const { // What each costs at the peak: an observation's arrays, and its group CSR until the combine; a full's // arrays with its error-model sample, or with the surface terms; an hkl's raw run and its ASU group's // sums. Fitted as an upper bound to the engine's device buffers measured at their peak on fourteen - // engines of five large data sets (23-123 M observations): up to 18% over. The events are counted + // engines of five large data sets (23-123 M observations): up to 14% over. The events are counted // on the device path's ingest; elsewhere every observation may be an event of its own. - constexpr size_t BYTES_PER_OBS = 66, BYTES_PER_EVENT = 94, BYTES_PER_HKL = 52; + constexpr size_t BYTES_PER_OBS = 58, BYTES_PER_EVENT = 30, BYTES_PER_HKL = 86; const auto events = static_cast(resident_ingest ? rocking_events_at_ingest : n_partials_obs); const size_t bytes = BYTES_PER_OBS * static_cast(n_partials_obs) + BYTES_PER_EVENT * events + BYTES_PER_HKL * rawrun_start.size(); - // The resident ingest's own peak comes before any of that, over every reflection read and every run - // before the resolution cut, measured on the same engines. Beside the observation arrays and the - // integrated delta_phi (58 bytes), each radix sort holds two pairs of keys and values, 82 bytes per - // reflection in all; the cut (IngestSort) of the runs 74 per reflection and 56 per run, that of the - // observations after it 89 per reflection and 24 per run kept. The upload before them holds the - // observation arrays, the sort key and delta_phi (66 bytes) and the two 64 MB staging buffers. - constexpr size_t INGEST_BYTES_PER_REFLECTION = 89, INGEST_BYTES_PER_RUN = 56; - constexpr size_t INGEST_ARRAY_BYTES = 66, INGEST_STAGING_BYTES = size_t(160) << 20; - const size_t ingest = std::max(INGEST_BYTES_PER_REFLECTION * ingest_reflections + INGEST_BYTES_PER_RUN * ingest_runs, - INGEST_ARRAY_BYTES * ingest_reflections + INGEST_STAGING_BYTES); + // The resident ingest's own peak, over every reflection read and every run before the resolution + // cut, measured on the same engines: the observations taken from the reflections hold the + // observation arrays, the place the cut gave each reflection and the runs it kept, 60 to 72 bytes + // per reflection; the sort and the cut before them 32 per reflection and 56 per run. The two 64 MB + // staging buffers beside either. + constexpr size_t INGEST_BYTES_PER_REFLECTION = 60, INGEST_BYTES_PER_RUN = 42; + constexpr size_t INGEST_STAGING_BYTES = size_t(160) << 20; + const size_t ingest = INGEST_BYTES_PER_REFLECTION * ingest_reflections + INGEST_BYTES_PER_RUN * ingest_runs + + INGEST_STAGING_BYTES; return std::max(bytes + bytes / 5, ingest); } @@ -799,7 +823,7 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro using F = RotationScaleMergeGPU::ObsField; constexpr int STAGE_SLICE = 8 << 20; const int slice = std::min(std::max(n, 1), STAGE_SLICE); - constexpr int NFLOAT = 12; + constexpr int NFLOAT = 10; std::unique_ptr fbuf[NFLOAT]; for (auto &b : fbuf) b.reset(new float[slice]); std::unique_ptr frm(new int32_t[slice]); @@ -812,19 +836,25 @@ void RotationScaleMerge::IngestDevice(const RotationScaleMerge *observations_fro fbuf[0][i] = o.I; fbuf[1][i] = o.sigma; fbuf[2][i] = o.prescaling_corr; fbuf[3][i] = o.partiality; fbuf[4][i] = o.zeta; fbuf[5][i] = o.corr; fbuf[6][i] = o.bkg; fbuf[7][i] = o.var_bkg; fbuf[8][i] = o.image_number; - fbuf[9][i] = o.d; fbuf[10][i] = o.px; fbuf[11][i] = o.py; + fbuf[9][i] = o.d; frm[i] = o.frame; onice[i] = o.on_ice; clip[i] = o.clipped; } }); constexpr F fields[NFLOAT] = {F::I, F::Sigma, F::PrescalingCorr, F::Partiality, F::Zeta, F::Corr0, F::Bkg, F::VarBkg, F::ImageNumber, - F::D, F::Px, F::Py}; + F::D}; for (int k = 0; k < NFLOAT; ++k) gpu_->SetObsField(fields[k], at, cnt, fbuf[k].get()); gpu_->SetObsFrame(at, cnt, frm.get()); gpu_->SetObsOnIce(at, cnt, onice.get()); gpu_->SetObsClipped(at, cnt, clip.get()); } + auto pos = std::make_shared(); + pos->px.resize(n); pos->py.resize(n); + ParallelChunks(n, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) { pos->px[i] = partials[i].px; pos->py[i] = partials[i].py; } + }); + obs_position = std::move(pos); gpu_->SetRawRuns(static_cast(rawrun_start.size()), static_cast(perm.size()), perm.data(), rawrun_start.data(), rawrun_count.data(), rawrun_h.data(), rawrun_k.data(), rawrun_l.data()); @@ -6706,23 +6736,24 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool fulls.resize(nf); FullsStaging &st = fulls_staging; st.Resize(nf); - gpu_->GetFulls(st.h.data(), st.k.data(), st.l.data(), st.I.data(), st.sigma.data(), st.d.data(), + gpu_->GetFulls(st.run.data(), st.I.data(), st.sigma.data(), st.d.data(), st.image_number.data(), st.frame.data(), st.on_ice.data(), st.clipped.data(), st.group.data()); - gpu_->GetFullsPxPy(st.px.data(), st.py.data()); + gpu_->GetFullsPeak(st.peak.data()); + const ObsPosition &pos = *obs_position; gpu_->GetFullsVariance(st.var_bkg.data(), st.var_per_I.data(), st.capture.data()); if (scaled_fulls_on_gpu) gpu_->GetFullsCorr(st.corr.data()); gpu_->FullsDownloaded(); ParallelChunks(nf, ThreadsForWork(nf, nthreads), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) { Obs &o = fulls[i]; - o.h = st.h[i]; o.k = st.k[i]; o.l = st.l[i]; + o.h = rawrun_h[st.run[i]]; o.k = rawrun_k[st.run[i]]; o.l = rawrun_l[st.run[i]]; o.I = st.I[i]; o.sigma = st.sigma[i]; o.d = st.d[i]; o.var_bkg = st.var_bkg[i]; o.var_per_I = st.var_per_I[i]; o.capture = st.capture[i]; o.prescaling_corr = 1.0f; o.partiality = 1.0f; o.corr = scaled_fulls_on_gpu ? st.corr[i] : 1.0f; o.image_number = st.image_number[i]; o.frame = st.frame[i]; - o.px = st.px[i]; o.py = st.py[i]; + o.px = pos.px[st.peak[i]]; o.py = pos.py[st.peak[i]]; o.on_ice = st.on_ice[i]; o.clipped = st.clipped[i]; o.group = st.group[i]; // A full has no rocking geometry of its own and carries no background; the combine on // the device does not produce these, and a reused buffer would keep the last chain's. diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 61c13378c..514378909 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -173,8 +173,8 @@ public: // card. Call after IngestHost(). It counts what the ingest knows - the observations, the rocking // events (the combine makes at most one full of each) and the distinct hkl (at most one ASU group // each) - at what each costs at the engine's peak (the correction-surface terms, or the merge with - // its outlier test), plus a fifth. Measured on 11 data sets of 13-123 M observations, every - // engine's peak was at most 81% of this. Where the ingest on the card peaks higher, its peak. + // its outlier test), plus a fifth. Measured on 14 engines of 23-123 M observations, every engine's + // peak was at most 88% of this. Where the ingest on the card peaks higher, its peak. // Call after IngestDevice(). [[nodiscard]] size_t DeviceBytes() const; // The device scaling's estimated peak (see the .cpp). Call after the ingest. @@ -370,6 +370,10 @@ private: bool partials_released = false; int rocking_event_frames_at_ingest = 0; int64_t rocking_events_at_ingest = 0; + // Each partial's predicted detector position, for the fulls the device path combines (GetFullsPeak): + // only the host reads it, so the card does not hold it. Shared by the engines CloneIngest makes. + struct ObsPosition { std::vector px, py; }; + std::shared_ptr obs_position; // The reflections and raw-hkl runs the resident ingest read before the resolution cut (DeviceBytes). size_t ingest_reflections = 0, ingest_runs = 0; std::vector frame_start, frame_count; // CSR ranges of `partials` per frame @@ -446,13 +450,13 @@ private: // scale->combine->merge chain runs several times per run and these are a few hundred megabytes // between them, so as locals every chain allocates, faults in and zeroes the lot again. struct FullsStaging { - std::vector h, k, l, frame, group; - std::vector I, sigma, d, image_number, corr, px, py, var_bkg, var_per_I, capture; + std::vector run, peak, frame, group; + std::vector I, sigma, d, image_number, corr, var_bkg, var_per_I, capture; std::vector on_ice, clipped; void Resize(int n) { - h.resize(n); k.resize(n); l.resize(n); frame.resize(n); group.resize(n); + run.resize(n); peak.resize(n); frame.resize(n); group.resize(n); I.resize(n); sigma.resize(n); d.resize(n); image_number.resize(n); corr.resize(n); - px.resize(n); py.resize(n); var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); + var_bkg.resize(n); var_per_I.resize(n); capture.resize(n); on_ice.resize(n); clipped.resize(n); } }; @@ -605,6 +609,10 @@ private: // Where the scaling's estimated peak is over CPU_SCALING_ABOVE_CARD_FRACTION of the card: says so, // gives the device up (scale_on_cpu_) and returns true. Call once the ingest has counted. bool ScaleOnCPUWhereTooLarge(); + // The same for the ingest's sort of n_reflections, decided before anything is on the card. + bool ScaleOnCPUWhereSortTooLarge(size_t n_reflections); + // What the rule above says for `bytes` of the card, `what` naming the work in the warning. + bool ScaleOnCPUWhereOverCard(size_t bytes, const std::string &what); #endif // Compute the dense ASU-group id for the current space group by grouping the (pre-sorted) raw-hkl diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.cu b/image_analysis/scale_merge/RotationScaleMergeGPU.cu index e1c425eba..612072e94 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.cu +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.cu @@ -294,7 +294,7 @@ namespace { float max_frame_gap; const ObsRecord *__restrict__ obs; const float *__restrict__ corr, *__restrict__ bkg, *__restrict__ var_bkg, - *__restrict__ image_number, *__restrict__ d, *__restrict__ px, *__restrict__ py; + *__restrict__ image_number, *__restrict__ d; const int32_t *__restrict__ frame; const uint8_t *__restrict__ on_ice, *__restrict__ clipped; const int32_t *__restrict__ perm, *__restrict__ rr_start, *__restrict__ rr_count, @@ -303,8 +303,8 @@ namespace { int32_t *rr_nevents; // count pass output int32_t *rr_noverloaded; // count pass output: events dropped for a saturated pixel const int32_t *rr_offset; // emit pass: per-run base offset into the fulls arrays - int32_t *f_h, *f_k, *f_l, *f_frame, *f_group; - float *f_I, *f_sigma, *f_d, *f_img, *f_px, *f_py, *f_var_bkg, *f_var_per_I, *f_capture; + int32_t *f_run, *f_peak, *f_frame, *f_group; + float *f_I, *f_sigma, *f_d, *f_img, *f_var_bkg, *f_var_per_I, *f_capture; uint8_t *f_on_ice, *f_clipped; }; @@ -363,7 +363,7 @@ namespace { const int first = p.perm[ev_start]; int peak_outcome = p.frame[first]; float peak_frame = p.image_number[first]; - float peak_px = p.px[first], peak_py = p.py[first]; + int peak_obs = first; float peak_partiality = -1.0f; const bool on_ice = p.on_ice[first]; uint8_t flags = 0; @@ -383,7 +383,7 @@ namespace { peak_partiality = p.obs[i].partiality; peak_outcome = p.frame[i]; peak_frame = p.image_number[i]; - peak_px = p.px[i]; peak_py = p.py[i]; + peak_obs = i; } if (!isfinite(d) && isfinite(p.d[i]) && p.d[i] > 0.0f) d = p.d[i]; } @@ -432,7 +432,7 @@ namespace { if (Emit) { const int o = p.rr_offset[r] + n_emit; - p.f_h[o] = p.rr_h[r]; p.f_k[o] = p.rr_k[r]; p.f_l[o] = p.rr_l[r]; + p.f_run[o] = r; p.f_I[o] = float(F); p.f_sigma[o] = float(sigma_full); p.f_var_bkg[o] = float(var_bkg_full); @@ -440,7 +440,7 @@ namespace { p.f_capture[o] = float(capture); p.f_d[o] = d; p.f_img[o] = peak_frame; - p.f_px[o] = peak_px; p.f_py[o] = peak_py; + p.f_peak[o] = peak_obs; p.f_frame[o] = peak_outcome; p.f_on_ice[o] = on_ice ? 1 : 0; p.f_clipped[o] = clipped ? 1 : 0; @@ -672,7 +672,7 @@ namespace { double min_partiality, error_model_a, error_model_b, reject_nsigma; double reject_up, reject_down; // the multiplicative part of the outlier band (OutlierBand.h) int for_search, error_model_active, reject_outliers; - const float *I, *sigma, *corr, *partiality, *d, *reject_median, *var_bkg, *var_per_I, *capture; + const float *I, *sigma, *corr, *d, *reject_median, *var_bkg, *var_per_I, *capture; const float *reject_var_add; // per group: the shell's measured share of dI^2/4 const int32_t *group, *frame; const uint8_t *on_ice, *frame_cell_ok, *half; @@ -693,6 +693,9 @@ namespace { uint8_t *rejected_obs; // per-full: set by MergeAccum (outlier-rejected), read by MergeRmeas }; + // A full's partiality: the device merges the fulls it scaled, in the Unity model. + constexpr float FULL_PARTIALITY = 1.0f; + // A full passes the merge / error-model filter (mirrors MergeAndStats::usable_merge). __device__ __forceinline__ bool MergeUsable(int i, const MergeParams &p) { const int g = p.group[i]; @@ -701,7 +704,7 @@ namespace { const float c = p.corr[i]; if (!(c > 0.0f) || !isfinite(c)) return false; if (p.for_search && p.on_ice[i]) return false; - if (p.partiality[i] < p.min_partiality) return false; + if (FULL_PARTIALITY < p.min_partiality) return false; const float I_corr = p.I[i] * c, sigma_corr = p.sigma[i] * c; return isfinite(I_corr) && isfinite(sigma_corr) && sigma_corr > 0.0f; } @@ -713,7 +716,7 @@ namespace { if (!p.frame_cell_ok[p.frame[i]]) return false; const float c = p.corr[i]; if (!(c > 0.0f) || !isfinite(c)) return false; - if (p.partiality[i] < p.min_partiality) return false; + if (FULL_PARTIALITY < p.min_partiality) return false; const float I_corr = p.I[i] * c, sigma_corr = p.sigma[i] * c; return isfinite(I_corr) && isfinite(sigma_corr) && sigma_corr > 0.0f; } @@ -838,10 +841,12 @@ namespace { } // One thread per full: its observation key, and whether the merge can use it. - __global__ void ObservationKeyKernel(int n_obs, MergeParams p, const int32_t *h, const int32_t *k, - const int32_t *l, uint64_t *key, uint8_t *usable) { + // A full's raw hkl is its run's. + __global__ void ObservationKeyKernel(int n_obs, MergeParams p, const int32_t *run, const int32_t *h, + const int32_t *k, const int32_t *l, uint64_t *key, uint8_t *usable) { for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n_obs; i += gridDim.x * blockDim.x) { - key[i] = ObservationKeyDevice(h[i], k[i], l[i], p.frame[i]); + const int r = run[i]; + key[i] = ObservationKeyDevice(h[r], k[r], l[r], p.frame[i]); usable[i] = MergeUsable(i, p) ? 1 : 0; } } @@ -1257,19 +1262,20 @@ struct RotationScaleMergeGPU::Impl { int merge_em_active = 0; // combine: extra per-obs inputs + the one-time raw-hkl run layout - CudaDevicePtr bkg, var_bkg, image_number, d_obs, px_obs, py_obs; + CudaDevicePtr bkg, var_bkg, image_number, d_obs; int n_runs = 0, n_perm = 0; CudaDevicePtr perm, rr_start, rr_count, rr_h, rr_k, rr_l; // combine: resident fulls SoA (rebuilt each Combine) int n_fulls = 0; - CudaDevicePtr f_h, f_k, f_l, f_frame, f_group; - CudaDevicePtr f_I, f_sigma, f_d, f_img, f_px, f_py, f_var_bkg, f_var_per_I, f_capture; + // run: the raw-hkl run a full is of; peak: the partial its position is taken from + CudaDevicePtr f_run, f_peak, f_frame, f_group; + CudaDevicePtr f_I, f_sigma, f_d, f_img, f_var_bkg, f_var_per_I, f_capture; CudaDevicePtr f_on_ice, f_clipped; // scale-fulls (Unity model, kept resident): all-ones partiality/zeta and records with partiality and // prescaling_corr 1, so the shared scaling kernels yield coeff=mean, plus the working corr and the // fulls frame/group CSRs (built on the host from the small f_frame/f_group key arrays, over the // emit-ordered fulls). - CudaDevicePtr f_corr, f_partiality, f_zeta; + CudaDevicePtr f_corr, f_zeta; CudaDevicePtr f_obs; CudaDevicePtr f_frame_perm, f_frame_start, f_frame_count; CudaDevicePtr f_gperm, f_gstart, f_gcount; @@ -1285,12 +1291,22 @@ struct RotationScaleMergeGPU::Impl { CudaDevicePtr s_perm[3], s_seg_start[3]; CudaDevicePtr s_A, s_sw, s_swI, s_tcross, s_tref2, s_cross, s_ref2; CudaDevicePtr s_w_Is, s_w_Iref, s_Iref; // per term of a subset, in segment order - // the resident ingest (IngestUpload .. IngestFinish): each frame's first reflection - // (host, length n_frames + 1), per observation its raw hkl as the sort key until the sort has used - // it, and its delta_phi as integrated until the smoothing has + // the resident ingest (IngestUpload .. IngestFinish): each frame's first reflection (host, length + // n_frames + 1); per reflection what the sort and the cut read - its raw hkl as the sort key, its + // frame position, resolution and background - until they have, and the place the cut gives it (as + // an exclusive prefix over the kept ones, n_keep in all) until the observations are taken from the + // reflections; per observation its delta_phi as integrated until the smoothing has used it std::vector in_src_start; CudaDevicePtr in_key; + CudaDevicePtr in_img, in_d, in_bkg; + CudaDevicePtr in_new_idx; + int in_n_keep = 0; CudaDevicePtr in_dphi; + // The reflections of the ingest's frames, through two pinned buffers of whole frames, each filled on + // all threads while the other is on its way to the card; take(cnt, base, src) takes cnt of them, + // the first being reflection base, apart on the card. + template + void StageReflections(const Reflection *const *refl, const int32_t *count, size_t nthreads, Take take); RotationScaleMergeGPU::IngestLayout in_layout; }; @@ -1362,7 +1378,6 @@ void RotationScaleMergeGPU::Impl::AllocObservations(int n_in) { zeta = Alloc(n); corr = Alloc(n); bkg = Alloc(n); var_bkg = Alloc(n); image_number = Alloc(n); d_obs = Alloc(n); - px_obs = Alloc(n); py_obs = Alloc(n); frame = Alloc(n); on_ice = Alloc(n); clipped = Alloc(n); @@ -1386,10 +1401,10 @@ void RotationScaleMergeGPU::Impl::AllocFrames(const int32_t *fstart, const int32 void RotationScaleMergeGPU::Impl::ReleaseFulls() { n_fulls = 0; - f_h = {}; f_k = {}; f_l = {}; f_frame = {}; f_group = {}; - f_I = {}; f_sigma = {}; f_d = {}; f_img = {}; f_px = {}; f_py = {}; + f_run = {}; f_peak = {}; f_frame = {}; f_group = {}; + f_I = {}; f_sigma = {}; f_d = {}; f_img = {}; f_var_bkg = {}; f_var_per_I = {}; f_capture = {}; f_on_ice = {}; f_clipped = {}; - f_corr = {}; f_partiality = {}; f_zeta = {}; f_obs = {}; + f_corr = {}; f_zeta = {}; f_obs = {}; f_frame_perm = {}; f_frame_start = {}; f_frame_count = {}; f_gperm = {}; f_gstart = {}; f_gcount = {}; m_em_mean = {}; m_hand = {}; m_has_hands = {}; m_rejected = {}; @@ -1434,8 +1449,6 @@ void RotationScaleMergeGPU::SetObsField(ObsField f, int offset, int count, const case ObsField::VarBkg: UploadChunk(d.var_bkg, offset, count, v, impl_->s()); break; case ObsField::ImageNumber: UploadChunk(d.image_number, offset, count, v, impl_->s()); break; case ObsField::D: UploadChunk(d.d_obs, offset, count, v, impl_->s()); break; - case ObsField::Px: UploadChunk(d.px_obs, offset, count, v, impl_->s()); break; - case ObsField::Py: UploadChunk(d.py_obs, offset, count, v, impl_->s()); break; } } @@ -1457,7 +1470,7 @@ void RotationScaleMergeGPU::CopyObservationsFrom(const RotationScaleMergeGPU &fr const size_t n = d.n_obs; copy(d.obs, f.obs, n); copy(d.zeta, f.zeta, n); copy(d.corr, f.corr, n); copy(d.bkg, f.bkg, n); copy(d.var_bkg, f.var_bkg, n); copy(d.image_number, f.image_number, n); - copy(d.d_obs, f.d_obs, n); copy(d.px_obs, f.px_obs, n); copy(d.py_obs, f.py_obs, n); + copy(d.d_obs, f.d_obs, n); copy(d.frame, f.frame, n); copy(d.on_ice, f.on_ice, n); copy(d.clipped, f.clipped, n); const size_t nr = std::max(1, d.n_runs); d.perm = d.Alloc(std::max(1, d.n_perm)); @@ -1619,7 +1632,7 @@ void RotationScaleMergeGPU::MergeEmSamples(bool for_search, double min_partialit MergeParams p{}; p.n_groups = ng; p.min_partiality = min_partiality; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1702,7 +1715,7 @@ void RotationScaleMergeGPU::MergeAccum(double error_model_a, double error_model_ p.reject_down = -std::expm1(-reject_nsigma * reject_mult_sd); p.reject_median = median.get(); p.reject_var_add = var_add.get(); - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.half = d_half.get(); @@ -1760,7 +1773,7 @@ void RotationScaleMergeGPU::MergeRmeas(const double *merged_I, double *absdev, d MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); p.merged_I = d_merged_I.get(); @@ -1797,7 +1810,7 @@ void RotationScaleMergeGPU::MergeBijvoetExcess(double error_model_a, double erro MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1829,7 +1842,7 @@ void RotationScaleMergeGPU::MergeGroupScatter(double *sum, double *sum_sq, doubl MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); @@ -1856,15 +1869,15 @@ void RotationScaleMergeGPU::MergeHalves(uint64_t *obs_key, uint8_t *half) { MergeParams p{}; p.n_groups = ng; p.min_partiality = d.merge_min_part; p.for_search = d.merge_for_search; - p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.partiality = d.f_partiality.get(); + p.I = d.f_I.get(); p.sigma = d.f_sigma.get(); p.corr = d.f_corr.get(); p.d = d.f_d.get(); p.group = d.f_group.get(); p.frame = d.f_frame.get(); p.on_ice = d.f_on_ice.get(); p.frame_cell_ok = d.frame_cell_ok.get(); p.gperm = d.f_gperm.get(); p.gstart = d.f_gstart.get(); p.gcount = d.f_gcount.get(); const int obs_blocks = std::min(65535, (nf + BLK - 1) / BLK); const int grp_blocks = std::min(65535, (ng + BLK - 1) / BLK); - ObservationKeyKernel<<s()>>>(nf, p, d.f_h.get(), d.f_k.get(), d.f_l.get(), key.get(), - usable.get()); + ObservationKeyKernel<<s()>>>(nf, p, d.f_run.get(), d.rr_h.get(), d.rr_k.get(), + d.rr_l.get(), key.get(), usable.get()); CudaCheck(cudaGetLastError(), "observation key launch"); CudaCheck(cudaMemsetAsync(h.get(), 0, size_t(nf), impl_->s()), "memset halves"); HalvesKernel<<s()>>>(p, key.get(), usable.get(), h.get()); @@ -1981,7 +1994,6 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti p.max_frame_gap = max_frame_gap; p.obs = d.obs.get(); p.corr = d.corr.get(); p.bkg = d.bkg.get(); p.var_bkg = d.var_bkg.get(); p.image_number = d.image_number.get(); p.d = d.d_obs.get(); - p.px = d.px_obs.get(); p.py = d.py_obs.get(); p.frame = d.frame.get(); p.on_ice = d.on_ice.get(); p.clipped = d.clipped.get(); p.perm = d.perm.get(); p.rr_start = d.rr_start.get(); p.rr_count = d.rr_count.get(); p.rr_h = d.rr_h.get(); p.rr_k = d.rr_k.get(); p.rr_l = d.rr_l.get(); p.rr_group = rr_group.get(); @@ -2011,23 +2023,21 @@ 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 = d.Alloc(nf); d.f_k = d.Alloc(nf); d.f_l = d.Alloc(nf); + d.f_run = d.Alloc(nf); d.f_peak = d.Alloc(nf); d.f_frame = d.Alloc(nf); d.f_group = d.Alloc(nf); d.f_I = d.Alloc(nf); d.f_sigma = d.Alloc(nf); d.f_d = d.Alloc(nf); d.f_img = d.Alloc(nf); - d.f_px = d.Alloc(nf); d.f_py = d.Alloc(nf); d.f_var_bkg = d.Alloc(nf); d.f_var_per_I = d.Alloc(nf); d.f_capture = d.Alloc(nf); d.f_on_ice = d.Alloc(nf); d.f_clipped = d.Alloc(nf); - d.f_corr = d.Alloc(nf); d.f_partiality = d.Alloc(nf); + d.f_corr = d.Alloc(nf); d.f_zeta = d.Alloc(nf); d.f_obs = d.Alloc(nf); d.Upload(rr_offset, offset.data(), d.n_runs); p.rr_offset = rr_offset.get(); - p.f_h = d.f_h.get(); p.f_k = d.f_k.get(); p.f_l = d.f_l.get(); + p.f_run = d.f_run.get(); p.f_peak = d.f_peak.get(); p.f_frame = d.f_frame.get(); p.f_group = d.f_group.get(); p.f_I = d.f_I.get(); p.f_sigma = d.f_sigma.get(); p.f_d = d.f_d.get(); p.f_img = d.f_img.get(); - p.f_px = d.f_px.get(); p.f_py = d.f_py.get(); p.f_var_bkg = d.f_var_bkg.get(); p.f_var_per_I = d.f_var_per_I.get(); p.f_capture = d.f_capture.get(); p.f_on_ice = d.f_on_ice.get(); p.f_clipped = d.f_clipped.get(); if (d.n_fulls > 0) { @@ -2038,7 +2048,7 @@ int RotationScaleMergeGPU::Combine(const int32_t *rawrun_group, double min_parti return d.n_fulls; } -void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, +void RotationScaleMergeGPU::GetFulls(int32_t *run, float *I, float *sigma, float *d, float *image_number, int32_t *frame, uint8_t *on_ice, uint8_t *clipped, int32_t *group) const { DeviceGuard guard(impl_->device, impl_->available); @@ -2048,8 +2058,7 @@ void RotationScaleMergeGPU::GetFulls(int32_t *h, int32_t *k, int32_t *l, float * auto dl = [&](void *dst, const void *src, size_t bytes) { CopyAndWait(dst, src, bytes, cudaMemcpyDeviceToHost, impl_->s(), "download fulls"); }; - dl(h, dd.f_h.get(), n * sizeof(int32_t)); dl(k, dd.f_k.get(), n * sizeof(int32_t)); - dl(l, dd.f_l.get(), n * sizeof(int32_t)); dl(frame, dd.f_frame.get(), n * sizeof(int32_t)); + dl(run, dd.f_run.get(), n * sizeof(int32_t)); dl(frame, dd.f_frame.get(), n * sizeof(int32_t)); dl(group, dd.f_group.get(), n * sizeof(int32_t)); dl(I, dd.f_I.get(), n * sizeof(float)); dl(sigma, dd.f_sigma.get(), n * sizeof(float)); dl(d, dd.f_d.get(), n * sizeof(float)); dl(image_number, dd.f_img.get(), n * sizeof(float)); @@ -2066,8 +2075,6 @@ void RotationScaleMergeGPU::ResetFullsScale() { // Unity model: partiality/prescaling_corr/zeta = 1 so coeff = mean; corr starts at 1. FillKernel<<s()>>>(d.f_corr.get(), nf, 1.0f); CudaCheck(cudaGetLastError(), "FillKernel launch"); - FillKernel<<s()>>>(d.f_partiality.get(), nf, 1.0f); - CudaCheck(cudaGetLastError(), "FillKernel launch"); FullsRecordKernel<<s()>>>(nf, d.f_I.get(), d.f_sigma.get(), d.f_obs.get()); CudaCheck(cudaGetLastError(), "FullsRecordKernel launch"); FillKernel<<s()>>>(d.f_zeta.get(), nf, 1.0f); @@ -2105,13 +2112,12 @@ void RotationScaleMergeGPU::GetFullsCorr(float *corr) const { cudaMemcpyDeviceToHost, impl_->s(), "download f_corr"); } -void RotationScaleMergeGPU::GetFullsPxPy(float *px, float *py) const { +void RotationScaleMergeGPU::GetFullsPeak(int32_t *peak) const { DeviceGuard guard(impl_->device, impl_->available); const auto &d = *impl_; if (d.n_fulls == 0) return; - const size_t bytes = size_t(d.n_fulls) * sizeof(float); - CopyAndWait(px, d.f_px.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_px"); - CopyAndWait(py, d.f_py.get(), bytes, cudaMemcpyDeviceToHost, impl_->s(), "download f_py"); + CopyAndWait(peak, d.f_peak.get(), size_t(d.n_fulls) * sizeof(int32_t), cudaMemcpyDeviceToHost, impl_->s(), + "download f_peak"); } void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const { @@ -2128,7 +2134,7 @@ void RotationScaleMergeGPU::GetFullsVariance(float *var_bkg, float *var_per_I, f void RotationScaleMergeGPU::FullsDownloaded() { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; - d.f_img = {}; d.f_px = {}; d.f_py = {}; d.f_clipped = {}; + d.f_img = {}; d.f_peak = {}; d.f_clipped = {}; d.f_obs = {}; d.f_zeta = {}; d.f_frame_perm = {}; d.f_frame_start = {}; d.f_frame_count = {}; d.group_mean = {}; } @@ -2456,19 +2462,10 @@ namespace { return (u & 0x80000000u) ? ~u : (u | 0x80000000u); } - struct IngestArrays { - ObsRecord *obs; - float *zeta, *corr, *bkg, *var_bkg, *image_number, *d, *px, *py, *dphi; - int32_t *frame; - uint8_t *on_ice, *clipped; - uint64_t *key; - }; - // One block per BLK reflections: the block reads its reflections as whole words (coalesced) into - // shared memory and each thread takes its own apart. The deterministic correction is the product of - // its three factors, in the host's order; the frame's flux is divided out once it is known. - __global__ void IngestExtractKernel(int n, int base, const Reflection *__restrict__ refl, IngestArrays a) { - __shared__ uint32_t buf[BLK * REFL_WORDS]; + // shared memory, and each thread then has its own at hand. + __device__ __forceinline__ const Reflection *StageBlock(int n, const Reflection *__restrict__ refl, + uint32_t *buf) { const int first = blockIdx.x * BLK; const int cnt = min(BLK, n - first); const uint32_t *src = reinterpret_cast(refl + first); @@ -2476,9 +2473,48 @@ namespace { buf[w] = src[w]; __syncthreads(); if (int(threadIdx.x) >= cnt) + return nullptr; + return reinterpret_cast(buf + threadIdx.x * REFL_WORDS); + } + + // What the sort and the cut read of every reflection. + struct IngestKeys { + float *image_number, *d, *bkg; + uint64_t *key; + }; + __global__ void IngestKeyKernel(int n, int base, const Reflection *__restrict__ refl, IngestKeys a) { + __shared__ uint32_t buf[BLK * REFL_WORDS]; + const Reflection *r = StageBlock(n, refl, buf); + if (!r) + return; + const int i = base + blockIdx.x * BLK + threadIdx.x; + a.image_number[i] = r->image_number; + a.d[i] = r->d; + a.bkg[i] = r->bkg; + a.key[i] = RawHklKey(r->h, r->k, r->l); + } + + struct IngestArrays { + ObsRecord *obs; + float *zeta, *corr, *bkg, *var_bkg, *image_number, *d, *dphi; + uint8_t *on_ice, *clipped; + const int32_t *new_idx; // the cut's place of each reflection (an exclusive prefix over the kept) + int n, n_keep; // reflections and observations + }; + + // The observations, each at the place the cut gave its reflection (a reflection the cut dropped has + // the place of the next kept one). The deterministic correction is the product of its three + // factors, in the host's order; the frame's flux is divided out once the frames are laid out. + __global__ void IngestExtractKernel(int cnt, int base, const Reflection *__restrict__ refl, IngestArrays a) { + __shared__ uint32_t buf[BLK * REFL_WORDS]; + const Reflection *rp = StageBlock(cnt, refl, buf); + if (!rp) + return; + const Reflection &r = *rp; + const int src = base + blockIdx.x * BLK + threadIdx.x; + const int i = a.new_idx[src]; + if ((src + 1 < a.n ? a.new_idx[src + 1] : a.n_keep) == i) return; - const Reflection &r = *reinterpret_cast(buf + threadIdx.x * REFL_WORDS); - const int i = base + first + threadIdx.x; a.obs[i] = ObsRecord{r.I, r.sigma, r.partiality, __fmul_rn(__fmul_rn(r.prescaling_corr, r.qe_corr), r.flight_corr)}; a.zeta[i] = r.zeta; @@ -2487,12 +2523,9 @@ namespace { a.var_bkg[i] = r.var_bkg; a.image_number[i] = r.image_number; a.d[i] = r.d; - a.px[i] = r.predicted_x; - a.py[i] = r.predicted_y; a.dphi[i] = r.delta_phi_deg; a.on_ice[i] = r.on_ice_ring ? 1 : 0; a.clipped[i] = (r.clipped ? 1 : 0) | (r.overloaded ? 2 : 0); // OBS_CLIPPED | OBS_OVERLOADED - a.key[i] = RawHklKey(r.h, r.k, r.l); } __global__ void IngestFrameKernel(int n_frames, const int32_t *__restrict__ start, int32_t *__restrict__ frame) { @@ -2606,12 +2639,6 @@ namespace { kept[r] = keep_run[r] ? count[r] : 0; } - template - __global__ void IngestCompactKernel(int n, const int32_t *__restrict__ keep, const int32_t *__restrict__ new_idx, - const T *__restrict__ in, T *__restrict__ out) { - for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x) - if (keep[i]) out[new_idx[i]] = in[i]; - } int Blocks(int n) { return std::max(1, std::min(65535, (n + BLK - 1) / BLK)); } @@ -2796,83 +2823,85 @@ namespace { } } -void RotationScaleMergeGPU::IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, - size_t nthreads) { - DeviceGuard guard(impl_->device, impl_->available); - auto &d = *impl_; - d.n_frames = n_frames; - d.in_src_start.assign(n_frames + 1, 0); +template +void RotationScaleMergeGPU::Impl::StageReflections(const Reflection *const *refl, const int32_t *count, + size_t nthreads, Take take) { int max_count = 0; - for (int o = 0; o < n_frames; ++o) { - d.in_src_start[o + 1] = d.in_src_start[o] + count[o]; + for (int o = 0; o < n_frames; ++o) max_count = std::max(max_count, count[o]); - } - const int n = d.in_src_start[n_frames]; - d.n_obs = n; - size_t free_bytes = 0, total_bytes = 0; - CudaCheck(cudaMemGetInfo(&free_bytes, &total_bytes), "memory info"); - if (size_t(n) * OBS_BYTES > total_bytes) - ThrowTooLargeForGPU(n); - d.AllocObservations(n); - d.in_key = d.Alloc(std::max(1, n)); - d.in_dphi = d.Alloc(std::max(1, n)); - - CudaDevicePtr src_start; - d.Upload(src_start, d.in_src_start.data(), n_frames + 1); - IngestFrameKernel<<>>(n_frames, src_start.get(), - d.frame.get()); - CudaCheck(cudaGetLastError(), "IngestFrameKernel launch"); - - // Through two pinned buffers of whole frames, each filled on all threads while the other is on - // its way to the card and taken apart there. Where the card reads host memory directly - // (gpu_direct_host_read()), the kernel takes the pinned buffer apart in place, without the upload. + // Where the card reads host memory directly (gpu_direct_host_read()), the kernel takes the pinned + // buffer apart in place, without the upload. constexpr size_t STAGE_BYTES = 64 << 20; const bool direct = gpu_direct_host_read(); const int cap = std::max(static_cast(STAGE_BYTES / sizeof(Reflection)), max_count); CudaHostPtr stage[2] = {CudaHostPtr(cap), CudaHostPtr(cap)}; CudaDevicePtr staged[2]; if (!direct) { - staged[0] = d.Alloc(cap); - staged[1] = d.Alloc(cap); + staged[0] = Alloc(cap); + staged[1] = Alloc(cap); } CudaEvent copied[2]; bool in_flight[2] = {false, false}; - IngestArrays a{d.obs.get(), d.zeta.get(), d.corr.get(), d.bkg.get(), d.var_bkg.get(), d.image_number.get(), - d.d_obs.get(), d.px_obs.get(), d.py_obs.get(), d.in_dphi.get(), d.frame.get(), - d.on_ice.get(), d.clipped.get(), d.in_key.get()}; int b = 0; for (int f0 = 0; f0 < n_frames; ) { int f1 = f0 + 1; - while (f1 < n_frames && d.in_src_start[f1 + 1] - d.in_src_start[f0] <= cap) + while (f1 < n_frames && in_src_start[f1 + 1] - in_src_start[f0] <= cap) ++f1; - const int base = d.in_src_start[f0]; - const int cnt = d.in_src_start[f1] - base; + const int base = in_src_start[f0]; + const int cnt = in_src_start[f1] - base; if (in_flight[b]) CudaCheck(cudaEventSynchronize(copied[b]), "staging buffer"); Reflection *dst = stage[b].get(); ParallelFor(f1 - f0, nthreads, [&](int fo) { const int o = f0 + fo; - std::copy(refl[o], refl[o] + count[o], dst + (d.in_src_start[o] - base)); + std::copy(refl[o], refl[o] + count[o], dst + (in_src_start[o] - base)); }); if (cnt > 0) { const Reflection *src = dst; if (!direct) { CudaCheck(cudaMemcpyAsync(staged[b].get(), dst, size_t(cnt) * sizeof(Reflection), - cudaMemcpyHostToDevice, d.s()), "upload reflections"); - CudaCheck(cudaEventRecord(copied[b], d.s()), "upload reflections"); + cudaMemcpyHostToDevice, s()), "upload reflections"); + CudaCheck(cudaEventRecord(copied[b], s()), "upload reflections"); src = staged[b].get(); } in_flight[b] = true; - IngestExtractKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); - CudaCheck(cudaGetLastError(), "IngestExtractKernel launch"); + take(cnt, base, src); // Read in place, the buffer is free again only once the kernel has run. if (direct) - CudaCheck(cudaEventRecord(copied[b], d.s()), "extract reflections"); + CudaCheck(cudaEventRecord(copied[b], s()), "extract reflections"); } b ^= 1; f0 = f1; } - CudaCheck(cudaStreamSynchronize(d.s()), "upload reflections"); + CudaCheck(cudaStreamSynchronize(s()), "upload reflections"); +} + +void RotationScaleMergeGPU::IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, + size_t nthreads) { + DeviceGuard guard(impl_->device, impl_->available); + auto &d = *impl_; + d.n_frames = n_frames; + d.in_src_start.assign(n_frames + 1, 0); + for (int o = 0; o < n_frames; ++o) + d.in_src_start[o + 1] = d.in_src_start[o] + count[o]; + const int n = d.in_src_start[n_frames]; + d.n_obs = n; + size_t free_bytes = 0, total_bytes = 0; + CudaCheck(cudaMemGetInfo(&free_bytes, &total_bytes), "memory info"); + if (size_t(n) * OBS_BYTES > total_bytes) + ThrowTooLargeForGPU(n); + // See AllocObservations: these are the first of the scaling's buffers. + cuda_trim_allocation_pool(); + // Only what the sort and the cut read; the observations are taken from the reflections once the cut + // has decided which they are (IngestObservations), so the two are never on the card together. + const int na = std::max(1, n); + d.in_img = d.Alloc(na); d.in_d = d.Alloc(na); d.in_bkg = d.Alloc(na); + d.in_key = d.Alloc(na); + const IngestKeys a{d.in_img.get(), d.in_d.get(), d.in_bkg.get(), d.in_key.get()}; + d.StageReflections(refl, count, nthreads, [&](int cnt, int base, const Reflection *src) { + IngestKeyKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); + CudaCheck(cudaGetLastError(), "IngestKeyKernel launch"); + }); } void RotationScaleMergeGPU::IngestMeanBackground(double *mean_bkg) { @@ -2881,26 +2910,19 @@ void RotationScaleMergeGPU::IngestMeanBackground(double *mean_bkg) { CudaDevicePtr start; d.Upload(start, d.in_src_start.data(), d.n_frames + 1); CudaDevicePtr mean = d.Alloc(std::max(1, d.n_frames)); - IngestMeanBkgKernel<<>>(d.n_frames, start.get(), d.bkg.get(), mean.get()); + IngestMeanBkgKernel<<>>(d.n_frames, start.get(), d.in_bkg.get(), mean.get()); CudaCheck(cudaGetLastError(), "IngestMeanBkgKernel launch"); if (d.n_frames > 0) CopyAndWait(mean_bkg, mean.get(), size_t(d.n_frames) * sizeof(double), cudaMemcpyDeviceToHost, d.s(), "download mean background"); + d.in_bkg = {}; } -int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional d_min, - std::optional d_max) { +int RotationScaleMergeGPU::IngestSort(std::optional d_min, std::optional d_max) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; const int n = d.n_obs; const int na = std::max(1, n); - { - CudaDevicePtr flux; - d.Upload(flux, frame_flux, d.n_frames); - IngestFluxKernel<<>>(n, d.frame.get(), flux.get(), d.obs.get()); - CudaCheck(cudaGetLastError(), "IngestFluxKernel launch"); - CudaCheck(cudaStreamSynchronize(d.s()), "flux"); - } // (h, k, l, image_number, index) is a total order, so two stable radix passes - by the frame position // with the index as the value, then by the raw hkl - give the one permutation std::sort gives on the @@ -2926,10 +2948,11 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional ikey = d.Alloc(na), ikey_sorted = d.Alloc(na); CudaDevicePtr idx = d.Alloc(na); - IngestImageKeyKernel<<>>(n, d.image_number.get(), ikey.get(), idx.get()); + IngestImageKeyKernel<<>>(n, d.in_img.get(), ikey.get(), idx.get()); CudaCheck(cudaGetLastError(), "IngestImageKeyKernel launch"); if (n > 0) sort_pairs(ikey, ikey_sorted, idx, by_image, 32); } + d.in_img = {}; CudaDevicePtr key = d.Alloc(na); IngestGatherKeyKernel<<>>(n, by_image.get(), d.in_key.get(), key.get()); CudaCheck(cudaGetLastError(), "IngestGatherKeyKernel launch"); @@ -2984,7 +3007,7 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional keep_run = d.Alloc(nra); const bool limited = d_min.has_value() || d_max.has_value(); IngestRunKernel<<>>(n_runs, run_start.get(), run_count.get(), perm.get(), - d.d_obs.get(), d_min.has_value(), d_min.value_or(0.0), + d.in_d.get(), d_min.has_value(), d_min.value_or(0.0), d_max.has_value(), d_max.value_or(0.0), run_d.get(), keep_run.get()); CudaCheck(cudaGetLastError(), "IngestRunKernel launch"); @@ -3043,6 +3066,7 @@ int RotationScaleMergeGPU::IngestSort(const double *frame_flux, std::optional; - CudaDevicePtr out = d.Alloc(nka); - IngestCompactKernel<<>>(n, keep.get(), new_idx.get(), field.get(), out.get()); - CudaCheck(cudaGetLastError(), "IngestCompactKernel launch"); - CudaCheck(cudaStreamSynchronize(d.s()), "cut"); - field = std::move(out); - }; - cut(d.obs); cut(d.zeta); cut(d.corr); cut(d.bkg); cut(d.var_bkg); cut(d.image_number); cut(d.d_obs); - cut(d.px_obs); cut(d.py_obs); cut(d.frame); cut(d.on_ice); cut(d.clipped); cut(d.in_dphi); - } + // Where each observation goes, for IngestObservations to take them from the reflections. + d.in_new_idx = std::move(new_idx); + d.in_n_keep = n_keep; d.n_obs = n_keep; d.AllocFrames(L.frame_start.data(), L.frame_count.data()); @@ -3094,6 +3108,51 @@ void RotationScaleMergeGPU::IngestGetLayout(IngestLayout &out) const { out = impl_->in_layout; } +void RotationScaleMergeGPU::IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, + const double *frame_flux, float *px, float *py) { + DeviceGuard guard(impl_->device, impl_->available); + auto &d = *impl_; + const int n_keep = d.in_n_keep; + d.AllocObservations(n_keep); + d.in_dphi = d.Alloc(std::max(1, n_keep)); + const IngestArrays a{d.obs.get(), d.zeta.get(), d.corr.get(), d.bkg.get(), d.var_bkg.get(), d.image_number.get(), + d.d_obs.get(), d.in_dphi.get(), d.on_ice.get(), + d.clipped.get(), d.in_new_idx.get(), d.in_src_start[d.n_frames], n_keep}; + d.StageReflections(refl, count, nthreads, [&](int cnt, int base, const Reflection *src) { + IngestExtractKernel<<<(cnt + BLK - 1) / BLK, BLK, 0, d.s()>>>(cnt, base, src, a); + CudaCheck(cudaGetLastError(), "IngestExtractKernel launch"); + }); + // The predicted positions, which only the fulls' and only on the host read, stay there. + const int n = d.in_src_start[d.n_frames]; + std::vector new_idx(n); + if (n > 0) + CopyAndWait(new_idx.data(), d.in_new_idx.get(), size_t(n) * sizeof(int32_t), cudaMemcpyDeviceToHost, d.s(), + "download places"); + d.in_new_idx = {}; + ParallelFor(d.n_frames, nthreads, [&](int o) { + for (int j = 0; j < count[o]; ++j) { + const int src = d.in_src_start[o] + j, i = new_idx[src]; + if ((src + 1 < n ? new_idx[src + 1] : n_keep) == i) continue; + px[i] = refl[o][j].predicted_x; + py[i] = refl[o][j].predicted_y; + } + }); + + // The kept observations are in frame order, as the reflections were. + std::vector start(d.n_frames + 1, n_keep); + for (int o = 0; o < d.n_frames; ++o) start[o] = d.in_layout.frame_start[o]; + CudaDevicePtr d_start; + d.Upload(d_start, start.data(), d.n_frames + 1); + IngestFrameKernel<<>>(d.n_frames, d_start.get(), + d.frame.get()); + CudaCheck(cudaGetLastError(), "IngestFrameKernel launch"); + CudaDevicePtr flux; + d.Upload(flux, frame_flux, d.n_frames); + IngestFluxKernel<<>>(n_keep, d.frame.get(), flux.get(), d.obs.get()); + CudaCheck(cudaGetLastError(), "IngestFluxKernel launch"); + CudaCheck(cudaStreamSynchronize(d.s()), "flux"); +} + int64_t RotationScaleMergeGPU::IngestDeltaPhi(const DeltaPhiFrame &g) { DeviceGuard guard(impl_->device, impl_->available); auto &d = *impl_; @@ -3173,14 +3232,19 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector &key_in, CudaDevicePtr &key_out, + CudaDevicePtr &val_in, CudaDevicePtr &val_out) { + if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, val_in.get()); + CudaCheck(cudaGetLastError(), "IotaKernel launch"); + cub::DoubleBuffer keys(key_in.get(), key_out.get()), values(val_in.get(), val_out.get()); size_t bytes = 0; - CudaCheck(cub::DeviceRadixSort::SortPairs(nullptr, bytes, key_in, key_out, val_in, val_out, nf, 0, 32, d.s()), - "sort temp size"); + CudaCheck(cub::DeviceRadixSort::SortPairs(nullptr, bytes, keys, values, nf, 0, 32, d.s()), "sort temp size"); CudaDevicePtr temp = d.Alloc(std::max(1, bytes)); - CudaCheck(cub::DeviceRadixSort::SortPairs(temp.get(), bytes, key_in, key_out, val_in, val_out, nf, 0, 32, d.s()), - "sort"); + CudaCheck(cub::DeviceRadixSort::SortPairs(temp.get(), bytes, keys, values, nf, 0, 32, d.s()), "sort"); CudaCheck(cudaStreamSynchronize(d.s()), "sort"); + if (values.Current() != val_out.get()) std::swap(val_in, val_out); }; // The counts come back for the starts, a prefix the host forms in order. const auto count = [&](const int32_t *bin, int n_bins, CudaDevicePtr &out, std::vector &host) { @@ -3199,8 +3263,6 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector idx = d.Alloc(nfa), key = d.Alloc(nfa), sorted = d.Alloc(nfa); - if (nf > 0) IotaKernel<<<(nf + BLK - 1) / BLK, BLK, 0, d.s()>>>(nf, idx.get()); - CudaCheck(cudaGetLastError(), "IotaKernel launch"); std::vector fcount; count(d.f_frame.get(), nfr, d.f_frame_count, fcount); @@ -3208,7 +3270,11 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector fstart = prefix(fcount); d.Upload(d.f_frame_start, fstart.data(), nfr); d.f_frame_perm = d.Alloc(nfa); - if (nf > 0) stable_sort(d.f_frame.get(), key.get(), idx.get(), d.f_frame_perm.get()); + if (nf > 0) { + CopyAndWait(key.get(), d.f_frame.get(), size_t(nf) * sizeof(int32_t), cudaMemcpyDeviceToDevice, d.s(), + "frame key"); + stable_sort(key, sorted, idx, d.f_frame_perm); + } std::vector gcount; count(d.f_group.get(), ng, d.f_gcount, gcount); @@ -3218,7 +3284,7 @@ void RotationScaleMergeGPU::BuildFullsCSRs(int32_t *frame_count, std::vector>>(nf, d.f_group.get(), key.get()); CudaCheck(cudaGetLastError(), "FullGroupKeyKernel launch"); CudaDevicePtr by_group = d.Alloc(nfa); - if (nf > 0) stable_sort(key.get(), sorted.get(), idx.get(), by_group.get()); + if (nf > 0) stable_sort(key, sorted, idx, by_group); d.f_gperm = d.Alloc(std::max(1, n_grouped)); gperm.resize(n_grouped); if (n_grouped > 0) { diff --git a/image_analysis/scale_merge/RotationScaleMergeGPU.h b/image_analysis/scale_merge/RotationScaleMergeGPU.h index b5d0c64b8..054439b13 100644 --- a/image_analysis/scale_merge/RotationScaleMergeGPU.h +++ b/image_analysis/scale_merge/RotationScaleMergeGPU.h @@ -40,12 +40,12 @@ public: // One slice ([offset, offset+count) of n_obs) of one immutable per-observation field, so the // host can stage a bounded chunk at a time out of its AoS record: staged whole and all at once, // the fourteen arrays are more than half the entire observation payload held live again. - // Corr0 seeds the resident, mutable corr (SetCorr refreshes it). Bkg/VarBkg/ImageNumber/D/Px/Py + // Corr0 seeds the resident, mutable corr (SetCorr refreshes it). Bkg/VarBkg/ImageNumber/D // are the combine inputs (image-local background, fractional frame position for event - // contiguity, resolution, and the predicted detector position carried to the full for the - // absorption surface). + // contiguity and resolution; the predicted detector position carried to the full for the + // absorption surface stays on the host, GetFullsPeak). enum class ObsField { I, Sigma, PrescalingCorr, Partiality, Zeta, Corr0, - Bkg, VarBkg, ImageNumber, D, Px, Py }; + Bkg, VarBkg, ImageNumber, D }; void SetObsField(ObsField f, int offset, int count, const float *v); void SetObsFrame(int offset, int count, const int32_t *frame); void SetObsOnIce(int offset, int count, const uint8_t *on_ice); @@ -60,19 +60,20 @@ public: // --- the resident ingest (RotationScaleMerge::IngestDevice with no observation array on the host) --- // The steps below build, from the reflections, exactly the observations and raw-hkl runs that // RotationScaleMerge::BuildInRangeObservations builds on the host, and leave them where the uploads - // above and SetRawRuns would: the reflections cross the bus once, and are sorted and cut on the card. + // above and SetRawRuns would: the reflections are sorted and cut on the card, and cross the bus + // twice - once for what the sort and the cut read, once for the observations the cut kept. // The reflections of n_frames frames, frame o's count[o] of them at refl[o], uploaded in frame order; - // every field the ingest and the merge read is kept per reflection. + // what the sort, the cut and the background below read is kept per reflection. void IngestUpload(int n_frames, const Reflection *const *refl, const int32_t *count, size_t nthreads); // Each frame's mean background over its reflections with a finite, positive one, summed in their // order (NaN where there is none): the incident-flux meter. Length n_frames. void IngestMeanBackground(double *mean_bkg); - // Divide every reflection's deterministic correction by its frame's flux (length n_frames); order the - // reflections by (raw h,k,l, image_number, index) and cut that order into raw-hkl runs, each with the - // resolution of its first member that has a usable one; and drop the runs outside [d_min, d_max] - // (either may be absent) with their reflections. Returns how many observations are left. - int IngestSort(const double *frame_flux, std::optional d_min, std::optional d_max); + // Order the reflections by (raw h,k,l, image_number, index) and cut that order into raw-hkl runs, + // each with the resolution of its first member that has a usable one; and drop the runs outside + // [d_min, d_max] (either may be absent) with their reflections. Returns how many observations are + // left. Nothing per observation is on the card yet (IngestObservations). + int IngestSort(std::optional d_min, std::optional d_max); // The layout IngestSort left: the reflection and run counts before the resolution cut, the // observations' frame CSR, and the runs that survived it. struct IngestLayout { @@ -83,6 +84,12 @@ public: std::vector run_d; }; void IngestGetLayout(IngestLayout &out) const; + // The observations IngestSort kept, taken from the same reflections (as IngestUpload had them), in + // its order, each with its deterministic correction divided by its frame's flux (length n_frames). + // px/py (length IngestSort's count) receive each observation's predicted position, which is not + // kept on the card (GetFullsPeak). + void IngestObservations(const Reflection *const *refl, const int32_t *count, size_t nthreads, + const double *frame_flux, float *px, float *py); // The ingest's geometry smoothing, on the observations IngestSort left (RotationScaleMerge::IngestOnDevice): // each delta_phi re-solved against its frame's smoothed lattice where the frame has one - recip holds // a*, b*, c* per frame (9 floats), ok which frames, and the goniometer frame as SmoothedGeometry @@ -229,15 +236,16 @@ public: int Combine(const int32_t *rawrun_group, double min_partiality, double capture_uncertainty_coeff, double min_captured_fraction, float max_frame_gap, int64_t &n_overloaded); - // Download the combined fulls SoA (length = Combine()'s return). The working corr is downloaded - // separately by GetFullsCorr (it is only meaningful after the fulls scaling; otherwise the caller sets it). - void GetFulls(int32_t *h, int32_t *k, int32_t *l, float *I, float *sigma, float *d, + // Download the combined fulls SoA (length = Combine()'s return); run is the raw-hkl run each full is + // of. The working corr is downloaded separately by GetFullsCorr (it is only meaningful after the + // fulls scaling; otherwise the caller sets it). + void GetFulls(int32_t *run, float *I, float *sigma, float *d, float *image_number, int32_t *frame, uint8_t *on_ice, uint8_t *clipped, int32_t *group) const; - // Download the fulls' predicted detector position (peak partial's px/py), for the host absorption - // surface. Length = n_fulls. - void GetFullsPxPy(float *px, float *py) const; + // Download the observation each full's predicted detector position is that of (its peak partial), + // for the host absorption surface; the positions themselves stay on the host. Length = n_fulls. + void GetFullsPeak(int32_t *peak) const; // Download the fulls' variance model, var(I) = var_bkg + var_per_I * I + (capture * I)^2. Length = n_fulls. void GetFullsVariance(float *var_bkg, float *var_per_I, float *capture) const; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index e39a262e7..e9ce40980 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1162,7 +1162,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // free: the answer must not depend on the moment. What else stays on the card meanwhile - the // context, the run's indexers - was measured at up to 2.3 GB; a speculative pass running beside // holds more, but a merge that does not fit beside it waits for it to end. On a 16 GB card that - // is the data sets up to about 70 M observations at high multiplicity; the large ones merge + // is the data sets up to about 90 M observations at high multiplicity; the large ones merge // their arms one after the other. if (!geometry_prepass && !experiment_.GetGemmiSpaceGroup().has_value() && rsm->GetSearchMinZeta() > 0.0 && cards > 0 && rsm->OnGPU()) { diff --git a/tests/RotationIngestGPUTest.cpp b/tests/RotationIngestGPUTest.cpp index bb78e8e48..7b8588480 100644 --- a/tests/RotationIngestGPUTest.cpp +++ b/tests/RotationIngestGPUTest.cpp @@ -97,7 +97,7 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") { std::vector flux(n_frames); for (int f = 0; f < n_frames; ++f) flux[f] = 0.9 + 0.01 * f; - const int n_keep = gpu.IngestSort(flux.data(), d_min, d_max); + const int n_keep = gpu.IngestSort(d_min, d_max); // The host ingest: one total order, its raw-hkl runs, each with its first usable resolution, and // the runs the range keeps with their observations. @@ -160,6 +160,14 @@ TEST_CASE("RotationIngestGPU_SortAndCutMatchHost", "[RotationScale][gpu]") { s += c; } + std::vector px(n_keep), py(n_keep); + gpu.IngestObservations(ptr.data(), count.data(), 4, flux.data(), px.data(), py.data()); + for (int i = 0; i < n; ++i) + if (keep[i]) { + CHECK(px[new_idx[i]] == all[i]->predicted_x); + CHECK(py[new_idx[i]] == all[i]->predicted_y); + } + // The rocking events' spans, walked over each run's members in the sorted order, and corr as // integrated in the observations' order (no frame has a mosaicity to recompute partiality from). constexpr float GAP = 1.0f;