diff --git a/.gitignore b/.gitignore index e52ca6e..f7d7ea1 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,9 @@ spack* *.gdml *.log +cmake-build-debug +clueLib/cmake-build-debug +data/output build build.*.log NMake @@ -56,3 +59,6 @@ CTestTestfile.cmake CPackSourceConfig.cmake CPackConfig.cmake cmake_install.cmake + +# third party libraries +alpaka/ diff --git a/.gitmodules b/.gitmodules index 3d3651f..f3fdd18 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +1,4 @@ [submodule "alpaka"] - path = alpaka - url = https://github.com/psychocoderHPC/alpaka2.git - branch = dev + path = alpaka + url = git@github.com:alpaka-group/alpaka3.git + branch = dev diff --git a/CMakeLists.txt b/CMakeLists.txt index afe30d8..271da8d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,9 +9,6 @@ project( # Location of the ALPAKA set(ALPAKA_DIR ${CMAKE_CURRENT_SOURCE_DIR}/alpaka CACHE PATH "Path to ALPAKA (by default: submodule)") -# Location of TBB -#set(TBB_DIR "/cvmfs/cms.cern.ch/slc7_amd64_gcc820/external/tbb/2020_U2-ghbfee/cmake/TBB/") - # Activate VERBOSE output set(CMAKE_VERBOSE_MAKEFILE OFF) @@ -20,10 +17,12 @@ include(CheckLanguage) check_language(CUDA) if(CMAKE_CUDA_COMPILER) enable_language(CUDA) - set(CMAKE_CXX_STANDARD 17) - set(CMAKE_CUDA_ARCHITECTURES "60;70;75") + set(CMAKE_CXX_STANDARD 20) + if(NOT CMAKE_CUDA_ARCHITECTURES) + set(CMAKE_CUDA_ARCHITECTURES "60;70;75;80") + endif() else() - message(STATUS "No CUDA compiler found. Still, you can run the C++ version!") + message(STATUS "No CUDA compiler found. Still, you can run the C++ version and alpaka!") endif() if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) @@ -33,56 +32,60 @@ endif() include(GNUInstallDirs) # Set up C++ Standard -set(CMAKE_CXX_STANDARD 17 CACHE STRING "") +set(CMAKE_CXX_STANDARD 20 CACHE STRING "") if(NOT CMAKE_CXX_STANDARD MATCHES "17|20") message(FATAL_ERROR "Unsupported C++ standard: ${CMAKE_CXX_STANDARD}") endif() -# Find Boost -set(Boost_DEBUG 1) -find_package(Boost REQUIRED) - -# Include Boost headers -include_directories(${Boost_INCLUDE_DIRS}) - -if(Boost_FOUND) - message(STATUS "Boost package found!") -endif() - -#find_package(TBB REQUIRED) -#if(TBB_FOUND) -# message(STATUS "TBB package found!") -#endif() +add_subdirectory("${CMAKE_CURRENT_LIST_DIR}/alpaka") add_subdirectory(clueLib) add_subdirectory(src) # SOME TESTS enable_testing() + +# Shared implementation (run test [cpu/gpu]) +function(run_test_common target label repeats suffix extra_args) + file(GLOB inputs "${CMAKE_CURRENT_SOURCE_DIR}/data/input/*.csv") + foreach(input_file IN LISTS inputs) + get_filename_component(input_dir "${input_file}" DIRECTORY) + get_filename_component(input_name "${input_file}" NAME_WE) + get_filename_component(input_ext "${input_file}" EXT) + + set(test_name "${input_name}_${label}_${suffix}") + set(output_file "${input_dir}/${test_name}${input_ext}") + + add_test( + NAME "${test_name}" + COMMAND $ + -i "${input_file}" + -O "${output_file}" + -d 7.0 -r10.0 -o 2 -e "${repeats}" -v + ${extra_args} + ) + endforeach() +endfunction() + function(run_test_cpu target label repeats) - FILE(GLOB inputs "${CMAKE_CURRENT_SOURCE_DIR}/data/input/*.csv") - foreach(input_file IN LISTS inputs) - add_test(NAME ${input_file}_${label}_CPU COMMAND ${target} -i ${input_file} -d 7.0 -r10.0 -o 2 -e ${repeats} -v) - endforeach() + run_test_common("${target}" "${label}" "${repeats}" "CPU" "") endfunction() -function(run_test_gpu target label repeats) - FILE(GLOB inputs "${CMAKE_CURRENT_SOURCE_DIR}/data/input/*.csv") - foreach(input_file IN LISTS inputs) - add_test(NAME ${input_file}_${label}_GPU COMMAND ${target} -i ${input_file} -d 7.0 -r10.0 -o 2 -e ${repeats} -v -u) - endforeach() +function(run_test_gpu target label repeats executor) + # Pass the executor-specific args through to the shared function + run_test_common("${target}" "${label}" "${repeats}" "${executor}" "-u;${executor}") endfunction() -if(CMAKE_CUDA_COMPILER) - run_test_cpu(./src/clue/main main 4) - run_test_cpu(./src/clue_cuda_alpaka/mainAlpakaCUDA Alpaka 4) - run_test_gpu(./src/clue/main CUDA 100) - run_test_gpu(./src/clue_cuda_alpaka/mainAlpakaCUDA AlpakaCUDA 100) -#else() - #run_test_cpu(./src/clue_tbb_alpaka/mainAlpakaTBB TBB 4) +if(alpaka_EXEC_GpuCuda) + run_test_gpu(mainAlpaka Alpaka 100 GpuCuda) endif() +if(alpaka_EXEC_CpuOmpBlocks) + run_test_gpu(mainAlpaka Alpaka 100 CpuOmpBlocks) +endif() + + #--- add version files --------------------------------------------------------- configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cmake/CLUEVersion.h ${CMAKE_CURRENT_BINARY_DIR}/CLUEVersion.h ) diff --git a/alpaka b/alpaka index a4142d3..3a40935 160000 --- a/alpaka +++ b/alpaka @@ -1 +1 @@ -Subproject commit a4142d3feb7686d803e1ec5f25d7b2278337f455 +Subproject commit 3a4093538f3c640d46c087e34adde503e1b35b00 diff --git a/clueLib/include/CLUEAlgo.h b/clueLib/include/CLUEAlgo.h index 5bdfd54..29f0b07 100644 --- a/clueLib/include/CLUEAlgo.h +++ b/clueLib/include/CLUEAlgo.h @@ -1,10 +1,5 @@ -#ifndef CLUEAlgo_h -#define CLUEAlgo_h - -// C/C++ headers -#include "Points.h" -#include "Tiles.h" - +#pragma once +// clang-format off #include #include #include @@ -12,6 +7,12 @@ #include #include +// C/C++ headers +#include "Points.h" +#include "Tiles.h" + +// clang-format on + // The type T is used to pass the number of bins in each dimension and the // allowed ranges spanned. Ancillary quantities, like the inverse of the bin // width should also be provided. Code will not compile if any such information @@ -139,7 +140,6 @@ class CLUEAlgo points_.weight.resize(n); std::copy(std::begin(weight), std::end(weight), std::begin(points_.weight)); - points_.p_x = points_.x.data(); points_.p_y = points_.y.data(); points_.p_layer = points_.layer.data(); @@ -394,7 +394,8 @@ void CLUEAlgo::makeClusters() // std::cout << "STANDALONE: before prepare" << std::endl; prepareDataStructures(allLayerTiles); - // std::cout << "STANDALONE: after prepare datastructures makeClusters" << std::endl; + // std::cout << "STANDALONE: after prepare datastructures makeClusters" << + // std::endl; auto finish = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed = finish - start; std::cout << "--- prepareDataStructures: " << elapsed.count() * 1000 << " ms\n"; @@ -411,7 +412,6 @@ void CLUEAlgo::makeClusters() elapsed = finish - start; std::cout << "--- calculateDistanceToHigher: " << elapsed.count() * 1000 << " ms\n"; - findAndAssignClusters(); // std::cout << "STANDALONE: end makeClusters" << std::endl; } @@ -600,5 +600,3 @@ inline float CLUEAlgo::distance(int i, int j) const float const dy = points_.p_y[i] - points_.p_y[j]; return std::sqrt(dx * dx + dy * dy); } - -#endif diff --git a/clueLib/include/CLUEAlgoAlpaka.h b/clueLib/include/CLUEAlgoAlpaka.h index 6908d33..a61adeb 100644 --- a/clueLib/include/CLUEAlgoAlpaka.h +++ b/clueLib/include/CLUEAlgoAlpaka.h @@ -1,8 +1,11 @@ #pragma once +// clang-format off +#include + #include "CLUEAlgo.h" #include "TilesAlpaka.h" -#include +// clang-format on #include #include @@ -16,39 +19,6 @@ }; \ ALPAKA_FN_ACC void operator()(ACC const& acc, Kernel##NAME dummy, ##__VA_ARGS__) const -// Main interface to select the proper accelerator **at compile time**. -namespace alpaka -{ - //! Alias for the default accelerator used by examples. From a list of - //! all accelerators the first one which is enabled is chosen. - //! AccCpuSerial is selected last. - template -#if defined(ALPAKA_ACC_GPU_CUDA_ENABLED) - using SelectedAcc = alpaka::AccGpuCudaRt; -#elif defined(ALPAKA_ACC_GPU_HIP_ENABLED) - using SelectedAcc = alpaka::AccGpuHipRt; -#elif defined(ALPAKA_ACC_CPU_B_OMP2_T_SEQ_ENABLED) - using SelectedAcc = alpaka::AccCpuOmp2Blocks; -#elif defined(ALPAKA_ACC_CPU_B_TBB_T_SEQ_ENABLED) - using SelectedAcc = alpaka::AccCpuTbbBlocks; -#elif defined(ALPAKA_ACC_CPU_B_SEQ_T_FIBERS_ENABLED) - using SelectedAcc = alpaka::AccCpuFibers; -#elif defined(ALPAKA_ACC_CPU_B_SEQ_T_OMP2_ENABLED) - using SelectedAcc = alpaka::AccCpuOmp2Threads; -#elif defined(ALPAKA_ACC_CPU_B_SEQ_T_THREADS_ENABLED) - using SelectedAcc = alpaka::AccCpuThreads; -#elif defined(ALPAKA_ACC_ANY_BT_OMP5_ENABLED) - using SelectedAcc = alpaka::AccOmp5; -#elif defined(ALPAKA_ACC_ANY_BT_OACC_ENABLED) - using SelectedAcc = alpaka::AccOacc; -#elif defined(ALPAKA_ACC_CPU_B_SEQ_T_SEQ_ENABLED) - using SelectedAcc = alpaka::AccCpuSerial; -#else - class SelectedAcc; -# warning "No supported backend selected." -#endif -} // namespace alpaka - // Maximum number of uniques seeds that could be handled. A higher number of // potential seed will trigger an exception. static int const maxNSeeds = 262144; @@ -67,24 +37,27 @@ static int const localStackSizePerSeed = 128; // allowed ranges spanned. Anchillary quantitied, like the inverse of the bin // width should also be provided. Code will not compile if any such information // is missing. -template +template class CLUEAlgoAlpaka : public CLUEAlgo { public: - using Dim = alpaka::Dim; - using Idx = alpaka::Idx; + static constexpr uint32_t dim = 1u; + using Idx = uint32_t; template - using BufAccT = alpaka::Buf; + using BufAccT = ALPAKA_TYPEOF( + alpaka::onHost::alloc(std::declval(), std::declval>())); template - using ViewHostT = alpaka::ViewPlainPtr; + using ViewHostT + = ALPAKA_TYPEOF(alpaka::onHost::alloc(std::declval(), std::declval>())); - using LayerTilesAcc = TilesAlpaka; + using LayerTilesAcc = TilesAlpaka; using BufLayerTiles = BufAccT; using BufVecArrSeeds = BufAccT>; using BufVecArrFollowers = BufAccT>; + // // Bring base-class public variables into the scope of this template derived // class @@ -134,62 +107,62 @@ class CLUEAlgoAlpaka : public CLUEAlgo uint8_t* isSeed; // LayerTiles and utility data structures - TilesAlpaka* hist_; + TilesAlpaka* hist_; GPUAlpaka::VecArray* seeds_; GPUAlpaka::VecArray* followers_; }; // Kernel and KernelTask definitions - DECLARE_TASKTYPE_AND_KERNEL(TAcc, ComputeHistogram, unsigned int const num_elements); - DECLARE_TASKTYPE_AND_KERNEL(TAcc, SortHistogram); - DECLARE_TASKTYPE_AND_KERNEL(TAcc, ComputeLocalDensity, float dc, unsigned int const num_elements); + DECLARE_TASKTYPE_AND_KERNEL(auto, ComputeHistogram, unsigned int const num_elements); + DECLARE_TASKTYPE_AND_KERNEL(auto, SortHistogram); + DECLARE_TASKTYPE_AND_KERNEL(auto, ComputeLocalDensity, float dc, unsigned int const num_elements); DECLARE_TASKTYPE_AND_KERNEL( - TAcc, + auto, ComputeDistanceToHigherNoDetId, float outlierDeltaFactor, float dc, unsigned int const num_elements); DECLARE_TASKTYPE_AND_KERNEL( - TAcc, + auto, ComputeDistanceToHigher, float outlierDeltaFactor, float dc, unsigned int const num_elements); DECLARE_TASKTYPE_AND_KERNEL( - TAcc, + auto, FindClusters, float outlierDeltaFactor, float dc, float rhoc, unsigned int const num_elements); DECLARE_TASKTYPE_AND_KERNEL( - TAcc, + auto, FindClustersKappa, float outlierDeltaFactor, float dc, float kappa, unsigned int const num_elements); DECLARE_TASKTYPE_AND_KERNEL( - TAcc, + auto, AssignClusters, unsigned int* numberOfClustersScalar, bool writeOutNumClusters = true); DeviceRawPointers ptrs_; }; - CLUEAlgoAlpaka(TQueue& queue, float dc, float kappa, float outlierDeltaFactor, bool verbose) - : CLUEAlgo(dc, kappa, outlierDeltaFactor, verbose) - , device_(alpaka::getDev(queue)) - , queue_(queue) - , host_(alpaka::getDevByIdx(alpaka::Platform{}, 0u)) - { - } - - CLUEAlgoAlpaka(float dc, float kappa, float outlierDeltaFactor, bool verbose, bool useAbsoluteSigma = false) + CLUEAlgoAlpaka( + TComputeDevice& computeDevice, + TQueue& queue, + THostDevice& hostDevice, + float dc, + float kappa, + float outlierDeltaFactor, + bool verbose, + bool useAbsoluteSigma = false) : CLUEAlgo(dc, kappa, outlierDeltaFactor, verbose, useAbsoluteSigma) - , device_(alpaka::getDevByIdx(alpaka::Platform{}, 0u)) - , queue_(device_) - , host_(alpaka::getDevByIdx(alpaka::Platform{}, 0u)) + , device_(computeDevice) + , queue_(queue) + , host_(hostDevice) { init_device(); } @@ -219,10 +192,10 @@ class CLUEAlgoAlpaka : public CLUEAlgo DeviceRunner device_runner_; private: - alpaka::Dev device_; + TComputeDevice device_; // choose between Blocking and NonBlocking TQueue queue_; - alpaka::DevCpu host_; + THostDevice host_; // Memory management variables PointsBuf device_bufs_; @@ -235,54 +208,54 @@ class CLUEAlgoAlpaka : public CLUEAlgo Idx const reserve = 1'000'000; // If Dim is not 1, fail compilation. This is assumed to be a // mono-dimensional problem - static_assert(Dim::value == 1u); - alpaka::Vec const extents(reserve); + static_assert(dim == 1u); + alpaka::Vec const extents(reserve); // INPUT VARIABLES // Allocate device memory - device_bufs_.x = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.y = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.layer = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.weight = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.sigmaNoise = std::make_optional(alpaka::allocBuf(device_, extents)); + device_bufs_.x = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.y = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.layer = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.weight = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.sigmaNoise = std::make_optional(alpaka::onHost::alloc(device_, extents)); // RESULT VARIABLES - device_bufs_.rho = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.delta = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.nearestHigher = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.clusterIndex = std::make_optional(alpaka::allocBuf(device_, extents)); - device_bufs_.isSeed = std::make_optional(alpaka::allocBuf(device_, extents)); + device_bufs_.rho = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.delta = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.nearestHigher = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.clusterIndex = std::make_optional(alpaka::onHost::alloc(device_, extents)); + device_bufs_.isSeed = std::make_optional(alpaka::onHost::alloc(device_, extents)); // INTERNAL VARIABLES - alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); - device_hist_ = std::make_optional(alpaka::allocBuf(device_, layerTilesExtents)); + alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); + device_hist_ = std::make_optional(alpaka::onHost::alloc(device_, layerTilesExtents)); - alpaka::Vec const seedsExtents(1u); + alpaka::Vec const seedsExtents(1u); device_seeds_ - = std::make_optional(alpaka::allocBuf, Idx>(device_, seedsExtents)); + = std::make_optional(alpaka::onHost::alloc>(device_, seedsExtents)); device_followers_ - = std::make_optional(alpaka::allocBuf, Idx>(device_, extents)); + = std::make_optional(alpaka::onHost::alloc>(device_, extents)); // Update RAW device pointers, grouped in a struct for convenience - device_runner_.ptrs_.x = alpaka::getPtrNative(device_bufs_.x.value()); - device_runner_.ptrs_.y = alpaka::getPtrNative(device_bufs_.y.value()); - device_runner_.ptrs_.layer = alpaka::getPtrNative(device_bufs_.layer.value()); - device_runner_.ptrs_.weight = alpaka::getPtrNative(device_bufs_.weight.value()); + device_runner_.ptrs_.x = alpaka::onHost::data(device_bufs_.x.value()); + device_runner_.ptrs_.y = alpaka::onHost::data(device_bufs_.y.value()); + device_runner_.ptrs_.layer = alpaka::onHost::data(device_bufs_.layer.value()); + device_runner_.ptrs_.weight = alpaka::onHost::data(device_bufs_.weight.value()); if(useAbsoluteSigma_) - device_runner_.ptrs_.sigmaNoise = alpaka::getPtrNative(device_bufs_.sigmaNoise.value()); + device_runner_.ptrs_.sigmaNoise = alpaka::onHost::data(device_bufs_.sigmaNoise.value()); // RESULT VARIABLES - device_runner_.ptrs_.rho = alpaka::getPtrNative(device_bufs_.rho.value()); - device_runner_.ptrs_.delta = alpaka::getPtrNative(device_bufs_.delta.value()); - device_runner_.ptrs_.nearestHigher = alpaka::getPtrNative(device_bufs_.nearestHigher.value()); - device_runner_.ptrs_.clusterIndex = alpaka::getPtrNative(device_bufs_.clusterIndex.value()); - device_runner_.ptrs_.isSeed = alpaka::getPtrNative(device_bufs_.isSeed.value()); + device_runner_.ptrs_.rho = alpaka::onHost::data(device_bufs_.rho.value()); + device_runner_.ptrs_.delta = alpaka::onHost::data(device_bufs_.delta.value()); + device_runner_.ptrs_.nearestHigher = alpaka::onHost::data(device_bufs_.nearestHigher.value()); + device_runner_.ptrs_.clusterIndex = alpaka::onHost::data(device_bufs_.clusterIndex.value()); + device_runner_.ptrs_.isSeed = alpaka::onHost::data(device_bufs_.isSeed.value()); // UPDATE RAW POINTERS FOR INTERNATL DATA STRUCTURES - device_runner_.ptrs_.hist_ = alpaka::getPtrNative(device_hist_.value()); - device_runner_.ptrs_.seeds_ = alpaka::getPtrNative(device_seeds_.value()); - device_runner_.ptrs_.followers_ = alpaka::getPtrNative(device_followers_.value()); + device_runner_.ptrs_.hist_ = alpaka::onHost::data(device_hist_.value()); + device_runner_.ptrs_.seeds_ = alpaka::onHost::data(device_seeds_.value()); + device_runner_.ptrs_.followers_ = alpaka::onHost::data(device_followers_.value()); } void free_device() @@ -296,136 +269,148 @@ class CLUEAlgoAlpaka : public CLUEAlgo // This means the view will, possibly, become invalid if, in the meantime, // the vector re-allocated its underlying storage. template - auto getViewHost(TT& t) -> ViewHostT + auto getViewHost(std::vector& t) { - using type = typename TT::value_type; - using Dim1 = alpaka::DimInt<1ul>; - alpaka::Vec vectorSize(static_cast(t.size())); - ViewHostT tempHostView(t.data(), host_, vectorSize); - return tempHostView; + using type = typename std::vector::value_type; + using ExtentType = alpaka::Vec; + ExtentType vectorSize(static_cast(t.size())); + + auto deleter = [](type* ptr) {}; + auto pitches = ExtentType{sizeof(type)}; + /* + alpaka::onHost::data() + alpaka::makeMdSpan(t.data(),) + auto data = std::make_shared>( host_, t.data(), vectorSize, pitches, std::move(deleter));*/ + return alpaka::View(alpaka::api::host, t.data(), vectorSize, pitches); } template - auto getViewHost(TT* t, int size) -> ViewHostT + auto getViewHost(TT* t, int size) { - using Dim1 = alpaka::DimInt<1ul>; - alpaka::Vec vectorSize(static_cast(size)); - ViewHostT tempHostView(t, host_, vectorSize); - return tempHostView; + using ExtentType = alpaka::Vec; + ExtentType vectorSize(static_cast(size)); + + auto deleter = [](TT* ptr) {}; + auto pitches = ExtentType{sizeof(TT)}; + return alpaka::View(alpaka::api::host, t, vectorSize, pitches); } void copy_todevice() { // input variables - using Dim1 = alpaka::DimInt<1ul>; - alpaka::Vec const extentToTransfer(static_cast(points_.n)); - alpaka::memcpy(queue_, device_bufs_.x.value(), getViewHost(points_.p_x, points_.n), extentToTransfer); - alpaka::memcpy(queue_, device_bufs_.y.value(), getViewHost(points_.p_y, points_.n), extentToTransfer); - alpaka::memcpy(queue_, device_bufs_.layer.value(), getViewHost(points_.p_layer, points_.n), extentToTransfer); - alpaka::memcpy( + alpaka::Vec const extentToTransfer(static_cast(points_.n)); + alpaka::onHost::memcpy(queue_, device_bufs_.x.value(), getViewHost(points_.p_x, points_.n), extentToTransfer); + alpaka::onHost::memcpy(queue_, device_bufs_.y.value(), getViewHost(points_.p_y, points_.n), extentToTransfer); + alpaka::onHost::memcpy( + queue_, + device_bufs_.layer.value(), + getViewHost(points_.p_layer, points_.n), + extentToTransfer); + alpaka::onHost::memcpy( queue_, device_bufs_.weight.value(), getViewHost(points_.p_weight, points_.n), extentToTransfer); if(useAbsoluteSigma_) - alpaka::memcpy( + alpaka::onHost::memcpy( queue_, device_bufs_.sigmaNoise.value(), getViewHost(points_.p_sigmaNoise, points_.n), extentToTransfer); - alpaka::wait(queue_); + alpaka::onHost::wait(queue_); } void clear_internal_buffers() { // result variables - using Dim1 = alpaka::DimInt<1ul>; - alpaka::Vec extents(static_cast(points_.n)); - alpaka::memset(queue_, device_bufs_.rho.value(), 0x0, extents); - alpaka::memset(queue_, device_bufs_.delta.value(), 0x0, extents); - alpaka::memset(queue_, device_bufs_.nearestHigher.value(), 0x0, extents); - alpaka::memset(queue_, device_bufs_.clusterIndex.value(), 0x0, extents); - alpaka::memset(queue_, device_bufs_.isSeed.value(), 0x0, extents); + alpaka::Vec extents(static_cast(points_.n)); + alpaka::onHost::memset(queue_, device_bufs_.rho.value(), 0x0, extents); + alpaka::onHost::memset(queue_, device_bufs_.delta.value(), 0x0, extents); + alpaka::onHost::memset(queue_, device_bufs_.nearestHigher.value(), 0x0, extents); + alpaka::onHost::memset(queue_, device_bufs_.clusterIndex.value(), 0x0, extents); + alpaka::onHost::memset(queue_, device_bufs_.isSeed.value(), 0x0, extents); // algorithm internal variables // INTERNAL VARIABLES - alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); - alpaka::memset(queue_, device_hist_.value(), 0x0, layerTilesExtents); + alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); + alpaka::onHost::memset(queue_, device_hist_.value(), 0x0, layerTilesExtents); - alpaka::Vec const seedsExtents(1u); - alpaka::memset(queue_, device_seeds_.value(), 0x0, seedsExtents); + alpaka::Vec const seedsExtents(1u); + alpaka::onHost::memset(queue_, device_seeds_.value(), 0x0, seedsExtents); - alpaka::memset(queue_, device_followers_.value(), 0x0, extents); - alpaka::wait(queue_); + alpaka::onHost::memset(queue_, device_followers_.value(), 0x0, extents); + alpaka::onHost::wait(queue_); } void copy_tohost() { // result variables - using Dim1 = alpaka::DimInt<1ul>; - alpaka::Vec extents(static_cast(points_.n)); + alpaka::Vec extents(static_cast(points_.n)); auto clusterHV = getViewHost(points_.clusterIndex); - alpaka::memcpy(queue_, clusterHV, device_bufs_.clusterIndex.value(), extents); + alpaka::onHost::memcpy(queue_, clusterHV, device_bufs_.clusterIndex.value(), extents); if(verbose_) { // other variables, copy only when verbose_==True auto rhoHV = getViewHost(points_.rho); - alpaka::memcpy(queue_, rhoHV, device_bufs_.rho.value(), extents); + alpaka::onHost::memcpy(queue_, rhoHV, device_bufs_.rho.value(), extents); auto deltaHV = getViewHost(points_.delta); - alpaka::memcpy(queue_, deltaHV, device_bufs_.delta.value(), extents); + alpaka::onHost::memcpy(queue_, deltaHV, device_bufs_.delta.value(), extents); auto nearestHV = getViewHost(points_.nearestHigher); - alpaka::memcpy(queue_, nearestHV, device_bufs_.nearestHigher.value(), extents); + alpaka::onHost::memcpy(queue_, nearestHV, device_bufs_.nearestHigher.value(), extents); auto isSeedHV = getViewHost(points_.isSeed); - alpaka::memcpy(queue_, isSeedHV, device_bufs_.isSeed.value(), extents); + alpaka::onHost::memcpy(queue_, isSeedHV, device_bufs_.isSeed.value(), extents); } - alpaka::wait(queue_); + alpaka::onHost::wait(queue_); } }; -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelComputeHistogram dummy, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner::KernelComputeHistogram + dummy, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); - if(i < numberOfPoints) + for(auto [i] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{numberOfPoints})) { // push index of points into tiles - ptrs_.hist_[ptrs_.layer[i]].fill(ptrs_.x[i], ptrs_.y[i], i); + ptrs_.hist_[ptrs_.layer[i]].fill(acc, ptrs_.x[i], ptrs_.y[i], i); } } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelSortHistogram dummy) const -> void +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner::KernelSortHistogram + dummy) const -> void { - /* - const Idx layer(alpaka::getIdx(acc)[0u]); - ptrs_.hist_[layer].sort_unsafe(); - */ - Idx const i(alpaka::getIdx(acc)[0u]); - if(i < NLAYERS * T::nTiles) + for(auto [i] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{NLAYERS * T::nTiles})) { int layer = i / T::nTiles; int bin = i - layer * T::nTiles; - ptrs_.hist_[layer].sort_unsafe(bin); + ptrs_.hist_[layer].sort_unsafe(acc, bin); } } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelComputeLocalDensity dummy, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner::KernelComputeLocalDensity + dummy, float dc, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); - if(i < numberOfPoints) + for(auto [i] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{numberOfPoints})) { float rhoi{0.}; int layeri = ptrs_.layer[i]; @@ -465,18 +450,20 @@ ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::opera } } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelComputeDistanceToHigherNoDetId dummy, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeDistanceToHigherNoDetId dummy, float outlierDeltaFactor, float dc, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); float dm = outlierDeltaFactor * dc; - if(i < numberOfPoints) + for(auto [i] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{numberOfPoints})) { int layeri = ptrs_.layer[i]; @@ -520,25 +507,19 @@ ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::opera bool foundHigher = (ptrs_.rho[j] > rhoi); // in the rare case where rho is the same, use indices foundHigher = foundHigher || ((ptrs_.rho[j] == rhoi) && (j > i)); - if(foundHigher && dist_ij < deltai) - { - rho_max = ptrs_.rho[j]; - deltai = dist_ij; - nearestHigheri = j; - } - else if(foundHigher && dist_ij == deltai && ptrs_.rho[j] > rho_max) - { - rho_max = ptrs_.rho[j]; - deltai = dist_ij; - nearestHigheri = j; - } - else if(foundHigher && dist_ij == deltai && ptrs_.rho[j] == rho_max && j > i) + + // combine the three conditions to avoid branching into different code segments to speedup the + // memory write a little bit + bool condition + = foundHigher && dist_ij < deltai + || dist_ij == deltai && ((ptrs_.rho[j] > rho_max) || (ptrs_.rho[j] == rho_max && j > i)); + if(condition) { rho_max = ptrs_.rho[j]; deltai = dist_ij; nearestHigheri = j; } - } // end of interate inside this bin + } } } // end of loop over bins in search box ptrs_.delta[i] = std::sqrt(deltai); @@ -546,18 +527,20 @@ ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::opera } } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelComputeDistanceToHigher dummy, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeDistanceToHigher dummy, float outlierDeltaFactor, float dc, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); float dm = outlierDeltaFactor * dc; - if(i < numberOfPoints) + for(auto [i] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{numberOfPoints})) { int layeri = ptrs_.layer[i]; @@ -601,206 +584,336 @@ ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::opera bool foundHigher = (ptrs_.rho[j] > rhoi); // in the rare case where rho is the same, use detid foundHigher = foundHigher || ((ptrs_.rho[j] == rhoi) && (ptrs_.detid[j] > ptrs_.detid[i])); - if(foundHigher && dist_ij < deltai) - { - rho_max = ptrs_.rho[j]; - deltai = dist_ij; - nearestHigheri = j; - } - else if(foundHigher && dist_ij == deltai && ptrs_.rho[j] > rho_max) - { - rho_max = ptrs_.rho[j]; - deltai = dist_ij; - nearestHigheri = j; - } - else if( - foundHigher && dist_ij == deltai && ptrs_.rho[j] == rho_max && ptrs_.detid[j] > ptrs_.detid[i]) + + // combine the three conditions to avoid branching into different code segments to speedup the + // memory write a little bit + bool condition = foundHigher && dist_ij < deltai + || dist_ij == deltai + && ((ptrs_.rho[j] > rho_max) + || (ptrs_.rho[j] == rho_max && ptrs_.detid[j] > ptrs_.detid[i])); + if(condition) { rho_max = ptrs_.rho[j]; deltai = dist_ij; nearestHigheri = j; } - } // end of interate inside this bin - } + } + } // end of interate inside this bin } // end of loop over bins in search box ptrs_.delta[i] = std::sqrt(deltai); ptrs_.nearestHigher[i] = nearestHigheri; } } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, KernelFindClustersKappa dummy, float outlierDeltaFactor, float dc, float kappa, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); - - if(i < numberOfPoints) - { - // initialize clusterIndex - ptrs_.clusterIndex[i] = -1; - // determine seed or outlier - float deltai = ptrs_.delta[i]; - float rhoi = ptrs_.rho[i]; - float rhoc = ptrs_.sigmaNoise[i] * kappa; - bool isSeed = (deltai > dc) && (rhoi >= rhoc); - bool isOutlier = (deltai > outlierDeltaFactor * dc) && (rhoi < rhoc); - - if(isSeed) + /* This looks boring and is boring. + * Normally the MdSPan object would be passed as argument but since the original code only stores pointers I + * decided to create the MdSpan on the fly. The alignment is set to a safe value and could also be `sizeof(type) * + * simdWidth` because alpaka is always aligning allocated memory to simdWidth. MdSpan is required for the + * concurrent for each to generate SIMD code and/or at least code instruction level parallel code. + */ + auto clusterIndexSpan = makeMdSpan( + ptrs_.clusterIndex, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(int)}, + alpaka::Alignment{}); + + auto deltaSpan = makeMdSpan( + ptrs_.delta, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(float)}, + alpaka::Alignment{}); + + auto rohSpan = makeMdSpan( + ptrs_.rho, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(float)}, + alpaka::Alignment{}); + + auto sigmaNoiseSpan = makeMdSpan( + ptrs_.sigmaNoise, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(float)}, + alpaka::Alignment{}); + + auto simdGrid = alpaka::onAcc::SimdAlgo{alpaka::onAcc::worker::threadsInGrid}; + simdGrid.concurrent( + acc, + [&](auto const&, auto&& simdClusterIdx, auto&& simdDelta, auto&& simdRoh, auto&& simdSigmaNoise) constexpr { - // set isSeed as 1 - ptrs_.isSeed[i] = 1; - ptrs_.seeds_[0].push_back(acc, i); // head of device_seeds_ - } - else - { - if(!isOutlier) + // initialize clusterIndex + using T_SimdType = ALPAKA_TYPEOF(simdClusterIdx.load()); + simdClusterIdx = T_SimdType::fill(-1); + + // determine seed or outlier + alpaka::concepts::Simd auto deltai = simdDelta.load(); + alpaka::concepts::Simd auto rhoi = simdRoh.load(); + alpaka::concepts::Simd auto rhoc = simdSigmaNoise.load() * kappa; + alpaka::concepts::Simd auto isSeed = (deltai > dc) && (rhoi >= rhoc); + alpaka::concepts::Simd auto isOutlier = (deltai > outlierDeltaFactor * dc) && (rhoi < rhoc); + + auto idxOffset = simdDelta.getIdx().x(); + for(int sIdx = 0; sIdx < alpaka::getDim(isSeed); ++sIdx) { - assert(ptrs_.nearestHigher[i] < numberOfPoints); - // register as follower at its nearest higher - ptrs_.followers_[ptrs_.nearestHigher[i]].push_back(acc, i); + int i = idxOffset + sIdx; + if(isSeed[sIdx]) + { + // set isSeed as 1 + ptrs_.isSeed[i] = 1; + // head of device_seeds_ + ptrs_.seeds_[0].push_back(acc, i); + } + else + { + if(!isOutlier[sIdx]) + { + assert(ptrs_.nearestHigher[i] < numberOfPoints); + // register as follower at its nearest higher + ptrs_.followers_[ptrs_.nearestHigher[i]].push_back(acc, i); + } + } } - } - } + }, + clusterIndexSpan, + deltaSpan, + rohSpan, + sigmaNoiseSpan); } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, KernelFindClusters dummy, float outlierDeltaFactor, float dc, float rhoc, unsigned int const numberOfPoints) const -> void { - Idx const i(alpaka::getIdx(acc)[0u]); - - if(i < numberOfPoints) - { - // initialize clusterIndex - ptrs_.clusterIndex[i] = -1; - // determine seed or outlier - float deltai = ptrs_.delta[i]; - float rhoi = ptrs_.rho[i]; - bool isSeed = (deltai > dc) && (rhoi >= rhoc); - bool isOutlier = (deltai > outlierDeltaFactor * dc) && (rhoi < rhoc); - - if(isSeed) - { - // set isSeed as 1 - ptrs_.isSeed[i] = 1; - ptrs_.seeds_[0].push_back(acc, i); // head of device_seeds_ - } - else + // see the kernel above + auto clusterIndexSpan = makeMdSpan( + ptrs_.clusterIndex, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(int)}, + alpaka::Alignment{}); + + auto deltaSpan = makeMdSpan( + ptrs_.delta, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(float)}, + alpaka::Alignment{}); + + auto rohSpan = makeMdSpan( + ptrs_.rho, + alpaka::Vec{numberOfPoints}, + alpaka::Vec{sizeof(float)}, + alpaka::Alignment{}); + + auto simdGrid = alpaka::onAcc::SimdAlgo{alpaka::onAcc::worker::threadsInGrid}; + simdGrid.concurrent( + acc, + [&](auto const&, auto&& simdClusterIdx, auto&& simdDelta, auto&& simdRoh) constexpr { - if(!isOutlier) + // initialize clusterIndex + using T_SimdType = ALPAKA_TYPEOF(simdClusterIdx.load()); + simdClusterIdx = T_SimdType::fill(-1); + + // determine seed or outlier + alpaka::concepts::Simd auto deltai = simdDelta.load(); + alpaka::concepts::Simd auto rhoi = simdRoh.load(); + alpaka::concepts::Simd auto isSeed = (deltai > dc) && (rhoi >= rhoc); + alpaka::concepts::Simd auto isOutlier = (deltai > outlierDeltaFactor * dc) && (rhoi < rhoc); + + auto idxOffset = simdDelta.getIdx().x(); + for(int sIdx = 0; sIdx < alpaka::getDim(isSeed); ++sIdx) { - assert(ptrs_.nearestHigher[i] < numberOfPoints); - // register as follower at its nearest higher - ptrs_.followers_[ptrs_.nearestHigher[i]].push_back(acc, i); + int i = idxOffset + sIdx; + if(isSeed[sIdx]) + { + // set isSeed as 1 + ptrs_.isSeed[i] = 1; + // head of device_seeds_ + ptrs_.seeds_[0].push_back(acc, i); + } + else + { + if(!isOutlier[sIdx]) + { + assert(ptrs_.nearestHigher[i] < numberOfPoints); + // register as follower at its nearest higher + ptrs_.followers_[ptrs_.nearestHigher[i]].push_back(acc, i); + } + } } - } - } + }, + clusterIndexSpan, + deltaSpan, + rohSpan); } -template -ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner::operator()( - TAcc const& acc, - CLUEAlgoAlpaka::DeviceRunner::KernelAssignClusters dummy, +template +ALPAKA_FN_ACC auto CLUEAlgoAlpaka::DeviceRunner:: +operator()( + auto const& acc, + CLUEAlgoAlpaka::DeviceRunner::KernelAssignClusters + dummy, unsigned int* numberOfClustersScalar, bool writeOutNumClusters) const -> void { - Idx const idxCls(alpaka::getIdx(acc)[0u]); - if(idxCls == 0 && writeOutNumClusters) + if(writeOutNumClusters) { - *numberOfClustersScalar = ptrs_.seeds_[0].size(); + for(auto [idxCls] : alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInGrid, alpaka::IdxRange{1u})) + { + *numberOfClustersScalar = ptrs_.seeds_[0].size(); + } } - if(idxCls < (unsigned int) ptrs_.seeds_[0].size()) + +#if ALPAKA_LANG_CUDA && __CUDA_ARCH__ + /* alpaka has currently not implemented warps. Never the less it is possible to define any kind of group therefore + * we hard code a warp size of 32 if CLUE_USE_CUDA_WARP is defined and use later warps in a grid to iterate over + * the seed noods and all therads within a warp to iterate of the followers of the seed node. + * Currently the start parameters are not adjusted for warps and we use always the same frame extent. + * Native support for warps is coming soon! */ +# define CLUE_USE_CUDA_WARP 1 +#endif + + // iterate with thread blocks over the seed particles + for(auto [idxCls] : alpaka::onAcc::makeIdxMap( + acc, +#if CLUE_USE_CUDA_WARP + alpaka::onAcc::WorkerGroup{ + (acc[alpaka::layer::block].idx() * acc[alpaka::layer::thread].count() + + acc[alpaka::layer::thread].idx()) + / 32u, + (acc[alpaka::layer::thread].count() * acc[alpaka::layer::block].count()) / 32u}, +#else + alpaka::onAcc::worker::blocksInGrid, +#endif + + alpaka::IdxRange{(unsigned int) ptrs_.seeds_[0].size()})) { int localStack[localStackSizePerSeed] = {-1}; int localStackSize = 0; // assign cluster to seed[idxCls] int idxThisSeed = ptrs_.seeds_[0][idxCls]; - ptrs_.clusterIndex[idxThisSeed] = idxCls; - // push_back idThisSeed to localStack - assert(localStackSize < localStackSizePerSeed); - localStack[localStackSize] = idxThisSeed; - localStackSize++; - - // process all elements in localStack - while(localStackSize > 0) + for(auto [idxCls] : + alpaka::onAcc::makeIdxMap(acc, alpaka::onAcc::worker::threadsInBlock, alpaka::IdxRange{1u})) + { + ptrs_.clusterIndex[idxThisSeed] = idxCls; + } + // the first level of the hierarchy will be processed by all threads in a block + for(auto [stackIdx] : alpaka::onAcc::makeIdxMap( + acc, +#if CLUE_USE_CUDA_WARP + alpaka::onAcc::WorkerGroup{ + alpaka::Vec{acc[alpaka::layer::thread].idx() % 32u}, + alpaka::Vec{32u}}, +#else + alpaka::onAcc::worker::threadsInBlock, +#endif + alpaka::IdxRange{alpaka::Vec{ptrs_.followers_[idxThisSeed].size()}})) { - // get last element of localStack - assert(localStackSize - 1 < localStackSizePerSeed); - int idxEndOflocalStack = localStack[localStackSize - 1]; - - int temp_clusterIndex = ptrs_.clusterIndex[idxEndOflocalStack]; - // pop_back last element of localStack - assert(localStackSize - 1 < localStackSizePerSeed); - localStack[localStackSize - 1] = -1; - localStackSize--; - - // loop over followers of last element of localStack - for(int j : ptrs_.followers_[idxEndOflocalStack]) + int rootIdx = ptrs_.followers_[idxThisSeed][stackIdx]; + ptrs_.clusterIndex[rootIdx] = idxCls; + int currentRootIdx = rootIdx; + + if(ptrs_.followers_[currentRootIdx].size() > 0) { - // pass id to follower - ptrs_.clusterIndex[j] = temp_clusterIndex; - // push_back follower to localStack - assert(localStackSize < localStackSizePerSeed); - localStack[localStackSize] = j; - localStackSize++; + // process all elements in localStack + do + { + // during the first visit of this loop we do not need from the stack + if(currentRootIdx != rootIdx) + { + // get last element of localStack + assert(localStackSize - 1 < localStackSizePerSeed); + currentRootIdx = localStack[localStackSize - 1]; + + // pop_back last element of localStack + localStack[localStackSize - 1] = -1; + localStackSize--; + } + + // loop over followers of last element of localStack + for(int j : ptrs_.followers_[currentRootIdx]) + { + // pass id to follower + ptrs_.clusterIndex[j] = idxCls; + // push only to the stack of we have followers to avoid useless memory operations and reduce + // the number of elements on te stack + if(ptrs_.followers_[currentRootIdx].size() > 0) + { + // push_back follower to localStack + assert(localStackSize < localStackSizePerSeed); + localStack[localStackSize] = j; + localStackSize++; + } + } + // reset the current root to load next index from the stack + currentRootIdx = -1; + } while(localStackSize > 0); } } } } -template -void CLUEAlgoAlpaka::makeClusters() +inline auto fString(std::string s) +{ + s.resize(std::max(s.size(), size_t{30}), ' '); + return s; +} + +template +void CLUEAlgoAlpaka::makeClusters() { copy_todevice(); clear_internal_buffers(); // Dimension the grid for submission - alpaka::Vec const threadsPerBlock(1024u); - alpaka::Vec const blocksPerGrid(static_cast(ceil(points_.n / (float) threadsPerBlock[0]))); - alpaka::Vec const elementsPerThread(1u); - using WorkDiv = alpaka::WorkDivMembers; - auto const manualWorkDiv = WorkDiv{blocksPerGrid, threadsPerBlock, elementsPerThread}; + alpaka::Vec const threadsPerBlock(1024u); + alpaka::Vec const blocksPerGrid(alpaka::divExZero(static_cast(points_.n), threadsPerBlock[0])); + + auto const manualWorkDiv = alpaka::onHost::FrameSpec{blocksPerGrid, threadsPerBlock}; // Create the kernel execution tasks. - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeHistogram taskComputeHistogram; - auto const kernelComputeHistogram = (alpaka::createTaskKernel( - manualWorkDiv, - device_runner_, - taskComputeHistogram, - static_cast(points_.n))); + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeHistogram taskComputeHistogram; + auto const kernelComputeHistogram + = alpaka::KernelBundle(device_runner_, taskComputeHistogram, static_cast(points_.n)); - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeLocalDensity taskComputeLocalDensity; - auto const kernelComputeLocalDensity = (alpaka::createTaskKernel( - manualWorkDiv, - device_runner_, - taskComputeLocalDensity, - dc_, - static_cast(points_.n))); + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeLocalDensity taskComputeLocalDensity; + auto const kernelComputeLocalDensity + = alpaka::KernelBundle(device_runner_, taskComputeLocalDensity, dc_, static_cast(points_.n)); - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeDistanceToHigherNoDetId - taskComputeDistanceToHigherNoDetId; - auto const kernelComputeDistanceToHigherNoDetId = (alpaka::createTaskKernel( - manualWorkDiv, + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeDistanceToHigherNoDetId taskComputeDistanceToHigherNoDetId; + auto const kernelComputeDistanceToHigherNoDetId = (alpaka::KernelBundle( device_runner_, taskComputeDistanceToHigherNoDetId, outlierDeltaFactor_, dc_, static_cast(points_.n))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelFindClusters taskFindClusters; - auto const kernelFindClusters = (alpaka::createTaskKernel( - manualWorkDiv, + // use int as data type since we handle indecision and float value in the kernels + uint32_t elementsPerFrameItem = alpaka::getNumElemPerThread(queue_); + alpaka::Vec const blocksPerGridSimd( + alpaka::divExZero(static_cast(points_.n), (threadsPerBlock[0] * elementsPerFrameItem))); + auto const manualWorkDivSimd = alpaka::onHost::FrameSpec{blocksPerGridSimd, threadsPerBlock}; + + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelFindClusters taskFindClusters; + auto const kernelFindClusters = (alpaka::KernelBundle( device_runner_, taskFindClusters, outlierDeltaFactor_, @@ -808,9 +921,9 @@ void CLUEAlgoAlpaka::makeClusters() rhoc_, static_cast(points_.n))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelFindClustersKappa taskFindClustersKappa; - auto const kernelFindClustersKappa = (alpaka::createTaskKernel( - manualWorkDiv, + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelFindClustersKappa taskFindClustersKappa; + auto const kernelFindClustersKappa = (alpaka::KernelBundle( device_runner_, taskFindClustersKappa, outlierDeltaFactor_, @@ -818,67 +931,76 @@ void CLUEAlgoAlpaka::makeClusters() kappa_, static_cast(points_.n))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelAssignClusters taskAssignClusters; - auto const kernelAssignClusters - = (alpaka::createTaskKernel(manualWorkDiv, device_runner_, taskAssignClusters, nullptr, false)); + // Dimension the grid for submission + alpaka::Vec const threadsPerBlockX(64u); + + // This value is too large and should be substituted with something reasonable + alpaka::Vec const blocksPerGridX(static_cast(points_.n)); + + auto const manualWorkDivX = alpaka::onHost::FrameSpec{blocksPerGridX, threadsPerBlockX}; + + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelAssignClusters taskAssignClusters; + auto const kernelAssignClusters = (alpaka::KernelBundle(device_runner_, taskAssignClusters, nullptr, false)); // Enqueue the kernel execution task auto start = std::chrono::high_resolution_clock::now(); auto finish = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed; + alpaka::onHost::wait(queue_); start = std::chrono::high_resolution_clock::now(); - alpaka::enqueue(queue_, kernelComputeHistogram); - alpaka::wait(queue_); // wait in case we are using an asynchronous queue to + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeHistogram); + alpaka::onHost::wait(queue_); // wait in case we are using an asynchronous queue to // time actual kernel runtime finish = std::chrono::high_resolution_clock::now(); elapsed = finish - start; - std::cout << "--- computeHistogram: " << elapsed.count() * 1000 << "ms\n"; + std::cout << fString("--- computeHistogram:") << elapsed.count() * 1000 << "ms\n"; start = std::chrono::high_resolution_clock::now(); - alpaka::enqueue(queue_, kernelComputeLocalDensity); - alpaka::wait(queue_); // wait in case we are using an asynchronous queue to + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeLocalDensity); + alpaka::onHost::wait(queue_); // wait in case we are using an asynchronous queue to // time actual kernel runtime finish = std::chrono::high_resolution_clock::now(); elapsed = finish - start; - std::cout << "--- computeLocalDensity: " << elapsed.count() * 1000 << "ms\n"; + std::cout << fString("--- computeLocalDensity:") << elapsed.count() * 1000 << "ms\n"; start = std::chrono::high_resolution_clock::now(); - alpaka::enqueue(queue_, kernelComputeDistanceToHigherNoDetId); - alpaka::wait(queue_); // wait in case we are using an asynchronous queue to + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeDistanceToHigherNoDetId); + alpaka::onHost::wait(queue_); // wait in case we are using an asynchronous queue to // time actual kernel runtime finish = std::chrono::high_resolution_clock::now(); elapsed = finish - start; - std::cout << "--- computeDistanceToHigher: " << elapsed.count() * 1000 << "ms\n"; + std::cout << fString("--- computeDistanceToHigher:") << elapsed.count() * 1000 << "ms\n"; start = std::chrono::high_resolution_clock::now(); if(useAbsoluteSigma_) { - alpaka::enqueue(queue_, kernelFindClustersKappa); + queue_.enqueue(TExecutor{}, manualWorkDivSimd, kernelFindClustersKappa); } else { - alpaka::enqueue(queue_, kernelFindClusters); + queue_.enqueue(TExecutor{}, manualWorkDivSimd, kernelFindClusters); } - alpaka::wait(queue_); // wait in case we are using an asynchronous queue to + alpaka::onHost::wait(queue_); // wait in case we are using an asynchronous queue to // time actual kernel runtime finish = std::chrono::high_resolution_clock::now(); elapsed = finish - start; - std::cout << "--- findClusters: " << elapsed.count() * 1000 << "ms\n"; + std::cout << fString("--- findClusters:") << elapsed.count() * 1000 << "ms\n"; start = std::chrono::high_resolution_clock::now(); - alpaka::enqueue(queue_, kernelAssignClusters); - alpaka::wait(queue_); // wait in case we are using an asynchronous queue to + queue_.enqueue(TExecutor{}, manualWorkDivX, kernelAssignClusters); + alpaka::onHost::wait(queue_); // wait in case we are using an asynchronous queue to // time actual kernel runtime finish = std::chrono::high_resolution_clock::now(); elapsed = finish - start; - std::cout << "--- assignClusters: " << elapsed.count() * 1000 << "ms\n"; + std::cout << fString("--- assignClusters:") << elapsed.count() * 1000 << "ms\n"; copy_tohost(); } -template -void CLUEAlgoAlpaka::makeClustersCMSSW( +template +void CLUEAlgoAlpaka::makeClustersCMSSW( unsigned int const points, float const* x, float const* y, @@ -893,27 +1015,29 @@ void CLUEAlgoAlpaka::makeClustersCMSSW( uint8_t* isSeed, unsigned int* numberOfClustersScalar) { - // std::cout << "makeClustersCMSSW received " << points << " RecHits" << std::endl; + // std::cout << "makeClustersCMSSW received " << points << " RecHits" << + // std::endl; // INTERNAL TILES VARIABLES Idx const reserve = 1'000'000; // If Dim is not 1, fail compilation. This is assumed to be a // mono-dimensional problem - static_assert(Dim::value == 1u); - alpaka::Vec const extents(reserve); + static_assert(dim == 1u); + alpaka::Vec const extents(reserve); + + alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); - alpaka::Vec const layerTilesExtents(static_cast(NLAYERS)); - device_hist_ = std::make_optional(alpaka::allocAsyncBuf(queue_, layerTilesExtents)); - alpaka::Vec const seedsExtents(1u); + /** @todo use queue allocateor */ + device_hist_ = std::make_optional(alpaka::onHost::alloc(device_, layerTilesExtents)); + alpaka::Vec const seedsExtents(1u); device_seeds_ - = std::make_optional(alpaka::allocAsyncBuf, Idx>(queue_, seedsExtents)); + = std::make_optional(alpaka::onHost::alloc>(device_, seedsExtents)); device_followers_ - = std::make_optional(alpaka::allocAsyncBuf, Idx>(queue_, extents)); + = std::make_optional(alpaka::onHost::alloc>(device_, extents)); // INTERNAL VARIABLES RESETTING - alpaka::memset(queue_, device_hist_.value(), 0x0, layerTilesExtents); - alpaka::memset(queue_, device_seeds_.value(), 0x0, seedsExtents); - alpaka::memset(queue_, device_followers_.value(), 0x0, extents); - + alpaka::onHost::memset(queue_, device_hist_.value(), 0x0, layerTilesExtents); + alpaka::onHost::memset(queue_, device_seeds_.value(), 0x0, seedsExtents); + alpaka::onHost::memset(queue_, device_followers_.value(), 0x0, extents); // Set Device Raw Pointers using values from outsice and also internal buffers device_runner_.ptrs_.x = const_cast(x); @@ -931,27 +1055,23 @@ void CLUEAlgoAlpaka::makeClustersCMSSW( device_runner_.ptrs_.isSeed = isSeed; // UPDATE RAW POINTERS FOR INTERNATL DATA STRUCTURES - device_runner_.ptrs_.hist_ = alpaka::getPtrNative(device_hist_.value()); - device_runner_.ptrs_.seeds_ = alpaka::getPtrNative(device_seeds_.value()); - device_runner_.ptrs_.followers_ = alpaka::getPtrNative(device_followers_.value()); + device_runner_.ptrs_.hist_ = alpaka::onHost::data(device_hist_.value()); + device_runner_.ptrs_.seeds_ = alpaka::onHost::data(device_seeds_.value()); + device_runner_.ptrs_.followers_ = alpaka::onHost::data(device_followers_.value()); // Dimension the grid for submission Idx threads_per_block = 256u; - if constexpr(std::is_same_v, alpaka::DevCpu>) - { - threads_per_block = 1u; - } - alpaka::Vec const threadsPerBlock(threads_per_block); - alpaka::Vec const blocksPerGrid(static_cast(ceil(points / (float) threadsPerBlock[0]))); - alpaka::Vec const elementsPerThread(1u); - using WorkDiv = alpaka::WorkDivMembers; - auto const manualWorkDiv = WorkDiv{blocksPerGrid, threadsPerBlock, elementsPerThread}; + alpaka::Vec const threadsPerBlock(threads_per_block); + + alpaka::Vec const blocksPerGrid(alpaka::divExZero(static_cast(points), threadsPerBlock[0])); + + auto const manualWorkDiv = alpaka::onHost::FrameSpec{blocksPerGrid, threadsPerBlock}; // Create the kernel execution tasks. - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeHistogram taskComputeHistogram; - auto const kernelComputeHistogram - = alpaka::createTaskKernel(manualWorkDiv, device_runner_, taskComputeHistogram, points); + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeHistogram taskComputeHistogram; + auto const kernelComputeHistogram = alpaka::KernelBundle(device_runner_, taskComputeHistogram, points); #if ORDER_TILE // printf("Sorting tile for all layers.\n"); @@ -959,32 +1079,29 @@ void CLUEAlgoAlpaka::makeClustersCMSSW( alpaka::Vec const blocksPerGridTile(std::ceil((NLAYERS * T::nTiles) / (float) threadsPerBlockTile[0])); auto const manualWorkDivTile = WorkDiv{blocksPerGridTile, threadsPerBlockTile, elementsPerThread}; - typename CLUEAlgoAlpaka::DeviceRunner::KernelSortHistogram taskSortHistogram; + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelSortHistogram taskSortHistogram; auto const kernelSortHistogram = alpaka::createTaskKernel(manualWorkDivTile, device_runner_, taskSortHistogram); #endif - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeLocalDensity taskComputeLocalDensity; - auto const kernelComputeLocalDensity = (alpaka::createTaskKernel( - manualWorkDiv, - device_runner_, - taskComputeLocalDensity, - dc_, - static_cast(points))); + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeLocalDensity taskComputeLocalDensity; + auto const kernelComputeLocalDensity + = (alpaka::KernelBundle(device_runner_, taskComputeLocalDensity, dc_, static_cast(points))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelComputeDistanceToHigher - taskComputeDistanceToHigher; - auto const kernelComputeDistanceToHigher = (alpaka::createTaskKernel( - manualWorkDiv, + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelComputeDistanceToHigher taskComputeDistanceToHigher; + auto const kernelComputeDistanceToHigher = (alpaka::KernelBundle( device_runner_, taskComputeDistanceToHigher, outlierDeltaFactor_, dc_, static_cast(points))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelFindClustersKappa taskFindClustersKappa; - auto const kernelFindClustersKappa = (alpaka::createTaskKernel( - manualWorkDiv, + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelFindClustersKappa taskFindClustersKappa; + auto const kernelFindClustersKappa = (alpaka::KernelBundle( device_runner_, taskFindClustersKappa, outlierDeltaFactor_, @@ -992,18 +1109,19 @@ void CLUEAlgoAlpaka::makeClustersCMSSW( kappa_, static_cast(points))); - typename CLUEAlgoAlpaka::DeviceRunner::KernelAssignClusters taskAssignClusters; + typename CLUEAlgoAlpaka::DeviceRunner:: + KernelAssignClusters taskAssignClusters; auto const kernelAssignClusters - = (alpaka::createTaskKernel(manualWorkDiv, device_runner_, taskAssignClusters, numberOfClustersScalar)); + = (alpaka::KernelBundle(device_runner_, taskAssignClusters, numberOfClustersScalar)); // Enqueue the kernel execution task - alpaka::enqueue(queue_, kernelComputeHistogram); + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeHistogram); #if ORDER_TILE - alpaka::enqueue(queue_, kernelSortHistogram); + alpaka::onHost::enqueue(queue_, TExecutor{}, manualWorkDiv, kernelSortHistogram); #endif - alpaka::enqueue(queue_, kernelComputeLocalDensity); - alpaka::enqueue(queue_, kernelComputeDistanceToHigher); - alpaka::enqueue(queue_, kernelFindClustersKappa); - alpaka::enqueue(queue_, kernelAssignClusters); + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeLocalDensity); + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelComputeDistanceToHigher); + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelFindClustersKappa); + queue_.enqueue(TExecutor{}, manualWorkDiv, kernelAssignClusters); } diff --git a/clueLib/include/CLUEAlgoGPU.h b/clueLib/include/CLUEAlgoGPU.h index 7c6bd15..f5f872d 100644 --- a/clueLib/include/CLUEAlgoGPU.h +++ b/clueLib/include/CLUEAlgoGPU.h @@ -1,5 +1,4 @@ -#ifndef CLUEAlgoGPU_h -#define CLUEAlgoGPU_h +#pragma once #include #include @@ -363,7 +362,7 @@ __global__ void kernel_calculate_distanceToHigher( { int i = blockIdx.x * blockDim.x + threadIdx.x; - float constat dm = outlierDeltaFactor * dc; + float const dm = outlierDeltaFactor * dc; if(i < numberOfPoints) { @@ -432,7 +431,7 @@ __global__ void kernel_calculate_distanceToHigherTile( if(bin < d_hist[layeri][globalBinOnLayer].size()) { - float constat dm = outlierDeltaFactor * dc; + float const dm = outlierDeltaFactor * dc; int i = d_hist[layeri][globalBinOnLayer][bin]; float deltai = std::numeric_limits::max(); @@ -572,8 +571,8 @@ __global__ void kernel_assign_clusters( cudaStream_t) { int idxCls = blockIdx.x * blockDim.x + threadIdx.x; - auto constto& seeds = d_seeds[0]; - auto constto nSeeds = seeds.size(); + auto const& seeds = d_seeds[0]; + auto const nSeeds = seeds.size(); if(idxCls < nSeeds) { int localStack[localStackSizePerSeed] = {-1}; @@ -674,5 +673,3 @@ void CLUEAlgoGPU::makeClusters() copy_tohost(); CHECK_CUDA_ERROR(cudaStreamSynchronize(stream_)); } - -#endif diff --git a/clueLib/include/GPUVecArray.h b/clueLib/include/GPUVecArray.h index de4f9b0..0f7f7ef 100644 --- a/clueLib/include/GPUVecArray.h +++ b/clueLib/include/GPUVecArray.h @@ -1,5 +1,4 @@ -#ifndef GPUVecArray_h -#define GPUVecArray_h +#pragma once namespace GPU { @@ -165,5 +164,3 @@ namespace GPU }; } // end namespace GPU - -#endif // GPUVecArray_h diff --git a/clueLib/include/GPUVecArrayAlpaka.h b/clueLib/include/GPUVecArrayAlpaka.h index cf98240..b203de1 100644 --- a/clueLib/include/GPUVecArrayAlpaka.h +++ b/clueLib/include/GPUVecArrayAlpaka.h @@ -1,5 +1,5 @@ -#ifndef GPUVecArrayAlpaka_h -#define GPUVecArrayAlpaka_h +#pragma once +#include namespace GPUAlpaka { @@ -57,7 +57,7 @@ namespace GPUAlpaka template ALPAKA_FN_ACC int push_back(T_Acc const& acc, T const& element) { - auto previousSize = atomicAdd(acc, &m_size, 1, alpaka::hierarchy::Blocks{}); + auto previousSize = atomicAdd(acc, &m_size, 1, alpaka::onAcc::scope::block); if(previousSize < maxSize) { m_data[previousSize] = element; @@ -66,7 +66,7 @@ namespace GPUAlpaka else { assert(0); - atomicSub(acc, &m_size, 1, alpaka::hierarchy::Blocks{}); + atomicSub(acc, &m_size, 1, alpaka::onAcc::scope::block); // assert(("Too few elemets reserved", maxSize)); return -1; } @@ -75,7 +75,7 @@ namespace GPUAlpaka template ALPAKA_FN_ACC int emplace_back(T_Acc const& acc, Ts&&... args) { - auto previousSize = atomicAdd(acc, &m_size, 1, alpaka::hierarchy::Blocks{}); + auto previousSize = atomicAdd(acc, &m_size, 1, alpaka::onAcc::scope::block); if(previousSize < maxSize) { (new(&m_data[previousSize]) T(std::forward(args)...)); @@ -84,7 +84,7 @@ namespace GPUAlpaka else { assert(0); - atomicSub(acc, &m_size, 1, alpaka::hierarchy::Blocks{}); + atomicSub(acc, &m_size, 1, alpaka::onAcc::scope::block); return -1; } } @@ -193,5 +193,3 @@ namespace GPUAlpaka }; } // namespace GPUAlpaka - -#endif // GPUVecArray_h diff --git a/clueLib/include/Points.h b/clueLib/include/Points.h index f6746c5..9ea14d0 100644 --- a/clueLib/include/Points.h +++ b/clueLib/include/Points.h @@ -1,5 +1,4 @@ -#ifndef Points_h -#define Points_h +#pragma once struct Points { @@ -17,7 +16,6 @@ struct Points float const* p_weight; float const* p_sigmaNoise; - std::vector rho; std::vector delta; std::vector nearestHigher; @@ -49,4 +47,3 @@ struct Points n = 0; } }; -#endif diff --git a/clueLib/include/Tiles.h b/clueLib/include/Tiles.h index 7674de6..2d5584f 100644 --- a/clueLib/include/Tiles.h +++ b/clueLib/include/Tiles.h @@ -1,6 +1,4 @@ -#ifndef LayerTiles_h -#define LayerTiles_h - +#pragma once #include "TilesConstants.h" #include @@ -88,5 +86,3 @@ class Tiles private: std::vector> tiles_; }; - -#endif // LayerTiles_h diff --git a/clueLib/include/TilesAlpaka.h b/clueLib/include/TilesAlpaka.h index ca734a1..080927b 100644 --- a/clueLib/include/TilesAlpaka.h +++ b/clueLib/include/TilesAlpaka.h @@ -1,6 +1,4 @@ -#ifndef LayerTilesAlpaka_h -#define LayerTilesAlpaka_h - +#pragma once #include "GPUVecArrayAlpaka.h" #include "TilesConstants.h" @@ -9,29 +7,24 @@ #include #include -#if !defined(ALPAKA_ACC_GPU_CUDA_ENABLED) && !defined(ALPAKA_ACC_GPU_HIP_ENABLED) +#if !ALPAKA_LANG_CUDA && !ALPAKA_LANG_HIP struct int4 { int x, y, z, w; }; #endif -template +template class TilesAlpaka { public: using GPUVect = GPUAlpaka::VecArray; // constructor - TilesAlpaka(Acc const& acc) - { - acc_ = acc; - }; - ALPAKA_FN_ACC - void fill(float x, float y, int i) + void fill(auto const& acc, float x, float y, int i) { - tiles_[getGlobalBin(x, y)].push_back(acc_, i); + tiles_[getGlobalBin(x, y)].push_back(acc, i); } ALPAKA_FN_HOST_ACC int getDim1Bin(float x) const @@ -73,10 +66,10 @@ class TilesAlpaka t.reset(); } - ALPAKA_FN_HOST_ACC void sort_unsafe(int i) + ALPAKA_FN_HOST_ACC void sort_unsafe(auto const& acc, int i) { // for (int i = 0; i < T::nTiles; ++i) - tiles_[i].sort_unsafe(acc_); + tiles_[i].sort_unsafe(acc); } ALPAKA_FN_HOST_ACC GPUVect& operator[](int globalBinId) @@ -86,6 +79,4 @@ class TilesAlpaka private: GPUAlpaka::VecArray tiles_; - Acc const& acc_; }; -#endif diff --git a/clueLib/include/TilesConstants.h b/clueLib/include/TilesConstants.h index aa414e8..f3215d9 100644 --- a/clueLib/include/TilesConstants.h +++ b/clueLib/include/TilesConstants.h @@ -1,5 +1,4 @@ -#ifndef TilesConstants_h -#define TilesConstants_h +#pragma once namespace util { @@ -25,5 +24,3 @@ struct TilesConstants static constexpr int nTiles = nColumns * nRows; static constexpr int maxTileDepth = 64; // For accelerators. }; - -#endif // TilesConstants_h diff --git a/clueLib/include/TilesGPU.h b/clueLib/include/TilesGPU.h index 3bebe20..06ffaeb 100644 --- a/clueLib/include/TilesGPU.h +++ b/clueLib/include/TilesGPU.h @@ -1,6 +1,4 @@ -#ifndef LayerTilesGPU_h -#define LayerTilesGPU_h - +#pragma once #include #include #include @@ -75,4 +73,3 @@ class TilesGPU private: GPU::VecArray, T::nTiles> tiles_; }; -#endif diff --git a/plot.py b/plot.py new file mode 100755 index 0000000..1aa1523 --- /dev/null +++ b/plot.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python +import os +import sys +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt + + +def plot_particles(input_data): + df = pd.read_csv(input_data, dtype=np.float32) + df.columns = df.columns.str.strip() + max_clusterid = max(df["clusterId"]) + + plt.figure(dpi=200) + + df_out = df[df.clusterId == -1] # Outliers + plt.scatter(df_out.x, df_out.y, s=10, marker="x", color="0.4") + for i in range(0, int(max_clusterid) + 1): + dfi = df[df.clusterId == i] # ith cluster + plt.scatter(dfi.x, dfi.y, s=10, marker=".") + df_seed = df[df.isSeed == 1] # Only Seeds + plt.scatter(df_seed.x, df_seed.y, s=25, color="r", marker="*") + + plt.xlabel("x", fontsize=14) + plt.ylabel("y", fontsize=14) + + plt.grid() + + plt.savefig("result.svg", format="svg", dpi=1200) + + plt.show() + + +def main(): + if len(sys.argv) != 2: + print("Usage: python script_name.py ") + sys.exit(1) + + filename = sys.argv[1] + + if os.path.isfile(filename): + plot_particles(filename) + + +if __name__ == "__main__": + main() diff --git a/readme.md b/readme.md index be744da..2ef67e0 100644 --- a/readme.md +++ b/readme.md @@ -8,10 +8,10 @@ Z.Chen[1], A. Di Pilato[2,3], F. Pantaleo[4], M. Rovere[4], C. Seez[5] ## 1. Setup -The pre-requisite dependencies are `>=gcc7`, `<=gcc8.3`, `Boost`, `TBB`. Fork this repo if developers. +The pre-requisite dependencies are C++20 compiler. Fork this repo if developers. If CUDA/nvcc are found on the machine, the compilation is performed automatically also for the GPU case. -The path to the nvcc compiler will be automatically taken from the machine. In this case, `>=cuda10` and `<=nvcc11.2` are also required. +The path to the nvcc compiler will be automatically taken from the machine. In this case, `>=cuda12` is required. * **On a CERN machine with GPUs:** Source the LCG View containing GCC, Boost and CUDA: @@ -29,11 +29,8 @@ mkdir install cd build/ ; cmake .. -DCMAKE_INSTALL_PREFIX=../install; make install ``` -* **On an Ubuntu machine with GPUs:** Install Boost and TBB first. +* **On an Ubuntu machine with GPUs:** ```bash -sudo apt-get install libtbb-dev -sudo apt-get install libboost-all-dev - # then setup this project git clone --recurse-submodules https://gitlab.cern.ch/kalos/clue.git cd clue @@ -60,19 +57,28 @@ The test program accept the following parameter from the command line: * `-e sessions`: number of times the clustering algorithm has to run on the same input dataset. That's useful to have a more reliable measure of the timing performance. -* `-t number_TBB_threads`: set the number of TBB threads to be used (when this - makes sense) -* `-u use_accelerator`: enable the GPU version of the executable run. Every +* `-u accelerator`: run with the alpaka executor. You can find valid options with `-U` Every single executable, in fact, has both the CPU and the GPU version embedded. +* `-U`: list enabled alpaka executors. * `-v verbose`: activate verbose output. Among other things, this will also enable the saving of the results of the clustering steps in local text files. If the projects compiles without errors, you can go run the CLUE algorithm by ```bash -./build/src/clue/main -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v -u +# alpaka cpu serial +./build/src/clue_alpaka/mainAlpaka -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v -u CpuSerial + +# alpaka gpu CUDA +./build/src/clue_alpaka/mainAlpaka -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v -u GpuCuda -# in case of only CPU +# alpaka gpu OpenMP +./build/src/clue_alpaka/mainAlpaka -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v -u CpuOmpBlocks + +# in case of original CPU without alpaka ./build/src/clue/main -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v + +# in case of original CUDA without alpaka +./build/src/clue/main -i data/input/aniso_1000.csv -d 7.0 -r 10.0 -o 2 -e 10 -v -u cuda ``` The input files are `data/input/*.csv` with columns diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b34642..36b4d45 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,11 +1,7 @@ if(CMAKE_CUDA_COMPILER) # Native C++ CLUE and native CUDA CLUE add_subdirectory(clue) - - # Native C++ CLUE and CUDA CLUE using ALPAKA - add_subdirectory(clue_cuda_alpaka) - endif(CMAKE_CUDA_COMPILER) -# Native C++ CLUE and TBB CLUE using ALPAKA -#add_subdirectory(clue_tbb_alpaka) +# Native C++ CLUE and CLUE using ALPAKA +add_subdirectory(clue_alpaka) diff --git a/src/clue/CMakeLists.txt b/src/clue/CMakeLists.txt index 9341595..95f5ee3 100644 --- a/src/clue/CMakeLists.txt +++ b/src/clue/CMakeLists.txt @@ -3,22 +3,22 @@ # runtime via a flag. if(CMAKE_CUDA_COMPILER) - set_source_files_properties(${PROJECT_SOURCE_DIR}/src/main.cc PROPERTIES LANGUAGE CUDA) - add_executable(main ${PROJECT_SOURCE_DIR}/src/main.cc) + set_source_files_properties(${PROJECT_SOURCE_DIR}/src/main.cpp PROPERTIES LANGUAGE CUDA) + add_executable(main ${PROJECT_SOURCE_DIR}/src/main.cpp) target_compile_definitions(main PRIVATE) target_compile_options(main PRIVATE - $<$: - -Xcompiler; - -m64; - -expt-relaxed-constexpr; - -w; - > + $<$: + -Xcompiler; + -m64; + -expt-relaxed-constexpr; + -w; + > ) target_link_libraries(main PRIVATE CLUELib) target_include_directories(main PUBLIC - $ + $ ) endif(CMAKE_CUDA_COMPILER) diff --git a/src/clue_alpaka/CMakeLists.txt b/src/clue_alpaka/CMakeLists.txt new file mode 100644 index 0000000..6128124 --- /dev/null +++ b/src/clue_alpaka/CMakeLists.txt @@ -0,0 +1,25 @@ +# mainAlpaka: this will build the native C++ implementation of CLUE and its +# corresponding one built using ALPAKA. Which one to use must be +# selected at runtime via a flag. + +add_executable(mainAlpaka ${PROJECT_SOURCE_DIR}/src/main.cpp) +target_compile_definitions(mainAlpaka PRIVATE USE_ALPAKA) +target_compile_options(mainAlpaka PRIVATE + $<$: + -Xcompiler; + -m64; + -expt-relaxed-constexpr; + -w; +> +) +target_link_libraries(mainAlpaka PRIVATE CLUELib) +target_link_libraries(mainAlpaka PRIVATE alpaka::alpaka) + +target_include_directories(mainAlpaka PUBLIC + $ + PUBLIC + $ + PUBLIC + $ +) +alpaka_finalize(mainAlpaka) diff --git a/src/clue_cuda_alpaka/CMakeLists.txt b/src/clue_cuda_alpaka/CMakeLists.txt deleted file mode 100644 index cd93afd..0000000 --- a/src/clue_cuda_alpaka/CMakeLists.txt +++ /dev/null @@ -1,26 +0,0 @@ -# mainAlpakaCUDA: this will build the native C++ implementation of CLUE and its -# corresponding CUDA one built using ALPAKA. Which one to use must be -# selected at runtime via a flag. - -if(CMAKE_CUDA_COMPILER) - set_source_files_properties(${PROJECT_SOURCE_DIR}/src/main.cc PROPERTIES LANGUAGE CUDA) - add_executable(mainAlpakaCUDA ${PROJECT_SOURCE_DIR}/src/main.cc) - target_compile_definitions(mainAlpakaCUDA PRIVATE USE_ALPAKA FOR_CUDA ALPAKA_ACC_GPU_CUDA_ENABLED=1) - target_compile_options(mainAlpakaCUDA PRIVATE - $<$: - -Xcompiler; - -m64; - -expt-relaxed-constexpr; - -w; - > - ) - target_link_libraries(mainAlpakaCUDA PRIVATE CLUELib) - - target_include_directories(mainAlpakaCUDA PUBLIC - $ - PUBLIC - $ - PUBLIC - $ - ) -endif(CMAKE_CUDA_COMPILER) diff --git a/src/clue_tbb_alpaka/CMakeLists.txt b/src/clue_tbb_alpaka/CMakeLists.txt deleted file mode 100644 index 1ea632a..0000000 --- a/src/clue_tbb_alpaka/CMakeLists.txt +++ /dev/null @@ -1,16 +0,0 @@ -# mainAlpakaTBB: this will build the native C++ implementation of CLUE and its -# corresponding TBB one built using ALPAKA. Which one to use must be -# selected at runtime via a flag. - -set_source_files_properties(${PROJECT_SOURCE_DIR}/src/main.cc PROPERTIES LANGUAGE CXX) -add_executable(mainAlpakaTBB ${PROJECT_SOURCE_DIR}/src/main.cc) -target_compile_definitions(mainAlpakaTBB PRIVATE USE_ALPAKA ALPAKA_ACC_CPU_B_TBB_T_SEQ_ENABLED) -target_link_libraries(mainAlpakaTBB PRIVATE CLUELib TBB::tbb pthread) - -target_include_directories(mainAlpakaTBB PUBLIC - $ - PUBLIC - $ - PUBLIC - $ -) diff --git a/src/main.cc b/src/main.cpp similarity index 72% rename from src/main.cc rename to src/main.cpp index cc8b443..18ef6b0 100644 --- a/src/main.cc +++ b/src/main.cpp @@ -10,17 +10,15 @@ #include #include #include - #if defined(USE_ALPAKA) # include "CLUEAlgoAlpaka.h" + +# include +# include #else # include "CLUEAlgoGPU.h" #endif -#ifdef ALPAKA_ACC_CPU_B_TBB_T_SEQ_ENABLED -# include "tbb/global_control.h" -#endif - #define NLAYERS 100 using namespace std; @@ -58,7 +56,7 @@ pair stats(std::vector const& v) return {m, std::sqrt(sum / den)}; } -void printTimingReport(std::vector& vals, int repeats, const std::string label = "SUMMARY ") +void printTimingReport(std::vector& vals, int repeats, std::string const label = "SUMMARY ") { int precision = 2; float mean = 0.f; @@ -131,6 +129,7 @@ void mainRun( float const rhoc, float const outlierDeltaFactor, bool const use_accelerator, + auto alpakaCfg, int const repeats, bool const verbose) { @@ -209,15 +208,29 @@ void mainRun( #elif defined(USE_ALPAKA) std::cout << "ALPAKA 'Backend' selected" << std::endl; using namespace alpaka; + using namespace alpaka::onHost; // Define the index domain - using Dim = alpaka::DimInt<1u>; - using Idx = uint32_t; - using Acc = SelectedAcc; - CLUEAlgoAlpaka, TilesConstants, NLAYERS> clueAlgo( - dc, - rhoc, - outlierDeltaFactor, - verbose); + + auto deviceSpec = alpakaCfg[object::deviceSpec]; + auto exec = alpakaCfg[object::exec]; + + std::cout << deviceSpec.getApi().getName() << std::endl; + auto devSelector = onHost::makeDeviceSelector(deviceSpec); + onHost::Device computeDevice = devSelector.makeDevice(0); + std::cout << "Using alpaka accelerator: " << onHost::demangledName(exec) << " for " + << deviceSpec.getApi().getName() << std::endl; + Queue queue = computeDevice.makeQueue(); + + Device cpuDevice = makeHostDevice(); + + CLUEAlgoAlpaka< + ALPAKA_TYPEOF(exec), + ALPAKA_TYPEOF(computeDevice), + ALPAKA_TYPEOF(queue), + ALPAKA_TYPEOF(cpuDevice), + TilesConstants, + NLAYERS> + clueAlgo(computeDevice, queue, cpuDevice, dc, rhoc, outlierDeltaFactor, verbose); vals.clear(); for(unsigned r = 0; r < repeats; r++) { @@ -248,7 +261,7 @@ void mainRun( std::cout << "Native CPU(serial) Backend selected" << std::endl; CLUEAlgo clueAlgo(dc, rhoc, outlierDeltaFactor, verbose); vals.clear(); - for(int r = 0; r < repeats; r++) + for(unsigned r = 0; r < repeats; r++) { if(!clueAlgo.setPoints(x.size(), &x[0], &y[0], &layer[0], &weight[0])) exit(EXIT_FAILURE); @@ -257,7 +270,8 @@ void mainRun( clueAlgo.makeClusters(); auto finish = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed = finish - start; - std::cout << "Elapsed time: " << elapsed.count() * 1000 << " ms\n"; + std::cout << "Iteration " << r; + std::cout << " | Elapsed time: " << elapsed.count() * 1000 << " ms\n"; // Skip first event if(r != 0 or repeats == 1) { @@ -286,17 +300,25 @@ int main(int argc, char* argv[]) bool verbose = false; float dc = 20.f, rhoc = 80.f, outlierDeltaFactor = 2.f; int repeats = 10; - int TBBNumberOfThread = 1; int opt; + bool has_outFileName = false; std::string inputFileName; + std::string outputFileName; + std::string alpakaExecutor; + bool list_alpaka_executors = false; - while((opt = getopt(argc, argv, "i:d:r:o:e:t:uv")) != -1) + while((opt = getopt(argc, argv, "i:d:r:o:O:e:u:Uv")) != -1) { switch(opt) { case 'i': /* input filename */ inputFileName = string(optarg); break; + case 'O': /* output filename */ + std::cout << " reached that" << std::endl; + has_outFileName = true; + outputFileName = string(optarg); + break; case 'd': /* delta_c */ dc = stof(string(optarg)); break; @@ -309,46 +331,110 @@ int main(int argc, char* argv[]) case 'e': /* number of repeated session(s) a the selected input file */ repeats = stoi(string(optarg)); break; - case 't': /* number of TBB threads */ - TBBNumberOfThread = stoi(string(optarg)); - std::cout << "Using " << TBBNumberOfThread; - std::cout << " TBB Threads" << std::endl; - break; case 'u': /* Use accelerator */ use_accelerator = true; + alpakaExecutor = string(optarg); + break; + case 'U': /* Use accelerator */ + list_alpaka_executors = true; break; case 'v': /* Verbose output */ verbose = true; break; default: std::cout << "bin/main -i [fileName] -d [dc] -r [rhoc] -o " - "[outlierDeltaFactor] -e [repeats] -t " - "[NumTBBThreads] -u -v" + "[outlierDeltaFactor] -e [repeats] -u [executor] -v" << std::endl; exit(EXIT_FAILURE); } } -#ifdef ALPAKA_ACC_CPU_B_TBB_T_SEQ_ENABLED - if(verbose) - { - std::cout << "Setting up " << TBBNumberOfThread << " TBB Threads" << std::endl; - } - tbb::global_control init(tbb::global_control::max_allowed_parallelism, TBBNumberOfThread); -#endif - ////////////////////////////// // MARK -- set input and output files ////////////////////////////// std::cout << "Input file: " << inputFileName << std::endl; - - std::string outputFileName = create_outputfileName(inputFileName, dc, rhoc, outlierDeltaFactor); + if(has_outFileName) + { + outputFileName = create_outputfileName(outputFileName, dc, rhoc, outlierDeltaFactor); + } + else + { + outputFileName = create_outputfileName(inputFileName, dc, rhoc, outlierDeltaFactor); + } std::cout << "Output file: " << outputFileName << std::endl; ////////////////////////////// // MARK -- test run ////////////////////////////// - mainRun(inputFileName, outputFileName, dc, rhoc, outlierDeltaFactor, use_accelerator, repeats, verbose); + if(use_accelerator) + { +#if defined(USE_ALPAKA) + if(list_alpaka_executors) + { + std::cout << "alpaka executors" << std::endl; + alpaka::onHost::executeForEach( + [&](auto const& cfg) + { + std::cout << " " << alpaka::onHost::getStaticName(cfg[alpaka::object::exec]) << std::endl; + return 0; + }, + alpaka::onHost::allBackends(alpaka::onHost::enabledApis, alpaka::onHost::example::enabledExecutors)); + return 0; + } + + return alpaka::onHost::executeForEach( + [&](auto const& cfg) + { + if(alpakaExecutor == alpaka::onHost::getStaticName(cfg[alpaka::object::exec])) + mainRun( + inputFileName, + outputFileName, + dc, + rhoc, + outlierDeltaFactor, + use_accelerator, + cfg, + repeats, + verbose); + + return 0; + }, + alpaka::onHost::allBackends(alpaka::onHost::enabledApis, alpaka::onHost::example::enabledExecutors)); +#else + mainRun( + inputFileName, + outputFileName, + dc, + rhoc, + outlierDeltaFactor, + use_accelerator, + // dummy, not used if alpaka is disabled + std::make_tuple(1, 1), + repeats, + verbose); +#endif + } + else + { + mainRun( + inputFileName, + outputFileName, + dc, + rhoc, + outlierDeltaFactor, + use_accelerator, + // dummy, not used if alpaka is disabled +#if defined(USE_ALPAKA) + + // select the first valid accelerator, -u is not set therefor we need only a valid configuration + std::get<0>( + alpaka::onHost::allBackends(alpaka::onHost::enabledApis, alpaka::onHost::example::enabledExecutors)), +#else + std::make_tuple(1, 1), +#endif + repeats, + verbose); + } return 0; }