From b3532a9f4a48d577e4914a0d7e4f8146d9ae499e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sat, 7 Mar 2026 15:10:44 +0000 Subject: [PATCH 01/27] convert to parallelproj 2.0 (WIP) --- CMakeLists.txt | 8 +++----- src/recon_buildblock/CMakeLists.txt | 5 +---- .../BackProjectorByBinParallelproj.cxx | 14 ++++++++------ .../ForwardProjectorByBinParallelproj.cxx | 15 +++++++++------ 4 files changed, 21 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 773c3a5986..67c584dd0b 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -261,17 +261,15 @@ endif() # Parallelproj if(NOT DISABLE_Parallelproj_PROJECTOR) - find_package(parallelproj 1.3.4 CONFIG) + find_package(parallelproj 2.0 CONFIG) if (parallelproj_FOUND) set(STIR_WITH_Parallelproj_PROJECTOR ON) - if (parallelproj_built_with_CUDA) + 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/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index b0d6466df1..8e13d532bd 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -167,10 +167,7 @@ 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) - endif() + target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj) endif() if (HAVE_HDF5) diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index fd8ee5d9c6..bf3279beff 100644 --- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx @@ -30,11 +30,13 @@ #include "stir/ProjDataInMemory.h" #include "stir/LORCoordinates.h" #include "stir/recon_array_functions.h" -#ifdef parallelproj_built_with_CUDA -# include "parallelproj_cuda.h" -#else -# include "parallelproj_c.h" -#endif + +//#ifdef parallelproj_built_with_CUDA +//# include "parallelproj_cuda.h" +//#else +//# include "parallelproj_c.h" +//#endif +# include "parallelproj.h" // for debugging, remove later #include "stir/info.h" #include "stir/error.h" @@ -248,7 +250,7 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density TOF_transpose(mem_for_PP_back, STIR_mem, _helper, 0); - joseph3d_back_tof_sino(_helper->xend.data(), + joseph3d_tof_sino_back(_helper->xend.data(), _helper->xstart.data(), image_ptr, _helper->origin.data(), diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index e47862b277..2cf6961909 100644 --- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx @@ -30,11 +30,14 @@ #include "stir/utilities.h" #include "stir/TOF_conversions.h" #include -#ifdef parallelproj_built_with_CUDA -# include "parallelproj_cuda.h" -#else -# include "parallelproj_c.h" -#endif +//#if 0 +//#ifdef parallelproj_built_with_CUDA +//# include "parallelproj_cuda.h" +//#else +//# include "parallelproj_c.h" +//#endif +# include "parallelproj.h" +//#endif START_NAMESPACE_STIR @@ -257,7 +260,7 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& { std::vector mem_for_PP(_helper->num_lors * _helper->num_tof_bins); - joseph3d_fwd_tof_sino(_helper->xend.data(), + joseph3d_tof_sino_fwd(_helper->xend.data(), _helper->xstart.data(), image_ptr, _helper->origin.data(), From b357df693515da51ac47bc7c0a8f402c9a750b65 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Sat, 7 Mar 2026 16:51:35 +0000 Subject: [PATCH 02/27] [GHA] use parallelproj 2.0 --- .github/workflows/build-test.yml | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index aec0040802..4dfb046344 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -303,16 +303,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) From 87b207e0cee53de2ad01e2d10eebbbca5b99e907 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 11 Mar 2026 14:50:43 +0000 Subject: [PATCH 03/27] [GHA] MacOS clang21: disable OpenMP and enable parallelproj --- .github/workflows/build-test.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 4dfb046344..ae2edbce80 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -140,8 +140,8 @@ jobs: compiler: clang compiler_version: 21 cuda_version: "0" - BUILD_FLAGS: "-DSTIR_OPENMP=ON" - parallelproj: "OFF" + BUILD_FLAGS: "-DSTIR_OPENMP=OFF" + parallelproj: "ON" BUILD_TYPE: "Release" ROOT: "OFF" ITK: "OFF" From 87fa48cc255492266d0d393c20f027d4b120c821 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 11 Mar 2026 15:47:59 +0000 Subject: [PATCH 04/27] added cuvec submodule --- .gitmodules | 3 +++ cuvec | 1 + 2 files changed, 4 insertions(+) create mode 160000 cuvec diff --git a/.gitmodules b/.gitmodules index 5d1ac04565..3fb69c6576 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 = cuvec + url = https://github.com/AMYPAD/cuvec diff --git a/cuvec b/cuvec new file mode 160000 index 0000000000..f27ececc00 --- /dev/null +++ b/cuvec @@ -0,0 +1 @@ +Subproject commit f27ececc00088db944a30ea8080016b9bdaa8e39 From 858ddbf27009940e9f71be0a8362916fdd38a7eb Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 11 Mar 2026 15:59:02 +0000 Subject: [PATCH 05/27] add cuvec include (TEMP COMMIT) --- src/recon_buildblock/CMakeLists.txt | 3 +++ .../Parallelproj_projector/BackProjectorByBinParallelproj.cxx | 2 +- .../ForwardProjectorByBinParallelproj.cxx | 2 +- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 8e13d532bd..a96d2a6e90 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -182,3 +182,6 @@ endif() if (STIR_WITH_CUDA) target_link_libraries(${TARGET} PRIVATE CUDA::cudart) endif() + +target_link_libraries(${TARGET} PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) + diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index bf3279beff..025a245e04 100644 --- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx @@ -36,7 +36,7 @@ //#else //# include "parallelproj_c.h" //#endif -# include "parallelproj.h" +#include "parallelproj.h" // for debugging, remove later #include "stir/info.h" #include "stir/error.h" diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index 2cf6961909..b11219b2a5 100644 --- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx @@ -36,7 +36,7 @@ //#else //# include "parallelproj_c.h" //#endif -# include "parallelproj.h" +#include "parallelproj.h" //#endif START_NAMESPACE_STIR From 62c3e0314a4fd99c4facd5dc3d365f23b6832f78 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 11 Mar 2026 16:13:22 +0000 Subject: [PATCH 06/27] fix cuvec include --- src/recon_buildblock/CMakeLists.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index a96d2a6e90..82a09f9e51 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -183,5 +183,7 @@ if (STIR_WITH_CUDA) target_link_libraries(${TARGET} PRIVATE CUDA::cudart) endif() -target_link_libraries(${TARGET} PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) - +target_include_directories(${TARGET} PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) +if (NOT STIR_WITH_CUDA) + target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) +endif() From a4005a54acc45b69753d059a0d61969677aacf62 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 11 Mar 2026 18:20:04 +0000 Subject: [PATCH 07/27] Use cuvec in CudaGibbsPenalty --- src/include/stir/cuda_utilities.h | 21 +++++++++++++++++-- .../CUDA/CudaGibbsPenalty.cuh | 17 ++++++--------- .../recon_buildblock/CUDA/CudaGibbsPenalty.h | 3 ++- 3 files changed, 27 insertions(+), 14 deletions(-) diff --git a/src/include/stir/cuda_utilities.h b/src/include/stir/cuda_utilities.h index cf17eb4321..a8ff0eba20 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,7 +49,6 @@ struct cuda_int3 int x = 0, y = 0, z = 0; }; #else -# include typedef dim3 cuda_dim3; typedef int3 cuda_int3; #endif @@ -71,6 +74,13 @@ array_to_device(elemT* dev_data, const Array& stir_array) } } +template +inline void +array_to_device(CuVec& dev_data, const Array& stir_array) +{ + std::copy(stir_array.begin_all(), stir_array.end_all(), dev_data.begin()); +} + template inline void array_to_host(Array& stir_array, const elemT* dev_data) @@ -92,6 +102,13 @@ array_to_host(Array& stir_array, const elemT* dev_data) } } +template +inline void +array_to_host(Array& stir_array, const CuVec& dev_data) +{ + 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..647d1ead88 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -470,8 +470,6 @@ CudaGibbsPenalty::~CudaGibbsPenalty() 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) @@ -499,7 +497,7 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_value_kernel<<>>(d_scalar, - d_image_data, + d_image_data.data(), d_weights_data, do_kappa ? d_kappa_data : nullptr, do_kappa, @@ -542,7 +540,7 @@ CudaGibbsPenalty::compute_gradient(DiscretisedDensity<3, elem error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_kernel<<>>(d_output_data, - d_image_data, + d_image_data.data(), d_weights_data, do_kappa ? d_kappa_data : nullptr, do_kappa, @@ -588,7 +586,7 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti CudaGibbsPenalty_gradient_dot_input_kernel <<>>(d_scalar, d_input_data, - d_image_data, + d_image_data.data(), d_weights_data, do_kappa ? d_kappa_data : nullptr, do_kappa, @@ -630,7 +628,7 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data, - d_image_data, + d_image_data.data(), d_weights_data, do_kappa ? d_kappa_data : nullptr, do_kappa, @@ -673,7 +671,7 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data, - d_image_data, + d_image_data.data(), d_input_data, d_weights_data, do_kappa ? d_kappa_data : nullptr, @@ -722,10 +720,7 @@ 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) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h index 7154082bfa..90b5d1184b 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; + 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; From 572c8a3f854fecb8cbbd92c154961451f6ce2c54 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 12 Mar 2026 00:26:04 +0000 Subject: [PATCH 08/27] [CMake] add cuvec path to registries --- src/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 14db4ab2a1..72280f63f1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -269,6 +269,10 @@ 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 PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) +if (NOT STIR_WITH_CUDA) + target_compile_definitions(stir_registries PUBLIC CUVEC_DISABLE_CUDA) +endif() # go and look for CMakeLists.txt files in all those directories foreach(STIR_DIR ${STIR_DIRS} ${STIR_TEST_DIRS}) From e6cf2625ed2a3ba5cff7063b398666940cbefe21 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 12 Mar 2026 00:27:48 +0000 Subject: [PATCH 09/27] [CMake] target_include_directories for cuvec PUBLIC --- src/recon_buildblock/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 82a09f9e51..6a684eb642 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -183,7 +183,7 @@ if (STIR_WITH_CUDA) target_link_libraries(${TARGET} PRIVATE CUDA::cudart) endif() -target_include_directories(${TARGET} PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) +target_include_directories(${TARGET} PUBLIC ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) if (NOT STIR_WITH_CUDA) target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) endif() From 5ca02e39deabc30f89611c90f3aa57df1375427a Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 12 Mar 2026 09:45:29 +0000 Subject: [PATCH 10/27] [CMake] install cuvec.cuh and fix CMake includes still WIP --- src/CMakeLists.txt | 7 ++++++- src/recon_buildblock/CMakeLists.txt | 8 +++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 72280f63f1..96f262d935 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -262,6 +262,9 @@ if (STIR_ENABLE_EXPERIMENTAL) ADD_SUBDIRECTORY(${CMAKE_SOURCE_DIR}/src/experimental) endif() +set(CUVEC_INCLUDE_DIR "${CMAKE_SOURCE_DIR}/cuvec/cuvec/include") + +install(FILES "${CUVEC_INCLUDE_DIR}/cuvec.cuh" DESTINATION "${STIR_INCLUDE_INSTALL_DIR}/") # create an object library with all the "registries" # Warning: dependencies for object libraries are not transitive! add_library(stir_registries OBJECT ${STIR_REGISTRIES}) @@ -269,7 +272,9 @@ 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 PRIVATE ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) +target_include_directories(stir_registries PUBLIC + $ + $) if (NOT STIR_WITH_CUDA) target_compile_definitions(stir_registries PUBLIC CUVEC_DISABLE_CUDA) endif() diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 6a684eb642..a2b035898f 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -183,7 +183,13 @@ if (STIR_WITH_CUDA) target_link_libraries(${TARGET} PRIVATE CUDA::cudart) endif() -target_include_directories(${TARGET} PUBLIC ${CMAKE_SOURCE_DIR}/cuvec/cuvec/include) +get_target_property(CHECK_CUVEC ${TARGET} INTERFACE_INCLUDE_DIRECTORIES) +message(STATUS "INC ${CHECK_CUVEC}") + +# TODO really should be PUBLIC +target_include_directories(${TARGET} PRIVATE ${CUVEC_INCLUDE_DIR}) + +get_target_property(CHECK_CUVEC ${TARGET} INTERFACE_INCLUDE_DIRECTORIES) if (NOT STIR_WITH_CUDA) target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) endif() From 3eb37ac6eebe246e7e1efbdb8136fc69719cb16e Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 12 Mar 2026 10:03:10 +0000 Subject: [PATCH 11/27] make CuVeC member mutable for writing --- src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h index 90b5d1184b..3fcfe83d49 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h @@ -63,7 +63,7 @@ class CudaGibbsPenalty : public GibbsPenalty int threads_per_block; size_t shared_mem_bytes; - CuVec d_image_data; + 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; From cf7b7afac2a456c1dde26327fa517f65591c896b Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 12 Mar 2026 10:04:31 +0000 Subject: [PATCH 12/27] fix target_include_directories for cuvec --- src/recon_buildblock/CMakeLists.txt | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index a2b035898f..69e4489c17 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -183,13 +183,11 @@ if (STIR_WITH_CUDA) target_link_libraries(${TARGET} PRIVATE CUDA::cudart) endif() -get_target_property(CHECK_CUVEC ${TARGET} INTERFACE_INCLUDE_DIRECTORIES) -message(STATUS "INC ${CHECK_CUVEC}") +target_include_directories(${TARGET} PUBLIC #${CUVEC_INCLUDE_DIR}) + $ + $) -# TODO really should be PUBLIC -target_include_directories(${TARGET} PRIVATE ${CUVEC_INCLUDE_DIR}) -get_target_property(CHECK_CUVEC ${TARGET} INTERFACE_INCLUDE_DIRECTORIES) if (NOT STIR_WITH_CUDA) target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) endif() From 8bb779337f5f710231cbd5d37b9a31b797eb2aa7 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Thu, 12 Mar 2026 13:08:04 +0000 Subject: [PATCH 13/27] Make CUDA dependency public. --- src/CMakeLists.txt | 2 ++ src/recon_buildblock/CMakeLists.txt | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 96f262d935..dad8574fd5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -277,6 +277,8 @@ 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 diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 69e4489c17..b2dae8ba85 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -180,7 +180,7 @@ 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}) From 2dd762e6a7b0c96a98e1de88b21d12c66556c0f7 Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 13:02:57 +0000 Subject: [PATCH 14/27] extend CuVec use in CudaGibbsPenalty --- .../CUDA/CudaGibbsPenalty.cuh | 65 ++++++------------- .../recon_buildblock/CUDA/CudaGibbsPenalty.h | 4 +- 2 files changed, 21 insertions(+), 48 deletions(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 647d1ead88..3f193914d5 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -451,25 +451,17 @@ 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_scalar) cudaFree(d_scalar); if (d_output_data) @@ -493,13 +485,11 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e 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))) - error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_value_kernel<<>>(d_scalar, d_image_data.data(), - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, d_image_dim, d_image_max_indices, @@ -536,13 +526,11 @@ 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))) - error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_kernel<<>>(d_output_data, d_image_data.data(), - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -580,15 +568,13 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti 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))) - error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_dot_input_kernel <<>>(d_scalar, d_input_data, d_image_data.data(), - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, d_image_dim, d_image_max_indices, @@ -624,13 +610,11 @@ 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))) - error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data, d_image_data.data(), - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -667,14 +651,12 @@ 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))) - error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data, d_image_data.data(), d_input_data, - d_weights_data, - do_kappa ? d_kappa_data : nullptr, + d_weights_data.data(), + do_kappa ? d_kappa_data.data() : nullptr, do_kappa, this->penalisation_factor, d_image_dim, @@ -742,16 +724,10 @@ CudaGibbsPenalty::set_up(shared_ptrweights.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); @@ -759,7 +735,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); } @@ -772,9 +748,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); } @@ -783,15 +757,14 @@ 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 3fcfe83d49..1c3f0b06b2 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h @@ -72,8 +72,8 @@ 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; From a0e3f149ff56536b804dd0ab1ffe11b92a9a6a12 Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 13:31:01 +0000 Subject: [PATCH 15/27] restoring back the consistency check --- .../stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 3f193914d5..8046f08cdf 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -485,6 +485,8 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e array_to_device(d_image_data, current_image_estimate); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); + if (do_kappa != (!d_kappa_data.empty())) + error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_value_kernel<<>>(d_scalar, d_image_data.data(), @@ -526,6 +528,8 @@ 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 != (!d_kappa_data.empty())) + error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_kernel<<>>(d_output_data, d_image_data.data(), @@ -568,6 +572,8 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti array_to_device(d_input_data, input); const bool do_kappa = !is_null_ptr(this->get_kappa_sptr()); + if (do_kappa != (!d_kappa_data.empty())) + error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_dot_input_kernel <<>>(d_scalar, @@ -610,6 +616,8 @@ 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 != (!d_kappa_data.empty())) + error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data, d_image_data.data(), @@ -651,6 +659,8 @@ 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 != (!d_kappa_data.empty())) + error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data, d_image_data.data(), From 13741ef64799f1e2bab99012a1b51467294f344b Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 14:40:41 +0000 Subject: [PATCH 16/27] extending cuvec use to d_scalar --- .../CUDA/CudaGibbsPenalty.cuh | 19 +++++++------------ .../recon_buildblock/CUDA/CudaGibbsPenalty.h | 2 +- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 8046f08cdf..9a16020b15 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -462,8 +462,6 @@ CudaGibbsPenalty::CudaGibbsPenalty(const bool only_2D, float template CudaGibbsPenalty::~CudaGibbsPenalty() { - if (d_scalar) - cudaFree(d_scalar); if (d_output_data) cudaFree(d_output_data); if (d_input_data) @@ -481,14 +479,14 @@ 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 != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_value_kernel<<>>(d_scalar, + CudaGibbsPenalty_value_kernel<<>>(d_scalar.data(), d_image_data.data(), d_weights_data.data(), do_kappa ? d_kappa_data.data() : nullptr, @@ -503,7 +501,7 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e checkCudaError("compute_value kernel"); double prior_value; - cudaMemcpy(&prior_value, d_scalar, sizeof(double), cudaMemcpyDeviceToHost); + prior_value = d_scalar[0]; return this->penalisation_factor * prior_value; } @@ -566,7 +564,7 @@ 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); @@ -576,7 +574,7 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); CudaGibbsPenalty_gradient_dot_input_kernel - <<>>(d_scalar, + <<>>(d_scalar.data(), d_input_data, d_image_data.data(), d_weights_data.data(), @@ -592,7 +590,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; } @@ -721,10 +719,7 @@ CudaGibbsPenalty::set_up(shared_ptr 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; // d_input_data is used for storing input image for compute_gradient_times_input and accumulate_Hessian_times_input From 2139312c2b412b0849aaa84844cb69cc91382ec7 Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 14:50:58 +0000 Subject: [PATCH 17/27] extending cuvec use to d_input and d_output --- .../CUDA/CudaGibbsPenalty.cuh | 25 ++++++------------- .../recon_buildblock/CUDA/CudaGibbsPenalty.h | 4 +-- 2 files changed, 10 insertions(+), 19 deletions(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 9a16020b15..9ea6ba64f0 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -462,10 +462,7 @@ CudaGibbsPenalty::CudaGibbsPenalty(const bool only_2D, float template CudaGibbsPenalty::~CudaGibbsPenalty() { - if (d_output_data) - cudaFree(d_output_data); - if (d_input_data) - cudaFree(d_input_data); + // Device vectors will automatically free their memory when they go out of scope } template @@ -529,7 +526,7 @@ CudaGibbsPenalty::compute_gradient(DiscretisedDensity<3, elem if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_gradient_kernel<<>>(d_output_data, + CudaGibbsPenalty_gradient_kernel<<>>(d_output_data.data(), d_image_data.data(), d_weights_data.data(), do_kappa ? d_kappa_data.data() : nullptr, @@ -575,7 +572,7 @@ CudaGibbsPenalty::compute_gradient_times_input(const Discreti CudaGibbsPenalty_gradient_dot_input_kernel <<>>(d_scalar.data(), - d_input_data, + d_input_data.data(), d_image_data.data(), d_weights_data.data(), do_kappa ? d_kappa_data.data() : nullptr, @@ -617,7 +614,7 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data, + CudaGibbsPenalty_Hessian_diagonal_kernel<<>>(d_output_data.data(), d_image_data.data(), d_weights_data.data(), do_kappa ? d_kappa_data.data() : nullptr, @@ -660,9 +657,9 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD if (do_kappa != (!d_kappa_data.empty())) error("CudaGibbsPenalty internal error: inconsistent CPU and device kappa"); - CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data, + CudaGibbsPenalty_Hessian_Times_Input_kernel<<>>(d_output_data.data(), d_image_data.data(), - d_input_data, + d_input_data.data(), d_weights_data.data(), do_kappa ? d_kappa_data.data() : nullptr, do_kappa, @@ -713,18 +710,12 @@ CudaGibbsPenalty::set_up(shared_ptrsize_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 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) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h index 1ea897ad54..cc4985bc80 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h @@ -78,9 +78,9 @@ class CudaGibbsPenalty : public GibbsPenalty // Buffers for GPU input/output to avoid reallocating memory on each call see usage in set_up() and ~CudaGibbsPenalty() 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 From 59af8113475155475ee316620a88ce2de2a2ab09 Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 14:56:11 +0000 Subject: [PATCH 18/27] removing destructor as redundant after switch to cuvec --- src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh | 6 ------ src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h | 1 - 2 files changed, 7 deletions(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 9ea6ba64f0..82b9f148bb 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -459,12 +459,6 @@ CudaGibbsPenalty::CudaGibbsPenalty(const bool only_2D, float : base_type(only_2D, penalization_factor) // Call parent constructor {} -template -CudaGibbsPenalty::~CudaGibbsPenalty() -{ - // Device vectors will automatically free their memory when they go out of scope -} - template double CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, elemT>& current_image_estimate) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h index cc4985bc80..265b702550 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.h @@ -87,7 +87,6 @@ class CudaGibbsPenalty : public GibbsPenalty 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; From eb9235c7347bedb714d8d54f2f5f1e0fdb1032cc Mon Sep 17 00:00:00 2001 From: Denis Prokopenko <22414094+denproc@users.noreply.github.com> Date: Thu, 12 Mar 2026 15:13:44 +0000 Subject: [PATCH 19/27] fixing else part for CudaGibbsPenalty --- src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 82b9f148bb..594d683005 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -753,8 +753,10 @@ CudaGibbsPenalty::set_kappa_sptr(const shared_ptr Date: Thu, 12 Mar 2026 15:14:07 +0000 Subject: [PATCH 20/27] adding cudaDeviceSynchronize --- src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index 594d683005..ce7a373d2b 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -490,7 +490,7 @@ CudaGibbsPenalty::compute_value(const DiscretisedDensity<3, e this->potential); checkCudaError("compute_value kernel"); - + cudaDeviceSynchronize(); double prior_value; prior_value = d_scalar[0]; @@ -622,7 +622,7 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity this->potential); checkCudaError("compute_hessian_diagonal kernel"); - + cudaDeviceSynchronize(); array_to_host(Hessian_diag, d_output_data); } @@ -666,7 +666,7 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD this->potential); checkCudaError("accumulate_Hessian_times_input kernel"); - + cudaDeviceSynchronize(); array_to_host(output, d_output_data); } From fde9800def51e8fdbe7685045aa5aff5c1dbc723 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Tue, 17 Mar 2026 14:12:30 +0000 Subject: [PATCH 21/27] Implementation of the GPU chunks and using CuVec start and end points in ParallelprojHelper --- .../ParallelprojHelper.h | 5 +- .../BackProjectorByBinParallelproj.cxx | 127 +++++------------- .../ForwardProjectorByBinParallelproj.cxx | 126 ++++------------- 3 files changed, 63 insertions(+), 195 deletions(-) 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/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index 025a245e04..0e49d128a7 100644 --- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx @@ -31,14 +31,10 @@ #include "stir/LORCoordinates.h" #include "stir/recon_array_functions.h" -//#ifdef parallelproj_built_with_CUDA -//# include "parallelproj_cuda.h" -//#else -//# include "parallelproj_c.h" -//#endif #include "parallelproj.h" -// for debugging, remove later + #include "stir/info.h" +#include "stir/format.h" #include "stir/error.h" #include "stir/stream.h" #include @@ -122,11 +118,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]; } @@ -160,9 +157,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; @@ -170,9 +165,6 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density long long num_lors_per_chunk; long long offset = 0; - // send image to all visible CUDA devices - float** image_on_cuda_devices = copy_float_array_to_all_devices(image_ptr, _helper->num_image_voxel); - // do (chuck-wise) back projection on the CUDA devices for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { @@ -186,101 +178,46 @@ 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, + image_ptr, + _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 ); 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, + image_ptr, + _helper->origin.data(), + _helper->voxsize.data(), + p.get_const_data_ptr() + offset, + num_lors_per_chunk, + _helper->imgdim.data()); } info("done", 2); offset += num_lors_per_chunk; } - // sum backprojected images on the first CUDA device - sum_float_arrays_on_first_device(image_on_cuda_devices, _helper->num_image_voxel); - - // copy summed image back to host - get_float_array_from_device(image_on_cuda_devices, _helper->num_image_voxel, 0, image_ptr); - - // 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_tof_sino_back(_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); - + cudaDeviceSynchronize(); p.release_const_data_ptr(); // --------------------------------------------------------------- // diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index b11219b2a5..bd2c0dd054 100644 --- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx @@ -25,19 +25,14 @@ #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 -//#if 0 -//#ifdef parallelproj_built_with_CUDA -//# include "parallelproj_cuda.h" -//#else -//# include "parallelproj_c.h" -//#endif + #include "parallelproj.h" -//#endif START_NAMESPACE_STIR @@ -88,14 +83,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) @@ -170,14 +157,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 @@ -187,10 +167,7 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& long long num_lors_per_chunk; long long offset = 0; - // send image to all visible CUDA devices - float** image_on_cuda_devices = copy_float_array_to_all_devices(image_ptr, _helper->num_image_voxel); - - // do (chuck-wise) projection on the CUDA devices + // do (chunk-wise) projection on the CUDA devices for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { if (chunk_num < remainder) @@ -206,22 +183,21 @@ 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, + image_ptr, + _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 ); float* STIR_mem = _projected_data_sptr->get_data_ptr(); @@ -230,70 +206,24 @@ 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, + image_ptr, + _helper->origin.data(), + _helper->voxsize.data(), + _projected_data_sptr->get_data_ptr() + offset, + num_lors_per_chunk, + _helper->imgdim.data()); if (chunk_num != _num_gpu_chunks - 1) _projected_data_sptr->release_data_ptr(); } offset += num_lors_per_chunk; } - // 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_tof_sino_fwd(_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); - + cudaDeviceSynchronize(); if (_density_sptr->is_contiguous()) { _density_sptr->release_full_data_ptr(); From b837cef47b4ae69c1f0916a17e85b40ad1a9bf04 Mon Sep 17 00:00:00 2001 From: Markus Jehl Date: Tue, 17 Mar 2026 15:13:16 +0000 Subject: [PATCH 22/27] Remove explicit sync statements. --- .../Parallelproj_projector/BackProjectorByBinParallelproj.cxx | 1 - .../Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx | 1 - 2 files changed, 2 deletions(-) diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index 0e49d128a7..7a3a417ea1 100644 --- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx @@ -217,7 +217,6 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density offset += num_lors_per_chunk; } - cudaDeviceSynchronize(); p.release_const_data_ptr(); // --------------------------------------------------------------- // diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index bd2c0dd054..2a0491de7f 100644 --- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx +++ b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx @@ -223,7 +223,6 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& offset += num_lors_per_chunk; } - cudaDeviceSynchronize(); if (_density_sptr->is_contiguous()) { _density_sptr->release_full_data_ptr(); From a321793cc0f168a0791a78f375deb2d3b90d3f06 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 18 Mar 2026 08:52:08 +0000 Subject: [PATCH 23/27] enabled parallelproj 1.* compatibility --- CMakeLists.txt | 7 ++- src/recon_buildblock/CMakeLists.txt | 13 ++++- .../BackProjectorByBinParallelproj.cxx | 55 +++++++++++++++++-- .../ForwardProjectorByBinParallelproj.cxx | 49 +++++++++++++++-- 4 files changed, 113 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 67c584dd0b..f6f6d3c724 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -261,12 +261,15 @@ endif() # Parallelproj if(NOT DISABLE_Parallelproj_PROJECTOR) - find_package(parallelproj 2.0 CONFIG) + find_package(parallelproj CONFIG) if (parallelproj_FOUND) set(STIR_WITH_Parallelproj_PROJECTOR ON) + 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) + 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() diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index b2dae8ba85..3ddde2de2b 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -167,7 +167,18 @@ endif() endif() if (STIR_WITH_Parallelproj_PROJECTOR) - target_link_libraries(${TARGET} PRIVATE parallelproj::parallelproj) + 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() if (HAVE_HDF5) diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx index 7a3a417ea1..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 @@ -31,7 +33,19 @@ #include "stir/LORCoordinates.h" #include "stir/recon_array_functions.h" -#include "parallelproj.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.h" +#endif #include "stir/info.h" #include "stir/format.h" @@ -165,7 +179,12 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density long long num_lors_per_chunk; long long offset = 0; - // do (chuck-wise) back projection on the CUDA devices +#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 (chunck-wise) back projection for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { if (chunk_num < remainder) @@ -185,7 +204,11 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density 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(), @@ -198,6 +221,10 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density _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(); @@ -206,17 +233,37 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density { 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()); + _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); + + // copy summed image back to host + get_float_array_from_device(image_on_cuda_devices, _helper->num_image_voxel, 0, image_ptr); + + // free image array from CUDA devices + free_float_array_on_all_devices(image_on_cuda_devices); +#endif + p.release_const_data_ptr(); // --------------------------------------------------------------- // diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx index 2a0491de7f..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 @@ -32,7 +35,18 @@ #include "stir/TOF_conversions.h" #include -#include "parallelproj.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.h" +#endif START_NAMESPACE_STIR @@ -167,7 +181,12 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& long long num_lors_per_chunk; long long offset = 0; - // do (chunk-wise) projection on the CUDA devices +#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 (chunk-wise) projection for (int chunk_num = 0; chunk_num < _num_gpu_chunks; chunk_num++) { if (chunk_num < remainder) @@ -185,7 +204,11 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& std::vector mem_for_PP(num_lors_per_chunk * _helper->num_tof_bins); 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 @@ -198,6 +221,10 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& _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(); @@ -211,18 +238,32 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>& { 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()); + _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); +#endif + if (_density_sptr->is_contiguous()) { _density_sptr->release_full_data_ptr(); From 05ad9de12924bcb2cf8a3c6383a0e825149f07ba Mon Sep 17 00:00:00 2001 From: Dimitra Kyriakopoulou Date: Wed, 18 Mar 2026 18:27:57 +0200 Subject: [PATCH 24/27] Expose num_gpu_chunks on ProjectorByBinPairUsingParallelproj --- .../ProjectorByBinPairUsingParallelproj.h | 4 ++++ .../ProjectorByBinPairUsingParallelproj.cxx | 19 +++++++++++++++++++ 2 files changed, 23 insertions(+) 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/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 From cbc8e4fbc6691b738f0154be27a62d931246e8f9 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 18 Mar 2026 21:50:25 +0000 Subject: [PATCH 25/27] add sync argument to array_to_host --- src/include/stir/cuda_utilities.h | 25 +++++++++++++++++-- .../CUDA/CudaGibbsPenalty.cuh | 3 --- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/include/stir/cuda_utilities.h b/src/include/stir/cuda_utilities.h index a8ff0eba20..be7f11e88e 100644 --- a/src/include/stir/cuda_utilities.h +++ b/src/include/stir/cuda_utilities.h @@ -54,6 +54,11 @@ 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) @@ -74,6 +79,10 @@ array_to_device(elemT* dev_data, const Array& stir_array) } } +//! copy an `Array` to pre-allocated CuVec +/*! + \ingroup CUDA +*/ template inline void array_to_device(CuVec& dev_data, const Array& stir_array) @@ -81,9 +90,14 @@ array_to_device(CuVec& dev_data, const Array& stir 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) +array_to_host(Array& stir_array, const elemT* dev_data, bool /* sync */ = true) { if (stir_array.is_contiguous()) { @@ -102,10 +116,17 @@ 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) +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()); } diff --git a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh index ce7a373d2b..cd760448aa 100644 --- a/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh +++ b/src/include/stir/recon_buildblock/CUDA/CudaGibbsPenalty.cuh @@ -534,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 @@ -622,7 +621,6 @@ CudaGibbsPenalty::compute_Hessian_diagonal(DiscretisedDensity this->potential); checkCudaError("compute_hessian_diagonal kernel"); - cudaDeviceSynchronize(); array_to_host(Hessian_diag, d_output_data); } @@ -666,7 +664,6 @@ CudaGibbsPenalty::accumulate_Hessian_times_input(DiscretisedD this->potential); checkCudaError("accumulate_Hessian_times_input kernel"); - cudaDeviceSynchronize(); array_to_host(output, d_output_data); } From 59dc91fdc56e67db2bd914c5bd083a892fd5fb45 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Wed, 18 Mar 2026 23:26:10 +0000 Subject: [PATCH 26/27] moved cuvec submodule to external_helpers/CuVec --- .gitmodules | 4 ++-- external_helpers/CMakeLists.txt | 8 ++++++++ cuvec => external_helpers/CuVec | 0 src/CMakeLists.txt | 3 --- 4 files changed, 10 insertions(+), 5 deletions(-) rename cuvec => external_helpers/CuVec (100%) diff --git a/.gitmodules b/.gitmodules index 3fb69c6576..2e6f41619f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,5 +2,5 @@ path = external_helpers/fmt url = https://github.com/fmtlib/fmt.git [submodule "cuvec"] - path = cuvec - url = https://github.com/AMYPAD/cuvec + path = external_helpers/CuVec + url = https://github.com/AMYPAD/CuVec 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/cuvec b/external_helpers/CuVec similarity index 100% rename from cuvec rename to external_helpers/CuVec diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dad8574fd5..aecf58f1e5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -262,9 +262,6 @@ if (STIR_ENABLE_EXPERIMENTAL) ADD_SUBDIRECTORY(${CMAKE_SOURCE_DIR}/src/experimental) endif() -set(CUVEC_INCLUDE_DIR "${CMAKE_SOURCE_DIR}/cuvec/cuvec/include") - -install(FILES "${CUVEC_INCLUDE_DIR}/cuvec.cuh" DESTINATION "${STIR_INCLUDE_INSTALL_DIR}/") # create an object library with all the "registries" # Warning: dependencies for object libraries are not transitive! add_library(stir_registries OBJECT ${STIR_REGISTRIES}) From ba91372e7cb57d87fb1b393394dcb55f41db4685 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Thu, 19 Mar 2026 17:29:32 +0000 Subject: [PATCH 27/27] updated release notes [ci skip] --- documentation/release_6.4.htm | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) 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: +
      +
    • + Compatibility with libparallelproj 2.0 and usage of CUDA managed pointers via + CuVec for internal variables in our parallelproj interface and CudaGibbsPenalty. + This results in a ~20% speed-up, but also code simplification. +
      + PR #1689. +
      • +
    • + Extra constructors for array, image and projdata classes that allow std::move for input arrays. +
      + PR #1693 and PR #1694. +
      • +
    +
    • +

    Changes to examples