CudaMuonSegmentFitting/src/CudaMuonHoughMaximum.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/EventData/CudaMuonHoughMaximum.hpp"
0010
0011 #include <limits>
0012 #include <stdexcept>
0013 #include <utility>
0014
0015 #include <cuda_runtime_api.h>
0016
0017 #include "CudaUtilities.hpp"
0018
0019 namespace ActsExamples {
0020
0021 CudaHoughMaximumBatch::CudaHoughMaximumBatch(size_type nBuckets,
0022                                              size_type capacityPerBucket)
0023     : m_nBuckets{nBuckets}, m_capacityPerBucket{capacityPerBucket} {
0024   if (m_nBuckets == 0u) {
0025     throw std::invalid_argument(
0026         "CudaHoughMaximumBatch requires non-zero nBuckets");
0027   }
0028
0029   if (m_capacityPerBucket == 0u || m_capacityPerBucket > 64u) {
0030     throw std::invalid_argument(
0031         "CudaHoughMaximumBatch capacityPerBucket must be between 1 and 64");
0032   }
0033
0034   if (m_nBuckets >
0035       std::numeric_limits<std::uint32_t>::max() / m_capacityPerBucket) {
0036     throw std::overflow_error(
0037         "CudaHoughMaximumBatch total capacity must fit into std::uint32_t");
0038   }
0039
0040   const size_type capacity = totalCapacity();
0041
0042   m_hostTanBeta.resize(capacity, CoordType{0.0});
0043   m_hostInterceptY.resize(capacity, CoordType{0.0});
0044
0045   m_hostHits.resize(capacity, YieldType{0.0});
0046   m_hostLayers.resize(capacity, YieldType{0.0});
0047   m_hostLayerMask.resize(capacity, LayerMask{0ull});
0048
0049   m_hostXBin.resize(capacity, 0u);
0050   m_hostYBin.resize(capacity, 0u);
0051
0052   m_hostNMaxima.resize(m_nBuckets, 0u);
0053   m_hostNAssociatedHits.resize(capacity, 0u);
0054 }
0055
0056 CudaHoughMaximumBatch::CudaHoughMaximumBatch(
0057     CudaHoughMaximumBatch&& other) noexcept
0058     : m_nBuckets{std::exchange(other.m_nBuckets, 0u)},
0059       m_capacityPerBucket{std::exchange(other.m_capacityPerBucket, 0u)},
0060       m_hostTanBeta{std::move(other.m_hostTanBeta)},
0061       m_hostInterceptY{std::move(other.m_hostInterceptY)},
0062       m_hostHits{std::move(other.m_hostHits)},
0063       m_hostLayers{std::move(other.m_hostLayers)},
0064       m_hostLayerMask{std::move(other.m_hostLayerMask)},
0065       m_hostXBin{std::move(other.m_hostXBin)},
0066       m_hostYBin{std::move(other.m_hostYBin)},
0067       m_hostNMaxima{std::move(other.m_hostNMaxima)},
0068       m_hostNAssociatedHits{std::move(other.m_hostNAssociatedHits)},
0069       m_hostAssociatedHitOffsets{std::move(other.m_hostAssociatedHitOffsets)},
0070       m_hostAssociatedHitIndices{std::move(other.m_hostAssociatedHitIndices)},
0071       m_associationMetadataOnHost{
0072           std::exchange(other.m_associationMetadataOnHost, false)},
0073       m_associationStorageAllocated{
0074           std::exchange(other.m_associationStorageAllocated, false)},
0075       m_associatedHitIndicesOnHost{
0076           std::exchange(other.m_associatedHitIndicesOnHost, false)},
0077       m_device{std::exchange(other.m_device, CudaHoughMaximumBatchArrays{})},
0078       m_onDevice{std::exchange(other.m_onDevice, false)} {}
0079
0080 CudaHoughMaximumBatch& CudaHoughMaximumBatch::operator=(
0081     CudaHoughMaximumBatch&& other) noexcept {
0082   if (this == &other) {
0083     return *this;
0084   }
0085
0086   clearDevice();
0087
0088   m_nBuckets = std::exchange(other.m_nBuckets, 0u);
0089   m_capacityPerBucket = std::exchange(other.m_capacityPerBucket, 0u);
0090
0091   m_hostTanBeta = std::move(other.m_hostTanBeta);
0092   m_hostInterceptY = std::move(other.m_hostInterceptY);
0093
0094   m_hostHits = std::move(other.m_hostHits);
0095   m_hostLayers = std::move(other.m_hostLayers);
0096   m_hostLayerMask = std::move(other.m_hostLayerMask);
0097
0098   m_hostXBin = std::move(other.m_hostXBin);
0099   m_hostYBin = std::move(other.m_hostYBin);
0100
0101   m_hostNMaxima = std::move(other.m_hostNMaxima);
0102   m_hostNAssociatedHits = std::move(other.m_hostNAssociatedHits);
0103   m_hostAssociatedHitOffsets = std::move(other.m_hostAssociatedHitOffsets);
0104   m_hostAssociatedHitIndices = std::move(other.m_hostAssociatedHitIndices);
0105
0106   m_associationMetadataOnHost =
0107       std::exchange(other.m_associationMetadataOnHost, false);
0108   m_associationStorageAllocated =
0109       std::exchange(other.m_associationStorageAllocated, false);
0110   m_associatedHitIndicesOnHost =
0111       std::exchange(other.m_associatedHitIndicesOnHost, false);
0112
0113   m_device = std::exchange(other.m_device, CudaHoughMaximumBatchArrays{});
0114   m_onDevice = std::exchange(other.m_onDevice, false);
0115
0116   return *this;
0117 }
0118
0119 CudaHoughMaximumBatch::~CudaHoughMaximumBatch() noexcept {
0120   clearDevice();
0121 }
0122
0123 CudaHoughMaximumBatch::size_type CudaHoughMaximumBatch::nMaxima(
0124     size_type bucket) const {
0125   checkBucket(bucket);
0126
0127   const size_type count = static_cast<size_type>(m_hostNMaxima[bucket]);
0128
0129   if (count > capacityPerBucket()) {
0130     throw std::runtime_error(
0131         "CudaHoughMaximumBatch contains an invalid maximum count");
0132   }
0133
0134   return count;
0135 }
0136
0137 CoordType CudaHoughMaximumBatch::tanBeta(size_type bucket,
0138                                          size_type maximum) const {
0139   checkMaximum(bucket, maximum);
0140   return m_hostTanBeta[slotIndex(bucket, maximum)];
0141 }
0142
0143 CoordType CudaHoughMaximumBatch::interceptY(size_type bucket,
0144                                             size_type maximum) const {
0145   checkMaximum(bucket, maximum);
0146   return m_hostInterceptY[slotIndex(bucket, maximum)];
0147 }
0148
0149 YieldType CudaHoughMaximumBatch::nHits(size_type bucket,
0150                                        size_type maximum) const {
0151   checkMaximum(bucket, maximum);
0152   return m_hostHits[slotIndex(bucket, maximum)];
0153 }
0154
0155 YieldType CudaHoughMaximumBatch::nLayers(size_type bucket,
0156                                          size_type maximum) const {
0157   checkMaximum(bucket, maximum);
0158   return m_hostLayers[slotIndex(bucket, maximum)];
0159 }
0160
0161 LayerMask CudaHoughMaximumBatch::layerMask(size_type bucket,
0162                                            size_type maximum) const {
0163   checkMaximum(bucket, maximum);
0164   return m_hostLayerMask[slotIndex(bucket, maximum)];
0165 }
0166
0167 CudaHoughMaximumBatch::size_type CudaHoughMaximumBatch::xBin(
0168     size_type bucket, size_type maximum) const {
0169   checkMaximum(bucket, maximum);
0170   return static_cast<size_type>(m_hostXBin[slotIndex(bucket, maximum)]);
0171 }
0172
0173 CudaHoughMaximumBatch::size_type CudaHoughMaximumBatch::yBin(
0174     size_type bucket, size_type maximum) const {
0175   checkMaximum(bucket, maximum);
0176   return static_cast<size_type>(m_hostYBin[slotIndex(bucket, maximum)]);
0177 }
0178
0179 CudaHoughMaximumBatch::size_type CudaHoughMaximumBatch::nAssociatedHits(
0180     size_type bucket, size_type maximum) const {
0181   if (!m_associationMetadataOnHost) {
0182     throw std::logic_error("Association metadata is not available on the host");
0183   }
0184
0185   checkMaximum(bucket, maximum);
0186   return m_hostNAssociatedHits[slotIndex(bucket, maximum)];
0187 }
0188
0189 std::span<const std::uint32_t> CudaHoughMaximumBatch::associatedHitIndices(
0190     size_type bucket, size_type maximum) const {
0191   if (!m_associatedHitIndicesOnHost) {
0192     throw std::logic_error(
0193         "Associated hit indices are not available on the host");
0194   }
0195
0196   checkMaximum(bucket, maximum);
0197
0198   const size_type slot = slotIndex(bucket, maximum);
0199   const size_type begin = m_hostAssociatedHitOffsets[slot];
0200   const size_type end = m_hostAssociatedHitOffsets[slot + 1u];
0201
0202   return {m_hostAssociatedHitIndices.data() + begin, end - begin};
0203 }
0204
0205 void CudaHoughMaximumBatch::moveToDevice(cudaStream_t /*stream*/) {
0206   clearDevice();
0207
0208   m_hostAssociatedHitOffsets.clear();
0209   m_hostAssociatedHitIndices.clear();
0210
0211   m_associationMetadataOnHost = false;
0212   m_associationStorageAllocated = false;
0213   m_associatedHitIndicesOnHost = false;
0214
0215   m_device.nBuckets = static_cast<std::uint32_t>(nBuckets());
0216   m_device.capacityPerBucket = static_cast<std::uint32_t>(capacityPerBucket());
0217
0218   try {
0219     allocateDeviceColumn(m_device.tanBeta, totalCapacity());
0220     allocateDeviceColumn(m_device.interceptY, totalCapacity());
0221
0222     allocateDeviceColumn(m_device.nHits, totalCapacity());
0223     allocateDeviceColumn(m_device.nLayers, totalCapacity());
0224     allocateDeviceColumn(m_device.layerMask, totalCapacity());
0225
0226     allocateDeviceColumn(m_device.xBin, totalCapacity());
0227     allocateDeviceColumn(m_device.yBin, totalCapacity());
0228
0229     allocateDeviceColumn(m_device.nMaxima, nBuckets());
0230     allocateDeviceColumn(m_device.nAssociatedHits, totalCapacity());
0231   } catch (...) {
0232     clearDevice();
0233     throw;
0234   }
0235
0236   m_onDevice = true;
0237 }
0238
0239 void CudaHoughMaximumBatch::moveToHost(cudaStream_t stream) {
0240   if (!m_onDevice) {
0241     return;
0242   }
0243
0244   copyColumnToHost(m_hostTanBeta, m_device.tanBeta, stream);
0245   copyColumnToHost(m_hostInterceptY, m_device.interceptY, stream);
0246   copyColumnToHost(m_hostHits, m_device.nHits, stream);
0247   copyColumnToHost(m_hostLayers, m_device.nLayers, stream);
0248   copyColumnToHost(m_hostLayerMask, m_device.layerMask, stream);
0249   copyColumnToHost(m_hostXBin, m_device.xBin, stream);
0250   copyColumnToHost(m_hostYBin, m_device.yBin, stream);
0251   copyColumnToHost(m_hostNMaxima, m_device.nMaxima, stream);
0252   copyColumnToHost(m_hostNAssociatedHits, m_device.nAssociatedHits, stream);
0253   ACTS_CUDA_CHECK(cudaStreamSynchronize(stream));
0254
0255   m_associationMetadataOnHost = true;
0256 }
0257
0258 void CudaHoughMaximumBatch::copyAssociationMetadataToHost(cudaStream_t stream) {
0259   if (!m_onDevice) {
0260     throw std::logic_error("CudaHoughMaximumBatch is not on the device");
0261   }
0262
0263   copyColumnToHost(m_hostNMaxima, m_device.nMaxima, stream);
0264   copyColumnToHost(m_hostNAssociatedHits, m_device.nAssociatedHits, stream);
0265   ACTS_CUDA_CHECK(cudaStreamSynchronize(stream));
0266
0267   m_associationMetadataOnHost = true;
0268 }
0269
0270 void CudaHoughMaximumBatch::prepareAssociationStorageHost() {
0271   clearAssociationStorage();
0272   m_hostAssociatedHitOffsets.assign(totalCapacity() + 1u, 0u);
0273   m_hostAssociatedHitIndices.clear();
0274
0275   std::uint64_t totalAssociatedHits = 0u;
0276   for (size_type bucket = 0u; bucket < nBuckets(); ++bucket) {
0277     const size_type maximaInBucket = nMaxima(bucket);
0278     for (size_type maximum = 0u; maximum < capacityPerBucket(); ++maximum) {
0279       const size_type slot = slotIndex(bucket, maximum);
0280       const std::uint32_t count = m_hostNAssociatedHits[slot];
0281       if (maximum >= maximaInBucket && count != 0u) {
0282         throw std::runtime_error(
0283             "Unoccupied maximum slot contains associated hits");
0284       }
0285       if (maximum < maximaInBucket) {
0286         totalAssociatedHits += count;
0287       }
0288       if (totalAssociatedHits > std::numeric_limits<std::uint32_t>::max()) {
0289         throw std::overflow_error(
0290             "Total associated hit count must fit into std::uint32_t and be "
0291             "allocatable");
0292       }
0293       m_hostAssociatedHitOffsets[slot + 1u] =
0294           static_cast<std::uint32_t>(totalAssociatedHits);
0295     }
0296   }
0297   m_hostAssociatedHitIndices.resize(
0298       static_cast<size_type>(totalAssociatedHits));
0299 }
0300
0301 void CudaHoughMaximumBatch::allocateAssociationStorage(cudaStream_t stream) {
0302   if (!m_onDevice) {
0303     throw std::logic_error("CudaHoughMaximumBatch is not on the device");
0304   }
0305
0306   if (!m_associationMetadataOnHost) {
0307     throw std::logic_error(
0308         "Association metadata must be copied before allocation");
0309   }
0310
0311   prepareAssociationStorageHost();
0312
0313   try {
0314     allocateDeviceColumn(m_device.associatedHitOffsets,
0315                          m_hostAssociatedHitOffsets.size());
0316     copyColumnToDevice(m_device.associatedHitOffsets,
0317                        m_hostAssociatedHitOffsets, stream);
0318     ACTS_CUDA_CHECK(cudaStreamSynchronize(stream));
0319
0320     allocateDeviceColumn(m_device.associatedHitIndices,
0321                          m_hostAssociatedHitIndices.size());
0322
0323     m_device.totalAssociatedHits =
0324         static_cast<std::uint32_t>(m_hostAssociatedHitIndices.size());
0325
0326     m_associationStorageAllocated = true;
0327     m_associatedHitIndicesOnHost = false;
0328   } catch (...) {
0329     clearAssociationStorage();
0330     m_hostAssociatedHitOffsets.clear();
0331     m_hostAssociatedHitIndices.clear();
0332     throw;
0333   }
0334 }
0335
0336 void CudaHoughMaximumBatch::copyAssociatedHitIndicesToHost(
0337     cudaStream_t stream) {
0338   if (!m_onDevice) {
0339     throw std::logic_error("CudaHoughMaximumBatch is not on the device");
0340   }
0341
0342   if (!m_associationStorageAllocated) {
0343     throw std::logic_error("Association storage has not been allocated");
0344   }
0345
0346   copyColumnToHost(m_hostAssociatedHitIndices, m_device.associatedHitIndices,
0347                    stream);
0348   ACTS_CUDA_CHECK(cudaStreamSynchronize(stream));
0349
0350   m_associatedHitIndicesOnHost = true;
0351 }
0352
0353 void CudaHoughMaximumBatch::clearAssociationStorage() noexcept {
0354   freeDeviceColumn(m_device.associatedHitOffsets);
0355   freeDeviceColumn(m_device.associatedHitIndices);
0356
0357   m_device.totalAssociatedHits = 0u;
0358   m_associationStorageAllocated = false;
0359 }
0360
0361 void CudaHoughMaximumBatch::clearDevice() noexcept {
0362   freeDeviceColumn(m_device.tanBeta);
0363   freeDeviceColumn(m_device.interceptY);
0364
0365   freeDeviceColumn(m_device.nHits);
0366   freeDeviceColumn(m_device.nLayers);
0367   freeDeviceColumn(m_device.layerMask);
0368
0369   freeDeviceColumn(m_device.xBin);
0370   freeDeviceColumn(m_device.yBin);
0371
0372   freeDeviceColumn(m_device.nMaxima);
0373   freeDeviceColumn(m_device.nAssociatedHits);
0374
0375   clearAssociationStorage();
0376
0377   m_device = {};
0378   m_onDevice = false;
0379 }
0380
0381 void CudaHoughMaximumBatch::checkBucket(size_type bucket) const {
0382   if (bucket >= nBuckets()) {
0383     throw std::out_of_range("CudaHoughMaximumBatch bucket index out of range");
0384   }
0385 }
0386
0387 void CudaHoughMaximumBatch::checkMaximum(size_type bucket,
0388                                          size_type maximum) const {
0389   if (maximum >= nMaxima(bucket)) {
0390     throw std::out_of_range("CudaHoughMaximumBatch maximum index out of range");
0391   }
0392 }
0393
0394 }  // namespace ActsExamples