GPU: Add sorting of tracks of attached compressed clusters in deterministic mode

Author: davidrohr

GPU/GPUTracking/DataCompression/GPUTPCCompression.cxx

@@ -140,7 +140,7 @@ void GPUTPCCompression::DumpCompressedClusters(std::ostream& out)
   for (uint32_t i = 0; i < NSECTORS; i++) {
     out << "Sector " << i << ": ";
     for (uint32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
-      out << O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j] << ", ";
+      out << (O.nSliceRowClusters ? O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j] : 0) << ", ";
     }
     out << "\n";
   }
@@ -153,18 +153,20 @@ void GPUTPCCompression::DumpCompressedClusters(std::ostream& out)
   }
   out << "\n\nUnattached Clusters\n";
   uint32_t offset = 0;
-  for (uint32_t i = 0; i < NSECTORS; i++) {
-    for (uint32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
-      out << "Sector " << i << " Row " << j << ": ";
-      for (uint32_t k = 0; k < O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j]; k++) {
-        if (k && k % 10 == 0) {
-          out << "\n    ";
+  if (O.nSliceRowClusters) {
+    for (uint32_t i = 0; i < NSECTORS; i++) {
+      for (uint32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
+        out << "Sector " << i << " Row " << j << ": ";
+        for (uint32_t k = 0; k < O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j]; k++) {
+          if (k && k % 10 == 0) {
+            out << "\n    ";
+          }
+          const uint32_t l = k + offset;
+          out << "[" << (uint32_t)O.qTotU[l] << ", " << (uint32_t)O.qMaxU[l] << ", " << (uint32_t)O.flagsU[l] << ", " << (int32_t)O.padDiffU[l] << ", " << (int32_t)O.timeDiffU[l] << ", " << (uint32_t)O.sigmaPadU[l] << ", " << (uint32_t)O.sigmaTimeU[l] << "] ";
         }
-        const uint32_t l = k + offset;
-        out << "[" << (uint32_t)O.qTotU[l] << ", " << (uint32_t)O.qMaxU[l] << ", " << (uint32_t)O.flagsU[l] << ", " << (int32_t)O.padDiffU[l] << ", " << (int32_t)O.timeDiffU[l] << ", " << (uint32_t)O.sigmaPadU[l] << ", " << (uint32_t)O.sigmaTimeU[l] << "] ";
+        offset += O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j];
+        out << "\n";
       }
-      offset += O.nSliceRowClusters[i * GPUCA_ROW_COUNT + j];
-      out << "\n";
     }
   }
   out << "\n\nAttached Clusters\n";
@@ -175,7 +177,7 @@ void GPUTPCCompression::DumpCompressedClusters(std::ostream& out)
       if (k && k % 10 == 0) {
         out << "\n    ";
       }
-      const uint32_t l1 = k + offset, l2 = k + offset - i;
+      const uint32_t l1 = offset + k, l2 = offset - i + k - 1;
       out << "[";
       if (k) {
         out << (int32_t)O.rowDiffA[l2] << ", " << (int32_t)O.sliceLegDiffA[l2] << ", " << (uint32_t)O.padResA[l2] << ", " << (uint32_t)O.timeResA[l2] << ", ";

GPU/GPUTracking/Global/GPUChainTracking.h

@@ -235,6 +235,7 @@ class GPUChainTracking : public GPUChain
   void PrintDebugOutput();
   void PrintOutputStat();
   static void DumpClusters(std::ostream& out, const o2::tpc::ClusterNativeAccess* clusters);
+  static void DebugSortCompressedClusters(o2::tpc::CompressedClustersFlat* cls);
 
   bool ValidateSteps();
   bool ValidateSettings();

GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx

@@ -203,6 +203,10 @@ int32_t GPUChainTracking::RunTPCCompression()
     ((GPUChainTracking*)GetNextChainInQueue())->mRec->BlockStackedMemory(mRec);
   }
   mRec->PopNonPersistentMemory(RecoStep::TPCCompression, qStr2Tag("TPCCOMPR"));
+  if (GetProcessingSettings().deterministicGPUReconstruction) {
+    SynchronizeGPU();
+    DebugSortCompressedClusters(Compressor.mOutputFlat);
+  }
   DoDebugAndDump(RecoStep::TPCCompression, GPUChainTrackingDebugFlags::TPCCompressedClusters, Compressor, &GPUTPCCompression::DumpCompressedClusters, *mDebugFile);
   return 0;
 }

GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx

@@ -20,6 +20,7 @@
 #include <map>
 #include <memory>
 #include <string>
+#include <numeric>
 
 #ifdef GPUCA_TRACKLET_CONSTRUCTOR_DO_PROFILE
 #include "bitmapfile.h"
@@ -348,3 +349,44 @@ void GPUChainTracking::DumpClusters(std::ostream& out, const o2::tpc::ClusterNat
     }
   }
 }
+
+void GPUChainTracking::DebugSortCompressedClusters(o2::tpc::CompressedClustersFlat* cls)
+{
+  o2::tpc::CompressedClusters c = *cls;
+  std::vector<uint32_t> sorted(c.nTracks), offsets(c.nTracks);
+  std::iota(sorted.begin(), sorted.end(), 0);
+  auto sorter = [&c](const auto a, const auto b) {
+    return std::tie(c.sliceA[a], c.rowA[a], c.timeA[a], c.padA[a], c.qPtA[a]) <
+           std::tie(c.sliceA[b], c.rowA[b], c.timeA[b], c.padA[b], c.qPtA[b]);
+  };
+  std::sort(sorted.begin(), sorted.end(), sorter);
+  uint32_t offset = 0;
+  for (uint32_t i = 0; i < c.nTracks; i++) {
+    offsets[i] = offset;
+    offset += c.nTrackClusters[i];
+  }
+
+  auto sortArray = [&c, &sorted, &offsets](auto* src, size_t totalSize, auto getOffset, auto getSize) {
+    auto buf = std::make_unique<std::remove_reference_t<decltype(src[0])>[]>(totalSize);
+    memcpy(buf.get(), src, totalSize * sizeof(*src));
+    uint32_t targetOffset = 0;
+    for (uint32_t i = 0; i < c.nTracks; i++) {
+      const uint32_t j = sorted[i];
+      memcpy(src + targetOffset, buf.get() + getOffset(offsets[j], j), getSize(j) * sizeof(*src));
+      targetOffset += getSize(j);
+    }
+  };
+  auto sortMultiple = [&sortArray](size_t totalSize, auto getOffset, auto getSize, auto&&... arrays) {
+    (..., sortArray(std::forward<decltype(arrays)>(arrays), totalSize, getOffset, getSize));
+  };
+  auto getFullOffset = [](uint32_t off, uint32_t ind) { return off; };
+  auto getReducedOffset = [](uint32_t off, uint32_t ind) { return off - ind; };
+  auto getIndex = [](uint32_t off, uint32_t ind) { return ind; };
+  auto getN = [&c](uint32_t j) { return c.nTrackClusters[j]; };
+  auto getN1 = [&c](uint32_t j) { return c.nTrackClusters[j] - 1; };
+  auto get1 = [](uint32_t j) { return 1; };
+
+  sortMultiple(c.nAttachedClusters, getFullOffset, getN, c.qTotA, c.qMaxA, c.flagsA, c.sigmaPadA, c.sigmaTimeA);
+  sortMultiple(c.nAttachedClustersReduced, getReducedOffset, getN1, c.rowDiffA, c.sliceLegDiffA, c.padResA, c.timeResA);
+  sortMultiple(c.nTracks, getIndex, get1, c.qPtA, c.rowA, c.sliceA, c.timeA, c.padA, c.nTrackClusters); // NOTE: This must be last, since nTrackClusters is used for handling the arrays above!
+}