From cd4616ff24b2502066c77d4ddeccb6086556983e Mon Sep 17 00:00:00 2001 From: Giulio Eulisse <10544+ktf@users.noreply.github.com> Date: Sat, 7 Nov 2026 07:27:23 +0100 Subject: [PATCH 1/4] GPU: put the orderings of the merger's sorts in a header Move the comparators of the TPC merger's sorts, and MergeLooperParam they sort, from GPUTPCGMMerger.cxx and GPUTPCGMO2Output.cxx into GPUTPCGMMergerSortComp.h, unchanged. The CUDA and HIP backends see them because they compile the kernel sources into one unit; a backend that runs these sorts in host code needs them from a header. --- GPU/GPUTracking/CMakeLists.txt | 1 + GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx | 79 +------------ .../Merger/GPUTPCGMMergerSortComp.h | 110 ++++++++++++++++++ GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx | 12 +- 4 files changed, 113 insertions(+), 89 deletions(-) create mode 100644 GPU/GPUTracking/Merger/GPUTPCGMMergerSortComp.h diff --git a/GPU/GPUTracking/CMakeLists.txt b/GPU/GPUTracking/CMakeLists.txt index 6f2e643fc98f6..b5c314fcee9d0 100644 --- a/GPU/GPUTracking/CMakeLists.txt +++ b/GPU/GPUTracking/CMakeLists.txt @@ -145,6 +145,7 @@ set(HDRS_INSTALL Global/GPUChainTrackingGetters.inc Global/GPUErrorCodes.h Merger/GPUTPCGMMergerTypes.h + Merger/GPUTPCGMMergerSortComp.h qa/GPUQAHelper.h qconfigoptions.h Refit/GPUTrackParamConvert.h diff --git a/GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx b/GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx index 267239a4bbab9..94acb7bfe9cf9 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx @@ -63,84 +63,7 @@ using namespace o2::gpu; using namespace o2::tpc; using namespace gputpcgmmergertypes; -namespace o2::gpu::internal -{ -struct MergeLooperParam { - float refz; - float x; - float y; - uint32_t id; -}; - -struct MergeBorderTracks_compMax { - GPUd() bool operator()(const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b) - { - return GPUCA_DETERMINISTIC_CODE((a.fMax != b.fMax) ? (a.fMax < b.fMax) : (a.fId < b.fId), a.fMax < b.fMax); - } -}; -struct MergeBorderTracks_compMin { - GPUd() bool operator()(const GPUTPCGMBorderRange& a, const GPUTPCGMBorderRange& b) - { - return GPUCA_DETERMINISTIC_CODE((a.fMin != b.fMin) ? (a.fMin < b.fMin) : (a.fId < b.fId), a.fMin < b.fMin); - } -}; - -struct GPUTPCGMMergerSortTracks_comp { - const GPUTPCGMMergedTrack* const mCmp; - GPUhd() GPUTPCGMMergerSortTracks_comp(GPUTPCGMMergedTrack* cmp) : mCmp(cmp) {} - GPUd() bool operator()(const int32_t aa, const int32_t bb) - { - const GPUTPCGMMergedTrack& GPUrestrict() a = mCmp[aa]; - const GPUTPCGMMergedTrack& GPUrestrict() b = mCmp[bb]; - if (a.OK() != b.OK()) { - return a.OK(); - } - if (a.CCE() != b.CCE()) { - return a.CCE() > b.CCE(); - } - GPUCA_DETERMINISTIC_CODE( // clang-format off - if (a.NClusters() != b.NClusters()) { - return a.NClusters() > b.NClusters(); - } if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) { - return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); - } if (a.GetParam().GetY() != b.GetParam().GetY()) { - return a.GetParam().GetY() > b.GetParam().GetY(); - } - return aa > bb; - , // !GPUCA_DETERMINISTIC_CODE - return a.NClusters() > b.NClusters(); - ) // clang-format on - } -}; - -struct GPUTPCGMMergerSortTracksQPt_comp { - const GPUTPCGMMergedTrack* const mCmp; - GPUhd() GPUTPCGMMergerSortTracksQPt_comp(GPUTPCGMMergedTrack* cmp) : mCmp(cmp) {} - GPUd() bool operator()(const int32_t aa, const int32_t bb) - { - const GPUTPCGMMergedTrack& GPUrestrict() a = mCmp[aa]; - const GPUTPCGMMergedTrack& GPUrestrict() b = mCmp[bb]; - GPUCA_DETERMINISTIC_CODE( // clang-format off - if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) { - return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); - } if (a.GetParam().GetY() != b.GetParam().GetY()) { - return a.GetParam().GetY() > b.GetParam().GetY(); - } - return a.GetParam().GetZ() > b.GetParam().GetZ(); - , // !GPUCA_DETERMINISTIC_CODE - return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); - ) // clang-format on - } -}; - -struct GPUTPCGMMergerMergeLoopers_comp { - GPUd() bool operator()(const MergeLooperParam& a, const MergeLooperParam& b) - { - return GPUCA_DETERMINISTIC_CODE(CAMath::Abs(a.refz) != CAMath::Abs(b.refz) ? CAMath::Abs(a.refz) < CAMath::Abs(b.refz) : a.id < b.id, CAMath::Abs(a.refz) < CAMath::Abs(b.refz)); - } -}; - -} // namespace o2::gpu::internal +#include "GPUTPCGMMergerSortComp.h" using namespace o2::gpu::internal; diff --git a/GPU/GPUTracking/Merger/GPUTPCGMMergerSortComp.h b/GPU/GPUTracking/Merger/GPUTPCGMMergerSortComp.h new file mode 100644 index 0000000000000..ea239a726a71a --- /dev/null +++ b/GPU/GPUTracking/Merger/GPUTPCGMMergerSortComp.h @@ -0,0 +1,110 @@ +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file GPUTPCGMMergerSortComp.h +/// \brief Orderings of the merger's sorts, for the kernels and for the backends that run these sorts themselves + +#ifndef GPUTPCGMMERGERSORTCOMP_H +#define GPUTPCGMMERGERSORTCOMP_H + +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" +#include "GPUTPCGMMergedTrack.h" +#include "GPUTPCGMMergerTypes.h" +#include "GPUTPCGMMerger.h" + +namespace o2::gpu::internal +{ +struct MergeLooperParam { + float refz; + float x; + float y; + uint32_t id; +}; + +struct MergeBorderTracks_compMax { + GPUd() bool operator()(const gputpcgmmergertypes::GPUTPCGMBorderRange& a, const gputpcgmmergertypes::GPUTPCGMBorderRange& b) + { + return GPUCA_DETERMINISTIC_CODE((a.fMax != b.fMax) ? (a.fMax < b.fMax) : (a.fId < b.fId), a.fMax < b.fMax); + } +}; +struct MergeBorderTracks_compMin { + GPUd() bool operator()(const gputpcgmmergertypes::GPUTPCGMBorderRange& a, const gputpcgmmergertypes::GPUTPCGMBorderRange& b) + { + return GPUCA_DETERMINISTIC_CODE((a.fMin != b.fMin) ? (a.fMin < b.fMin) : (a.fId < b.fId), a.fMin < b.fMin); + } +}; + +struct GPUTPCGMMergerSortTracks_comp { + const GPUTPCGMMergedTrack* const mCmp; + GPUhd() GPUTPCGMMergerSortTracks_comp(GPUTPCGMMergedTrack* cmp) : mCmp(cmp) {} + GPUd() bool operator()(const int32_t aa, const int32_t bb) + { + const GPUTPCGMMergedTrack& GPUrestrict() a = mCmp[aa]; + const GPUTPCGMMergedTrack& GPUrestrict() b = mCmp[bb]; + if (a.OK() != b.OK()) { + return a.OK(); + } + if (a.CCE() != b.CCE()) { + return a.CCE() > b.CCE(); + } + GPUCA_DETERMINISTIC_CODE( // clang-format off + if (a.NClusters() != b.NClusters()) { + return a.NClusters() > b.NClusters(); + } if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) { + return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); + } if (a.GetParam().GetY() != b.GetParam().GetY()) { + return a.GetParam().GetY() > b.GetParam().GetY(); + } + return aa > bb; + , // !GPUCA_DETERMINISTIC_CODE + return a.NClusters() > b.NClusters(); + ) // clang-format on + } +}; + +struct GPUTPCGMMergerSortTracksQPt_comp { + const GPUTPCGMMergedTrack* const mCmp; + GPUhd() GPUTPCGMMergerSortTracksQPt_comp(GPUTPCGMMergedTrack* cmp) : mCmp(cmp) {} + GPUd() bool operator()(const int32_t aa, const int32_t bb) + { + const GPUTPCGMMergedTrack& GPUrestrict() a = mCmp[aa]; + const GPUTPCGMMergedTrack& GPUrestrict() b = mCmp[bb]; + GPUCA_DETERMINISTIC_CODE( // clang-format off + if (CAMath::Abs(a.GetParam().GetQPt()) != CAMath::Abs(b.GetParam().GetQPt())) { + return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); + } if (a.GetParam().GetY() != b.GetParam().GetY()) { + return a.GetParam().GetY() > b.GetParam().GetY(); + } + return a.GetParam().GetZ() > b.GetParam().GetZ(); + , // !GPUCA_DETERMINISTIC_CODE + return CAMath::Abs(a.GetParam().GetQPt()) > CAMath::Abs(b.GetParam().GetQPt()); + ) // clang-format on + } +}; + +struct GPUTPCGMMergerMergeLoopers_comp { + GPUd() bool operator()(const MergeLooperParam& a, const MergeLooperParam& b) + { + return GPUCA_DETERMINISTIC_CODE(CAMath::Abs(a.refz) != CAMath::Abs(b.refz) ? CAMath::Abs(a.refz) < CAMath::Abs(b.refz) : a.id < b.id, CAMath::Abs(a.refz) < CAMath::Abs(b.refz)); + } +}; + +struct GPUTPCGMO2OutputSort_comp { + GPUd() bool operator()(const GPUTPCGMMerger::tmpSort& a, const GPUTPCGMMerger::tmpSort& b) + { + return GPUCA_DETERMINISTIC_CODE(a.y != b.y ? a.y > b.y : a.x > b.x, a.y > b.y); + } +}; + +} // namespace o2::gpu::internal + +#endif diff --git a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx b/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx index 1ba55d98f5ca1..33b658ae01285 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx @@ -33,17 +33,7 @@ using namespace o2::tpc::constants; GPUdi() static constexpr uint8_t getFlagsReject() { return GPUTPCGMMergedTrackHit::flagReject | GPUTPCGMMergedTrackHit::flagHighIncl; } -namespace o2::gpu::internal -{ - -struct GPUTPCGMO2OutputSort_comp { - GPUd() bool operator()(const GPUTPCGMMerger::tmpSort& a, const GPUTPCGMMerger::tmpSort& b) - { - return GPUCA_DETERMINISTIC_CODE(a.y != b.y ? a.y > b.y : a.x > b.x, a.y > b.y); - } -}; - -} // namespace o2::gpu::internal +#include "GPUTPCGMMergerSortComp.h" template <> GPUdii() void GPUTPCGMO2Output::Thread(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUsharedref() GPUSharedMemory& smem, processorType& GPUrestrict() merger) From 4bbec03b29626bf7152a043528d85ac2ded1fdb4 Mon Sep 17 00:00:00 2001 From: Giulio Eulisse <10544+ktf@users.noreply.github.com> Date: Sat, 7 Nov 2026 07:27:23 +0100 Subject: [PATCH 2/4] GPU: make Metal's cluster widths round as on the CPU Needs review. On a 30-collision Pb-Pb time frame, 7 of 4,416,506 clusters differed between the CPU and the Metal backend, each by one step of a packed width, sigmaTime or sigmaPad. Those widths were all within one float ulp below a quantisation step: sigma * 32 was 11.4999994, 14.4999992, 29.4999996 or 31.4999988, where the packing rounds at n + 0.5, so any difference in the last bit decides between n and n + 1. The CPU build compiles this with -ffast-math. Its arm64 code for finalize() computes 1 / Q once, the means as products with it, the variance as one fused multiply-add, fma(SS, 1 / Q, -mean * mean), and an IEEE sqrt. Metal's fast-math division and sqrt are not correctly rounded. Do the same operations on Metal explicitly, with precise::divide, fma and precise::sqrt: the pbpb clusters then match the CPU bit for bit, with no measurable cost; the rest of the Metal code keeps fast math. This matches the CPU build that the compiler produces on arm64 with clang. An x86 build, or another compiler, may choose a different sequence under -ffast-math. Writing the sequence explicitly for all backends would make the widths independent of those choices, but can change the CPU results on some platforms. --- .../TPCClusterFinder/ClusterAccumulator.cxx | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx b/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx index bcf04d2a6de1b..6a17b8ec49c2c 100644 --- a/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx +++ b/GPU/GPUTracking/TPCClusterFinder/ClusterAccumulator.cxx @@ -62,6 +62,18 @@ GPUd() void ClusterAccumulator::finalize(const CfChargePos& pos, const Charge q, { mQtot += q; +#ifdef __METAL__ + // Metal's fast-math division and sqrt are not correctly rounded, which moves + // widths that sit within an ulp of a packing step. Do the operations the + // arm64 -ffast-math CPU build does here: reciprocal, products, fused variance. + const float invQ = metal::precise::divide(1.f, mQtot); + mPadMean *= invQ; + mTimeMean *= invQ; + const float padMean2 = mPadMean * mPadMean; + const float timeMean2 = mTimeMean * mTimeMean; + mPadSigma = metal::precise::sqrt(metal::fma(mPadSigma, invQ, -padMean2)); + mTimeSigma = metal::precise::sqrt(metal::fma(mTimeSigma, invQ, -timeMean2)); +#else mPadMean /= mQtot; mTimeMean /= mQtot; mPadSigma /= mQtot; @@ -69,6 +81,7 @@ GPUd() void ClusterAccumulator::finalize(const CfChargePos& pos, const Charge q, mPadSigma = CAMath::Sqrt(mPadSigma - mPadMean * mPadMean); mTimeSigma = CAMath::Sqrt(mTimeSigma - mTimeMean * mTimeMean); +#endif Pad pad = pos.pad(); mPadMean += pad; From aebb281f6bf1cafd47d403961221626d94af23ad Mon Sep 17 00:00:00 2001 From: Giulio Eulisse <10544+ktf@users.noreply.github.com> Date: Sat, 7 Nov 2026 07:27:23 +0100 Subject: [PATCH 3/4] GPU Metal: run the merger's sorts on the host The sorts of the merged tracks by quality and by q/pT, of the looper candidates and of the O2 output tracks, which CUDA and HIP run with thrust, ran on a single GPU thread otherwise: about 730 ms per event on a 30-collision Pb-Pb time frame, at up to 140 times the CPU's time. The GPU memory is shared with the host, so sort there once the stream is drained, with the std::sort of the CPU backend, which also orders the elements that compare equal as on the CPU. The number of elements is read from the device copy of the counters, as the host copy of the O2 output count is still being transferred when the sort is launched. --- .../Base/metal/GPUReconstructionMetal.h | 4 ++ ...PUReconstructionMetalKernelsSpecialize.inc | 53 +++++++++++++++++++ 2 files changed, 57 insertions(+) diff --git a/GPU/GPUTracking/Base/metal/GPUReconstructionMetal.h b/GPU/GPUTracking/Base/metal/GPUReconstructionMetal.h index 247ae6657c01e..091fbbb415804 100644 --- a/GPU/GPUTracking/Base/metal/GPUReconstructionMetal.h +++ b/GPU/GPUTracking/Base/metal/GPUReconstructionMetal.h @@ -64,6 +64,10 @@ class GPUReconstructionMetal : public GPUReconstructionProcessing::KernelInterfa void CommitStream(int32_t stream); void DrainStream(int32_t stream); uint64_t EventValue(deviceEvent ev, bool next); + template + T* HostPointer(T* p); + template + void SortOnHost(const krnlSetupTime& _xyz, T* begin, const N& count, const S& comp); private: int32_t AddKernels(); diff --git a/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc b/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc index 5f456a7ce06bb..869ba7790526b 100644 --- a/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc +++ b/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc @@ -11,6 +11,8 @@ /// \file GPUReconstructionMetalKernelsSpecialize.inc +#include "GPUTPCGMMergerSortComp.h" + template <> inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz, void* const& ptr, uint64_t const& size) { @@ -30,3 +32,54 @@ inline void GPUReconstructionMetal::runKernelBackend(const krn RecordMarker(_xyz.z.ev, _xyz.x.stream); } } + +// The sorts that CUDA and HIP run with thrust (GPUCA_SPECIALIZE_THRUST_SORTS) +// run on a single GPU thread otherwise. The memory is shared with the host, so +// sort there, with the same std::sort as the CPU backend, which also orders +// the elements that compare equal as on the CPU. The number of elements is +// read once the stream is drained, from the device copy of the counters, +// since the host copy may still be on its way. +template +inline T* GPUReconstructionMetal::HostPointer(T* p) +{ + const size_t offset = (const char*)p - (const char*)mDeviceMemoryBase; + return offset < mDeviceMemorySize ? (T*)((char*)mInternals->mem_gpu.contents + offset) : p; +} + +template +inline void GPUReconstructionMetal::SortOnHost(const krnlSetupTime& _xyz, T* begin, const N& count, const S& comp) +{ + if (_xyz.z.evList) { + StreamWaitForEvents(_xyz.x.stream, _xyz.z.evList, _xyz.z.nEvents); + } + DrainStream(_xyz.x.stream); + T* hostBegin = HostPointer(begin); + GPUCommonAlgorithm::sortDeviceDynamic(hostBegin, hostBegin + count(), comp); + if (_xyz.z.ev) { + RecordMarker(_xyz.z.ev, _xyz.x.stream); + } +} + +template <> +inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz) +{ + SortOnHost(_xyz, mProcessorsShadow->tpcMerger.TrackOrderProcess(), [&] { return HostPointer(mProcessorsShadow->tpcMerger.Memory())->nMergedTracks; }, internal::GPUTPCGMMergerSortTracks_comp(HostPointer(mProcessorsShadow->tpcMerger.MergedTracks()))); +} + +template <> +inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz) +{ + SortOnHost(_xyz, mProcessorsShadow->tpcMerger.TrackSort(), [&] { return HostPointer(mProcessorsShadow->tpcMerger.Memory())->nMergedTracks; }, internal::GPUTPCGMMergerSortTracksQPt_comp(HostPointer(mProcessorsShadow->tpcMerger.MergedTracks()))); +} + +template <> +inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz) +{ + SortOnHost(_xyz, mProcessorsShadow->tpcMerger.LooperCandidates(), [&] { return HostPointer(mProcessorsShadow->tpcMerger.Memory())->nLooperMatchCandidates; }, internal::GPUTPCGMMergerMergeLoopers_comp()); +} + +template <> +inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz) +{ + SortOnHost(_xyz, mProcessorsShadow->tpcMerger.TrackSortO2(), [&] { return HostPointer(mProcessorsShadow->tpcMerger.Memory())->nO2Tracks; }, internal::GPUTPCGMO2OutputSort_comp()); +} From f058f6efd45b1251b32655139505eb0ecacdd32e Mon Sep 17 00:00:00 2001 From: Giulio Eulisse <10544+ktf@users.noreply.github.com> Date: Sat, 7 Nov 2026 07:27:23 +0100 Subject: [PATCH 4/4] GPU Metal: also run the border-range sorts on the host The border-range sorts of the track merging between sectors, 360 per time frame, took 312 ms per event on a 30-collision Pb-Pb time frame on a single GPU thread. On the host they take about 6 ms, plus 52 ms of waiting for the streams. --- .../metal/GPUReconstructionMetalKernelsSpecialize.inc | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc b/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc index 869ba7790526b..63b771565f80d 100644 --- a/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc +++ b/GPU/GPUTracking/Base/metal/GPUReconstructionMetalKernelsSpecialize.inc @@ -83,3 +83,13 @@ inline void GPUReconstructionMetal::runKernelBackendtpcMerger.TrackSortO2(), [&] { return HostPointer(mProcessorsShadow->tpcMerger.Memory())->nO2Tracks; }, internal::GPUTPCGMO2OutputSort_comp()); } + +template <> +inline void GPUReconstructionMetal::runKernelBackend(const krnlSetupTime& _xyz, gputpcgmmergertypes::GPUTPCGMBorderRange* const& range, int32_t const& N, int32_t const& cmpMax) +{ + if (cmpMax) { + SortOnHost(_xyz, range, [&] { return N; }, internal::MergeBorderTracks_compMax()); + } else { + SortOnHost(_xyz, range, [&] { return N; }, internal::MergeBorderTracks_compMin()); + } +}