diff --git a/.appveyor.yml b/.appveyor.yml index e0ca042a93..320acde8e1 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -51,7 +51,7 @@ build_script: - "set PATH=%MINICONDA%;%MINICONDA%\\Scripts;%MINICONDA%\\Library\\bin;%PATH%" # install parallelproj and Python stuff # don't do numpy here due to https://github.com/conda-forge/numpy-feedstock/issues/350 - - conda create --name stirbuild -c conda-forge -yq libparallelproj=1 swig pytest ccache ninja cmake + - conda create --name stirbuild -c conda-forge -yq libparallelproj=2 swig pytest ccache ninja cmake - CALL conda.bat activate stirbuild - python --version - pip install numpy matplotlib diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 1022d054ee..ad2ad0e87b 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -140,7 +140,9 @@ jobs: compiler: clang compiler_version: 21 cuda_version: "0" + # openMP off as problems again with brew BUILD_FLAGS: "-DSTIR_OPENMP=OFF" + # parallelproj off as it needs openMP parallelproj: "OFF" BUILD_TYPE: "Release" ROOT: "OFF" @@ -303,16 +305,20 @@ jobs: esac if test "${{matrix.parallelproj}}XX" == "ONXX"; then - git clone --depth 1 --branch v1.7.3 https://github.com/gschramm/parallelproj - mkdir parallelproj/build - cd parallelproj/build + git clone --depth 20 https://github.com/KUL-recon-lab/libparallelproj/ + cd libparallelproj + git checkout v2.0.3 + mkdir build + cd build if test "${{matrix.cuda_version}}" == "0"; then - extra_args="-DSKIP_CUDA_LIB:BOOL=ON" + extra_args="-DUSE_CUDA:BOOL=OFF" + else + extra_args="-DUSE_CUDA:BOOL=ON" fi cmake .. -DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX} ${extra_args} cmake --build . --target install --config Release cd ../.. - rm -rf parallelproj + rm -rf libparallelproj fi # Install ROOT (warning: brittle due to OS versions etc) diff --git a/.gitmodules b/.gitmodules index 5d1ac04565..2e6f41619f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "external_helpers/fmt"] path = external_helpers/fmt url = https://github.com/fmtlib/fmt.git +[submodule "cuvec"] + path = external_helpers/CuVec + url = https://github.com/AMYPAD/CuVec diff --git a/CMakeLists.txt b/CMakeLists.txt index 773c3a5986..f6f6d3c724 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -261,17 +261,18 @@ endif() # Parallelproj if(NOT DISABLE_Parallelproj_PROJECTOR) - find_package(parallelproj 1.3.4 CONFIG) + find_package(parallelproj CONFIG) if (parallelproj_FOUND) set(STIR_WITH_Parallelproj_PROJECTOR ON) - if (parallelproj_built_with_CUDA) + if (PARALLELPROJ_VERSION VERSION_LESS 2.0) + message(WARNING "Found parallelproj ${parallelproj_VERSION}. This will not be supported anymore in the near future.") + endif() + if (PARALLELPROJ_CUDA) message(STATUS "Found parallelproj ${parallelproj_VERSION} (will use its CUDA support)") + set(parallelproj_built_with_CUDA ON) else() message(STATUS "Found parallelproj ${parallelproj_VERSION} (but using its OpenMP version as it wasn't built with CUDA)") endif() - if (parallelproj_VERSION VERSION_LESS 1.0.1) - message(STATUS "If the above parallelproj info looks incorrect, upgrade it to at least 1.0.1 (but 1.2.13 or later is recommended)") - endif() endif() endif() if (STIR_WITH_Parallelproj_PROJECTOR) diff --git a/documentation/release_6.4.htm b/documentation/release_6.4.htm index 996ed7120e..dd36474707 100644 --- a/documentation/release_6.4.htm +++ b/documentation/release_6.4.htm @@ -12,15 +12,14 @@

Summary of changes in STIR release 6.4

Overall summary

- This version is 100% backwards compatible with STIR 6.3, aside from three important bug fixes, - which will change results (see below). + This version is 100% backwards compatible with STIR 6.3.

This is a release with many new features, as listed below. Highlights are

- But there are many other changes from other contibutors as well, see below. + But there are many other changes from other contributors as well, see below. Of course, there is also the usual code-cleanup and some improvements to the documentation. Overall overview and release management was by Kris Thielemans (UCL) with help from Daniel Deidda (NPL). @@ -29,8 +28,8 @@

Overall summary

Patch release info

@@ -61,7 +60,25 @@

Changed functionality


PR #1667. +
  • + Several optimisations were contributed as part of the SyneRBI AI-RBI hackathon: + +
    • +

    Changes to examples

    diff --git a/external_helpers/CMakeLists.txt b/external_helpers/CMakeLists.txt index 49a4ea2437..f03c498d27 100644 --- a/external_helpers/CMakeLists.txt +++ b/external_helpers/CMakeLists.txt @@ -39,3 +39,11 @@ else() set(FMT_INCLUDE_DIRS ${FMT_INCLUDE_DIRS} PARENT_SCOPE) endif() +################ CuVec +if(NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CuVec/Makefile") + message(FATAL_ERROR "The CuVec submodule was not downloaded! GIT_SUBMODULE was turned off or failed. Please update submodules and try again.") +endif() +set(CUVEC_INCLUDE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/CuVec/cuvec/include") +set(CUVEC_INCLUDE_DIR ${CUVEC_INCLUDE_DIR} PARENT_SCOPE) + +install(FILES "${CUVEC_INCLUDE_DIR}/cuvec.cuh" DESTINATION "${STIR_INCLUDE_INSTALL_DIR}/") diff --git a/external_helpers/CuVec b/external_helpers/CuVec new file mode 160000 index 0000000000..f27ececc00 --- /dev/null +++ b/external_helpers/CuVec @@ -0,0 +1 @@ +Subproject commit f27ececc00088db944a30ea8080016b9bdaa8e39 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 14db4ab2a1..aecf58f1e5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -269,6 +269,14 @@ target_link_libraries(stir_registries) # TODO, really should use stir_libs.cmake target_include_directories(stir_registries PRIVATE ${STIR_INCLUDE_DIR}) target_include_directories(stir_registries PRIVATE ${Boost_INCLUDE_DIR}) +target_include_directories(stir_registries PUBLIC + $ + $) +if (NOT STIR_WITH_CUDA) + target_compile_definitions(stir_registries PUBLIC CUVEC_DISABLE_CUDA) +else() + target_link_libraries(stir_registries PUBLIC CUDA::cudart) +endif() # go and look for CMakeLists.txt files in all those directories foreach(STIR_DIR ${STIR_DIRS} ${STIR_TEST_DIRS}) diff --git a/src/include/stir/cuda_utilities.h b/src/include/stir/cuda_utilities.h index cf17eb4321..be7f11e88e 100644 --- a/src/include/stir/cuda_utilities.h +++ b/src/include/stir/cuda_utilities.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2024, University College London + Copyright (C) 2024, 2026, University College London Copyright (C) 2025, University of Milano-Bicocca This file is part of STIR. @@ -22,6 +22,10 @@ #include "stir/Array.h" #include "stir/info.h" #include "stir/error.h" +#ifdef __CUDACC__ +# include +# include "cuvec.cuh" +#endif #include START_NAMESPACE_STIR @@ -45,12 +49,16 @@ struct cuda_int3 int x = 0, y = 0, z = 0; }; #else -# include typedef dim3 cuda_dim3; typedef int3 cuda_int3; #endif #ifdef __CUDACC__ + +//! copy an `Array` to pre-allocated device memory +/*! + \ingroup CUDA +*/ template inline void array_to_device(elemT* dev_data, const Array& stir_array) @@ -71,9 +79,25 @@ array_to_device(elemT* dev_data, const Array& stir_array) } } +//! copy an `Array` to pre-allocated CuVec +/*! + \ingroup CUDA +*/ template inline void -array_to_host(Array& stir_array, const elemT* dev_data) +array_to_device(CuVec& dev_data, const Array& stir_array) +{ + std::copy(stir_array.begin_all(), stir_array.end_all(), dev_data.begin()); +} + +//! copy CUDA pointer to `Array` +/*! + \ingroup CUDA + The third argument is ignored, as `cudaMemcpy` always syncs device and host. +*/ +template +inline void +array_to_host(Array& stir_array, const elemT* dev_data, bool /* sync */ = true) { if (stir_array.is_contiguous()) { @@ -92,6 +116,20 @@ array_to_host(Array& stir_array, const elemT* dev_data) } } +//! copy CuVec to `Array` +/*! + \ingroup CUDA + If \a sync = \c true, the function will call `cudaDeviceSynchronize()` before copying. +*/ +template +inline void +array_to_host(Array& stir_array, const CuVec& dev_data, bool sync = true) +{ + if (sync) + cudaDeviceSynchronize(); + std::copy(dev_data.begin(), dev_data.end(), stir_array.begin_all()); +} + //! \brief Performs a parallel reduction sum on shared memory within a CUDA thread block, final value stored in shared_mem[0]. template __device__ inline void diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 4f94800651..cd760448aa 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -451,35 +451,14 @@ CudaGibbsPenalty_Hessian_Times_Input_kernel(elemT* output, template CudaGibbsPenalty::CudaGibbsPenalty() - : base_type(), // Call parent constructor - d_weights_data(nullptr), - d_kappa_data(nullptr) + : base_type() // Call parent constructor {} template CudaGibbsPenalty::CudaGibbsPenalty(const bool only_2D, float penalization_factor) - : base_type(only_2D, penalization_factor), // Call parent constructor - d_weights_data(nullptr), - d_kappa_data(nullptr) + : base_type(only_2D, penalization_factor) // Call parent constructor {} -template -CudaGibbsPenalty::~CudaGibbsPenalty() -{ - if (d_weights_data) - cudaFree(d_weights_data); - if (d_kappa_data) - cudaFree(d_kappa_data); - if (d_image_data) - cudaFree(d_image_data); - if (d_scalar) - cudaFree(d_scalar); - if (d_output_data) - cudaFree(d_output_data); - if (d_input_data) - cudaFree(d_input_data); -} - template double CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, elemT>& current_image_estimate) @@ -491,17 +470,17 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e if (this->_already_set_up == false) error("CudaGibbsPenalty: set_up has not been called"); - cudaMemset(d_scalar, 0, sizeof(double)); + d_scalar[0] = 0.0; array_to_device(d_image_data, current_image_estimate); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); - if (do_kappa != (!is_null_ptr(d_kappa_data))) + if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_value_kernel<<>>(d_scalar, - d_image_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + CudaGibbsPenalty_value_kernel<<>>(d_scalar.data(), + d_image_data.data(), + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, d_image_dim, d_image_max_indices, @@ -511,9 +490,9 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e this->potential); checkCudaError("compute_value kernel"); - + cudaDeviceSynchronize(); double prior_value; - cudaMemcpy(&prior_value, d_scalar, sizeof(double), cudaMemcpyDeviceToHost); + prior_value = d_scalar[0]; return this->penalisation_factor * prior_value; } @@ -538,13 +517,13 @@ CudaGibbsPenalty::compute_gradient(DiscretisedDensity<3, elem array_to_device(d_image_data, current_image_estimate); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); - if (do_kappa != (!is_null_ptr(d_kappa_data))) + if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_gradient_kernel<<>>(d_output_data, - d_image_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + CudaGibbsPenalty_gradient_kernel<<>>(d_output_data.data(), + d_image_data.data(), + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -555,7 +534,6 @@ CudaGibbsPenalty::compute_gradient(DiscretisedDensity<3, elem this->potential); checkCudaError("compute_gradient kernel"); - cudaDeviceSynchronize(); array_to_host(prior_gradient, d_output_data); // Optional: write gradient to file @@ -576,21 +554,21 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti if (this->penalisation_factor == 0) return 0.; - cudaMemset(d_scalar, 0, sizeof(double)); + d_scalar[0] = 0.0; array_to_device(d_image_data, current_image_estimate); array_to_device(d_input_data, input); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); - if (do_kappa != (!is_null_ptr(this->d_kappa_data))) + if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_dot_input_kernel - <<>>(d_scalar, - d_input_data, - d_image_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + <<>>(d_scalar.data(), + d_input_data.data(), + d_image_data.data(), + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, d_image_dim, d_image_max_indices, @@ -602,7 +580,7 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti checkCudaError("compute_gradient_times_input kernel"); cudaDeviceSynchronize(); double result; - cudaMemcpy(&result, d_scalar, sizeof(double), cudaMemcpyDeviceToHost); + result = d_scalar[0]; return this->penalisation_factor * result; } @@ -626,13 +604,13 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity array_to_device(d_image_data, current_image_estimate); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); - if (do_kappa != (!is_null_ptr(this->d_kappa_data))) + if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data, - d_image_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data.data(), + d_image_data.data(), + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -643,7 +621,6 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity this->potential); checkCudaError("compute_hessian_diagonal kernel"); - array_to_host(Hessian_diag, d_output_data); } @@ -669,14 +646,14 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD array_to_device(d_output_data, output); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); - if (do_kappa != (!is_null_ptr(d_kappa_data))) + if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data, - d_image_data, - d_input_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data.data(), + d_image_data.data(), + d_input_data.data(), + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -687,7 +664,6 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD this->potential); checkCudaError("accumulate_Hessian_times_input kernel"); - array_to_host(output, d_output_data); } @@ -722,41 +698,23 @@ CudaGibbsPenalty::set_up(shared_ptrsize_all() * sizeof(elemT)); - checkCudaError("CudaGibbsPenalty: cudaMalloc for d_image_data"); + d_image_data.resize(target_sptr->size_all()); // Pre-allocate GPU memory for input data - if (d_input_data) - cudaFree(d_input_data); - cudaMalloc(&d_input_data, target_sptr->size_all() * sizeof(elemT)); - checkCudaError("CudaGibbsPenalty: cudaMalloc for d_input_data"); + d_input_data.resize(target_sptr->size_all()); // Pre-allocate GPU memory for outputs - if (d_scalar) - cudaFree(d_scalar); - cudaMalloc(&d_scalar, sizeof(double)); - checkCudaError("CudaGibbsPenalty: cudaMalloc for d_scalar"); + d_scalar.resize(1); - if (d_output_data) - cudaFree(d_output_data); - cudaMalloc(&d_output_data, target_sptr->size_all() * sizeof(elemT)); - checkCudaError("CudaGibbsPenalty: cudaMalloc for d_output_data"); + d_output_data.resize(target_sptr->size_all()); // Copy CPU weights to GPU (weights should already be set up by parent) if (this->weights.get_length() > 0) { - if (d_weights_data) - cudaFree(d_weights_data); - cudaMalloc(&d_weights_data, this->weights.size_all() * sizeof(float)); + d_weights_data.resize(this->weights.size_all()); array_to_device(d_weights_data, this->weights); - checkCudaError("CudaGibbsPenalty: cudaMalloc for d_weights_data"); } - - // Copy CPU kappa image to GPU - if (d_kappa_data) - cudaFree(d_kappa_data); + auto kappa_ptr = this->get_kappa_sptr(); const bool do_kappa = !is_null_ptr(kappa_ptr); @@ -764,7 +722,7 @@ CudaGibbsPenalty::set_up(shared_ptrhas_same_characteristics(*target_sptr)) error("CudaGibbsPenalty: kappa image does not have the same index range as the reconstructed image"); - cudaMalloc(&d_kappa_data, kappa_ptr->size_all() * sizeof(elemT)); + d_kappa_data.resize(kappa_ptr->size_all()); array_to_device(d_kappa_data, *kappa_ptr); } @@ -777,9 +735,7 @@ void CudaGibbsPenalty::set_weights(const Array<3, float>& w) { base_type::set_weights(w); - if (d_weights_data) - cudaFree(d_weights_data); - cudaMalloc(&d_weights_data, w.size_all() * sizeof(float)); + d_weights_data.resize(w.size_all()); array_to_device(d_weights_data, w); } @@ -788,15 +744,16 @@ void CudaGibbsPenalty::set_kappa_sptr(const shared_ptr>& k) { base_type::set_kappa_sptr(k); - if (d_kappa_data) - cudaFree(d_kappa_data); if (!is_null_ptr(k)) { - cudaMalloc(&d_kappa_data, k->size_all() * sizeof(elemT)); + d_kappa_data.resize(k->size_all()); array_to_device(d_kappa_data, *k); } else - d_kappa_data = nullptr; + { + d_kappa_data.resize(0); + d_kappa_data.shrink_to_fit(); + } } END_NAMESPACE_STIR diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h index 7154082bfa..265b702550 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h @@ -24,6 +24,7 @@ #include "stir/cuda_utilities.h" #include "stir/recon_buildblock/GibbsPenalty.h" +#include "cuvec.cuh" #include "stir/shared_ptr.h" #include @@ -62,7 +63,7 @@ class CudaGibbsPenalty : public GibbsPenalty int threads_per_block; size_t shared_mem_bytes; - elemT* d_image_data = nullptr; + mutable CuVec d_image_data; // Currently stir:CartesianCoordinate3D is not supported on GPU, we need a simple structure to store boundaries. cuda_int3 d_image_dim; cuda_int3 d_image_max_indices; @@ -71,22 +72,21 @@ class CudaGibbsPenalty : public GibbsPenalty cuda_int3 d_weight_min_indices; // GPU pointers to weights and kappa data - float* d_weights_data = nullptr; - elemT* d_kappa_data = nullptr; + CuVec d_weights_data; + CuVec d_kappa_data; // Buffers for GPU input/output to avoid reallocating memory on each call see usage in set_up() and ~CudaGibbsPenalty() - mutable double* d_scalar = nullptr; + mutable CuVec d_scalar; // d_scalar is used for compute_value and compute_gradient_times_input as output variable - mutable elemT* d_input_data = nullptr; + mutable CuVec d_input_data; // d_input_data is used for storing input image for compute_gradient_times_input and accumulate_Hessian_times_input - mutable elemT* d_output_data = nullptr; + mutable CuVec d_output_data; // d_output_data is used for storing output image for compute_gradient, accumulate_Hessian_times_input and // compute_Hessian_diagonal public: CudaGibbsPenalty(); CudaGibbsPenalty(const bool only_2D, float penalization_factor); - ~CudaGibbsPenalty(); //! Override CPU version to set up CUDA resources on GPU and call parent set_up Succeeded set_up(shared_ptr> const& target_sptr) override; diff --git a/src/include/stir/recon_buildblock/Parallelproj_projector/ParallelprojHelper.h b/src/include/stir/recon_buildblock/Parallelproj_projector/ParallelprojHelper.h index 47afa9a919..5b369e2603 100644 --- a/src/include/stir/recon_buildblock/Parallelproj_projector/ParallelprojHelper.h +++ b/src/include/stir/recon_buildblock/Parallelproj_projector/ParallelprojHelper.h @@ -21,6 +21,7 @@ #define __stir_recon_buildblock_ParallelprojHelper_h__ #include "stir/common.h" +#include "cuvec.cuh" #include #include @@ -47,8 +48,8 @@ class ParallelprojHelper std::array voxsize; std::array imgdim; std::array origin; - std::vector xstart; - std::vector xend; + CuVec xstart; + CuVec xend; long long num_image_voxel; long long num_lors; diff --git a/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h b/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h index bd83489c0f..47c928bca7 100644 --- a/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h +++ b/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h @@ -53,6 +53,9 @@ class ProjectorByBinPairUsingParallelproj /// Set verbosity void set_verbosity(const bool verbosity); + /// Set number of GPU chunks used by the forward and back projectors. + void set_num_gpu_chunks(const int num_gpu_chunks); + private: shared_ptr _helper; @@ -60,6 +63,7 @@ class ProjectorByBinPairUsingParallelproj void initialise_keymap() override; bool post_processing() override; bool _verbosity; + int _num_gpu_chunks; }; END_NAMESPACE_STIR diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index b0d6466df1..3ddde2de2b 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -167,9 +167,17 @@ endif() endif() if (STIR_WITH_Parallelproj_PROJECTOR) - target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj_c) - if (parallelproj_built_with_CUDA) - target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj_cuda) + if (PARALLELPROJ_VERSION VERSION_LESS 2.0) + # variable to enable work-arounds in .cxx + target_compile_definitions(${TARGET} PRIVATE parallelproj1) + target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj_c) + if (parallelproj_built_with_CUDA) + target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj_cuda) + # variable to enable work-arounds in .cxx + target_compile_definitions(${TARGET} PRIVATE parallelproj1_CUDA) + endif() + else() + target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj) endif() endif() @@ -183,5 +191,14 @@ if (HAVE_HDF5) endif() if (STIR_WITH_CUDA) - target_link_libraries(${TARGET} PRIVATE CUDA::cudart) + target_link_libraries(${TARGET} PUBLIC CUDA::cudart) +endif() + +target_include_directories(${TARGET} PUBLIC #${CUVEC_INCLUDE_DIR}) + $ + $) + + +if (NOT STIR_WITH_CUDA) + target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) endif() diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index fd8ee5d9c6..0107bf251f 100644 --- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx @@ -9,8 +9,10 @@ \author Richard Brown \author Kris Thielemans \author Nicole Jurjew + \author Markus Jehl - Copyright (C) 2019, 2021, 2024 University College London + Copyright (C) 2019, 2021, 2024, 2026 University College London + Copyright (C) 2026 Positrigo This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -30,13 +32,23 @@ #include "stir/ProjDataInMemory.h" #include "stir/LORCoordinates.h" #include "stir/recon_array_functions.h" -#ifdef parallelproj_built_with_CUDA -# include "parallelproj_cuda.h" + +#ifdef parallelproj1 +# ifdef parallelproj_built_with_CUDA +# include "parallelproj_cuda.h" +# define parallelproj1_CUDA +# define joseph3d_back joseph3d_back_cuda +# define joseph3d_tof_sino_back joseph3d_back_tof_sino_cuda +# else +# include "parallelproj_c.h" +# define joseph3d_tof_sino_back joseph3d_back_tof_sino +# endif #else -# include "parallelproj_c.h" +# include "parallelproj.h" #endif -// for debugging, remove later + #include "stir/info.h" +#include "stir/format.h" #include "stir/error.h" #include "stir/stream.h" #include @@ -120,11 +132,12 @@ static void TOF_transpose(std::vector& mem_for_PP_back, const float* STIR_mem, const shared_ptr _helper, - const long long offset) + const long long offset, + const long long num_lors_per_chunk) { const auto num_tof_bins = static_cast(_helper->num_tof_bins); for (unsigned tof_idx = 0; tof_idx < num_tof_bins; ++tof_idx) - for (long long lor_idx = 0; lor_idx < _helper->num_lors; ++lor_idx) + for (long long lor_idx = 0; lor_idx < num_lors_per_chunk; ++lor_idx) { mem_for_PP_back[lor_idx * num_tof_bins + tof_idx] = STIR_mem[offset + tof_idx * _helper->num_lors + lor_idx]; } @@ -158,9 +171,7 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density // create an alias for the projection data const ProjDataInMemory& p(*_proj_data_to_backproject_sptr); - info("Calling parallelproj backprojector", 2); - -#ifdef parallelproj_built_with_CUDA + info(format("Calling parallelproj backprojector with {} chunks", _num_gpu_chunks), 2); long long num_lors_per_chunk_floor = _helper->num_lors / _num_gpu_chunks; long long remainder = _helper->num_lors % _num_gpu_chunks; @@ -168,10 +179,12 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density long long num_lors_per_chunk; long long offset = 0; +#ifdef parallelproj1_CUDA // send image to all visible CUDA devices float** image_on_cuda_devices = copy_float_array_to_all_devices(image_ptr, _helper->num_image_voxel); +#endif - // do (chuck-wise) back projection on the CUDA devices + // do (chunck-wise) back projection for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { if (chunk_num < remainder) @@ -184,52 +197,63 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density } if (p.get_proj_data_info_sptr()->is_tof_data()) { - info("running the CUDA version of parallelproj, about to call function joseph3d_back_tof_sino_cuda for one chunk", 2); - std::vector mem_for_PP_back(num_lors_per_chunk * _helper->num_tof_bins); const float* STIR_mem = p.get_const_data_ptr(); - TOF_transpose(mem_for_PP_back, STIR_mem, _helper, offset); - - // info("created object mem_for_PP_img", 2); - joseph3d_back_tof_sino_cuda(_helper->xend.data() + 3 * offset, - _helper->xstart.data() + 3 * offset, - image_on_cuda_devices, - _helper->origin.data(), - _helper->voxsize.data(), - mem_for_PP_back.data(), // p.get_const_data_ptr() + offset* num_tof_bins, - num_lors_per_chunk, - _helper->imgdim.data(), - _helper->tofbin_width, - &_helper->sigma_tof, - &_helper->tofcenter_offset, - 4, // float n_sigmas - _helper->num_tof_bins, - 0, // unsigned char lor_dependent_sigma_tof - 0, // unsigned char lor_dependent_tofcenter_offset - 64 // threadsperblock + TOF_transpose(mem_for_PP_back, STIR_mem, _helper, offset, num_lors_per_chunk); + + joseph3d_tof_sino_back(_helper->xend.data() + 3 * offset, + _helper->xstart.data() + 3 * offset, +#ifdef parallelproj1_CUDA + image_on_cuda_devices, +#else + image_ptr, +#endif + _helper->origin.data(), + _helper->voxsize.data(), + mem_for_PP_back.data(), + num_lors_per_chunk, + _helper->imgdim.data(), + _helper->tofbin_width, + &_helper->sigma_tof, + &_helper->tofcenter_offset, + 4, // float n_sigmas + _helper->num_tof_bins, + 0, // unsigned char lor_dependent_sigma_tof + 0 // unsigned char lor_dependent_tofcenter_offset +#ifdef parallelproj1_CUDA + , + /* num_threads_per_block */ 64 +#endif ); if (chunk_num != _num_gpu_chunks - 1) p.release_const_data_ptr(); } else { - info("running the CUDA version of parallelproj, about to call function joseph3d_back_cuda for one chunk", 2); - - joseph3d_back_cuda(_helper->xstart.data() + 3 * offset, - _helper->xend.data() + 3 * offset, - image_on_cuda_devices, - _helper->origin.data(), - _helper->voxsize.data(), - p.get_const_data_ptr() + offset, - num_lors_per_chunk, - _helper->imgdim.data(), - /*threadsperblock*/ 64); + joseph3d_back(_helper->xstart.data() + 3 * offset, + _helper->xend.data() + 3 * offset, +#ifdef parallelproj1_CUDA + image_on_cuda_devices, +#else + image_ptr, +#endif + _helper->origin.data(), + _helper->voxsize.data(), + p.get_const_data_ptr() + offset, + num_lors_per_chunk, + _helper->imgdim.data() +#ifdef parallelproj1_CUDA + , + /* num_threads_per_block */ 64 +#endif + ); } info("done", 2); offset += num_lors_per_chunk; } +#ifdef parallelproj1_CUDA // sum backprojected images on the first CUDA device sum_float_arrays_on_first_device(image_on_cuda_devices, _helper->num_image_voxel); @@ -238,46 +262,7 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density // free image array from CUDA devices free_float_array_on_all_devices(image_on_cuda_devices); - -#else - info("Calling parallelproj backprojector (CPU)", 2); - if (this->_proj_data_info_sptr->is_tof_data() == 1) - { - std::vector mem_for_PP_back(_helper->num_lors * _helper->num_tof_bins); - const float* STIR_mem = p.get_const_data_ptr(); - - TOF_transpose(mem_for_PP_back, STIR_mem, _helper, 0); - - joseph3d_back_tof_sino(_helper->xend.data(), - _helper->xstart.data(), - image_ptr, - _helper->origin.data(), - _helper->voxsize.data(), - mem_for_PP_back.data(), - _helper->num_lors, - _helper->imgdim.data(), - _helper->tofbin_width, - &_helper->sigma_tof, - &_helper->tofcenter_offset, - 4, // float n_sigmas, - _helper->num_tof_bins, - 0, // unsigned char lor_dependent_sigma_tof - 0 // unsigned char lor_dependent_tofcenter_offset - ); - } - else - { - joseph3d_back(_helper->xstart.data(), - _helper->xend.data(), - image_ptr, - _helper->origin.data(), - _helper->voxsize.data(), - p.get_const_data_ptr(), - static_cast(p.get_proj_data_info_sptr()->size_all()), - _helper->imgdim.data()); - } #endif - info("PP backprojecting done. Finishing up", 2); p.release_const_data_ptr(); diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index e47862b277..d5a458e403 100644 --- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx @@ -10,7 +10,10 @@ \author Richard Brown \author Kris Thielemans \author Nicole Jurjew - Copyright (C) 2019, 2021, 2024 University College London + \author Markus Jehl + + Copyright (C) 2019, 2021, 2024, 2026 University College London + Copyright (C) 2026 Positrigo This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -25,15 +28,24 @@ #include "stir/ProjDataInfoCylindricalNoArcCorr.h" #include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h" #include "stir/info.h" +#include "stir/format.h" #include "stir/error.h" #include "stir/recon_array_functions.h" #include "stir/utilities.h" #include "stir/TOF_conversions.h" #include -#ifdef parallelproj_built_with_CUDA -# include "parallelproj_cuda.h" + +#ifdef parallelproj1 +# ifdef parallelproj_built_with_CUDA +# include "parallelproj_cuda.h" +# define joseph3d_fwd joseph3d_fwd_cuda +# define joseph3d_tof_sino_fwd joseph3d_fwd_tof_sino_cuda +# else +# include "parallelproj_c.h" +# define joseph3d_tof_sino_fwd joseph3d_fwd_tof_sino +# endif #else -# include "parallelproj_c.h" +# include "parallelproj.h" #endif START_NAMESPACE_STIR @@ -85,14 +97,6 @@ ForwardProjectorByBinParallelproj::set_up(const shared_ptr& check(*proj_data_info_sptr, *_density_sptr); _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr)); -#if 0 - shared_ptr - proj_data_info_cy_no_ar_cor_sptr( - dynamic_pointer_cast( - proj_data_info_sptr)); - if (is_null_ptr(proj_data_info_cy_no_ar_cor_sptr)) - error("ForwardProjectorByBinParallelproj: Failed casting to ProjDataInfoCylindricalNoArcCorr"); -#endif // Initialise projected_data_sptr from this->_proj_data_info_sptr _projected_data_sptr.reset(new ProjDataInMemory(this->_density_sptr->get_exam_info_sptr(), proj_data_info_sptr)); if (!this->_do_not_setup_helper) @@ -167,14 +171,7 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& image_ptr = image_vec.data(); } -#if 0 - // needed to set output to zero as parallelproj accumulates but is no longer the case - _projected_data_sptr->fill(0.F); -#endif - -#ifdef parallelproj_built_with_CUDA - - info("Calling parallelproj forward (CUDA)", 2); + info(format("Calling parallelproj forward with {} chunks", _num_gpu_chunks), 2); long long num_lors_per_chunk_floor = _helper->num_lors @@ -184,10 +181,12 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& long long num_lors_per_chunk; long long offset = 0; +#ifdef parallelproj1_CUDA // send image to all visible CUDA devices float** image_on_cuda_devices = copy_float_array_to_all_devices(image_ptr, _helper->num_image_voxel); +#endif - // do (chuck-wise) projection on the CUDA devices + // do (chunk-wise) projection for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { if (chunk_num < remainder) @@ -203,22 +202,29 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& { std::vector mem_for_PP(num_lors_per_chunk * _helper->num_tof_bins); - joseph3d_fwd_tof_sino_cuda(_helper->xend.data() + 3 * offset, - _helper->xstart.data() + 3 * offset, - image_on_cuda_devices, - _helper->origin.data(), - _helper->voxsize.data(), - mem_for_PP.data(), // this is where the data is written to - num_lors_per_chunk, // PP docu: "number of geometrical LORs"; - _helper->imgdim.data(), - _helper->tofbin_width, - &_helper->sigma_tof, - &_helper->tofcenter_offset, - 4, // float n_sigmas - _helper->num_tof_bins, // short n_tofbins - 0, // unsigned char lor_dependent_sigma_tof - 0, // unsigned char lor_dependent_tofcenter_offset - 64 // threadsperblock + joseph3d_tof_sino_fwd(_helper->xend.data() + 3 * offset, + _helper->xstart.data() + 3 * offset, +#ifdef parallelproj1_CUDA + image_on_cuda_devices, +#else + image_ptr, +#endif + _helper->origin.data(), + _helper->voxsize.data(), + mem_for_PP.data(), // this is where the data is written to + num_lors_per_chunk, // PP docu: "number of geometrical LORs"; + _helper->imgdim.data(), + _helper->tofbin_width, + &_helper->sigma_tof, + &_helper->tofcenter_offset, + 4, // float n_sigmas + _helper->num_tof_bins, // short n_tofbins + 0, // unsigned char lor_dependent_sigma_tof + 0 // unsigned char lor_dependent_tofcenter_offset +#ifdef parallelproj1_CUDA + , + 64 // threadsperblock +#endif ); float* STIR_mem = _projected_data_sptr->get_data_ptr(); @@ -227,69 +233,36 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& if (chunk_num != _num_gpu_chunks - 1) _projected_data_sptr->release_data_ptr(); - info("current proj max: " - + std::to_string(*std::max_element(_projected_data_sptr->begin(), _projected_data_sptr->end()))); } else { - joseph3d_fwd_cuda(_helper->xstart.data() + 3 * offset, - _helper->xend.data() + 3 * offset, - image_on_cuda_devices, - _helper->origin.data(), - _helper->voxsize.data(), - _projected_data_sptr->get_data_ptr() + offset, - num_lors_per_chunk, - _helper->imgdim.data(), - /*threadsperblock*/ 64); + joseph3d_fwd(_helper->xstart.data() + 3 * offset, + _helper->xend.data() + 3 * offset, +#ifdef parallelproj1_CUDA + image_on_cuda_devices, +#else + image_ptr, +#endif + _helper->origin.data(), + _helper->voxsize.data(), + _projected_data_sptr->get_data_ptr() + offset, + num_lors_per_chunk, + _helper->imgdim.data() +#ifdef parallelproj1_CUDA + , + 64 // threadsperblock +#endif + ); if (chunk_num != _num_gpu_chunks - 1) _projected_data_sptr->release_data_ptr(); } offset += num_lors_per_chunk; } +#ifdef parallelproj1_CUDA // free image array from CUDA devices free_float_array_on_all_devices(image_on_cuda_devices); - -#else - info("Calling parallelproj forward (CPU)", 2); - - if (this->_proj_data_info_sptr->is_tof_data() == 1) - { - - std::vector mem_for_PP(_helper->num_lors * _helper->num_tof_bins); - joseph3d_fwd_tof_sino(_helper->xend.data(), - _helper->xstart.data(), - image_ptr, - _helper->origin.data(), - _helper->voxsize.data(), - mem_for_PP.data(), - _helper->num_lors, - _helper->imgdim.data(), - _helper->tofbin_width, - &_helper->sigma_tof, - &_helper->tofcenter_offset, - 4, // float n_sigmas, - _helper->num_tof_bins, - 0, // unsigned char lor_dependent_sigma_tof - 0 // unsigned char lor_dependent_tofcenter_offset - ); - - float* STIR_mem = _projected_data_sptr->get_data_ptr(); - TOF_transpose(STIR_mem, mem_for_PP, _helper, 0, _helper->num_lors); - } - else - { - joseph3d_fwd(_helper->xstart.data(), - _helper->xend.data(), - image_ptr, - _helper->origin.data(), - _helper->voxsize.data(), - _projected_data_sptr->get_data_ptr(), - static_cast(_projected_data_sptr->get_proj_data_info_sptr()->size_all()), - _helper->imgdim.data()); - } #endif - info("done", 2); if (_density_sptr->is_contiguous()) { diff --git a/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx index ea35717f35..c656ef3cf3 100644 --- a/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx @@ -34,6 +34,7 @@ ProjectorByBinPairUsingParallelproj::initialise_keymap() parser.add_start_key("Projector Pair Using Parallelproj Parameters"); parser.add_stop_key("End Projector Pair Using Parallelproj Parameters"); parser.add_key("verbosity", &_verbosity); + parser.add_key("num_gpu_chunks", &_num_gpu_chunks); } void @@ -41,12 +42,14 @@ ProjectorByBinPairUsingParallelproj::set_defaults() { base_type::set_defaults(); this->set_verbosity(true); + this->set_num_gpu_chunks(1); } bool ProjectorByBinPairUsingParallelproj::post_processing() { this->set_verbosity(this->_verbosity); + this->set_num_gpu_chunks(this->_num_gpu_chunks); if (base_type::post_processing()) return true; @@ -98,4 +101,20 @@ ProjectorByBinPairUsingParallelproj::set_verbosity(const bool verbosity) bck_prj_downcast_sptr->set_verbosity(_verbosity); } +void +ProjectorByBinPairUsingParallelproj::set_num_gpu_chunks(const int num_gpu_chunks) +{ + _num_gpu_chunks = num_gpu_chunks; + + shared_ptr fwd_prj_downcast_sptr + = dynamic_pointer_cast(this->forward_projector_sptr); + if (fwd_prj_downcast_sptr) + fwd_prj_downcast_sptr->set_num_gpu_chunks(_num_gpu_chunks); + + shared_ptr bck_prj_downcast_sptr + = dynamic_pointer_cast(this->back_projector_sptr); + if (bck_prj_downcast_sptr) + bck_prj_downcast_sptr->set_num_gpu_chunks(_num_gpu_chunks); +} + END_NAMESPACE_STIR