Merge pull request #1629 from ColomboMatte0/GPU_gibbs_priors

GPU and CPU gibbs priors refactoring

Author: ColomboMatte0

CMakeLists.txt

@@ -222,7 +222,9 @@ if(NOT DISABLE_STIR_CUDA)
       if("${CMAKE_VERSION}" VERSION_LESS "3.23")
         message(FATAL_ERROR "CMake 3.23 or newer is required to use CMAKE_CUDA_ARCHITECTURES set to 'all'. Upgrade your CMake or set DISABLE_STIR_CUDA to ON.")
       else()
-        set(CMAKE_CUDA_ARCHITECTURES "all")
+        if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
+          set(CMAKE_CUDA_ARCHITECTURES "all")
+        endif()
       endif()
       message(STATUS "STIR CUDA support enabled. Using CUDA version ${CUDAToolkit_VERSION}.")
       set(STIR_WITH_CUDA ON)

documentation/release_6.3.htm

@@ -46,6 +46,13 @@ <h4>General</h4>
       is disabled for backwards compatibility.<br>
       <a href=https://github.com/UCL/STIR/pull/1291>PR #1291</a>
     </li>
+    <li>
+      The priors code has been refactored to provide two common, parallelized base classes for CPU and GPU implementations (<code>GibbsPenalty</code> and <code>CudaGibbsPenalty</code>).<br>
+      On CPU, we introduce <code>GibbsQuadraticPenalty</code> and <code>GibbsRelativeDifferencePenalty</code>, which inherit from the parallelized <code>GibbsPenalty</code> base class (now parallelized using OpenMP).<br>
+      GPU implementations are also provided as <code>CudaGibbsQuadraticPenalty</code> and <code>CudaGibbsRelativeDifferencePenalty</code>.<br>
+      Currently, the new Gibbs penalties are still missing the paraboloidal surrogates related method implementations. They provide two new methods: <code>compute_gradient_times_input</code> and <code>compute_Hessian_diagonal</code>.<br>
+      <a href="https://github.com/UCL/STIR/pull/1629">PR #1629</a>
+    </li>
     <li>
       Data from GE Discovery MI systems in RDF9 should now be readable. TOF information on these scanners has also been
       added.
@@ -245,7 +252,7 @@ <h3>Bug fixes</h3>
     </li>
   </ul>
 
-  <h3>New deprecations</h3>
+  <h3>New Deprecations and renames</h3>
   <ul>
     <li>
       <code>truncate_end_planes</code> will be removed in v7.0
@@ -255,6 +262,11 @@ <h3>New deprecations</h3>
       CMake options <tt>STIR_PROJECTORS_AS_V3</tt>, <tt>STIR_LEGACY_IGNORE_VIEW_OFFSET</tt> and
       <tt>STIR_ROOT_ROTATION_AS_V4</tt> will therefore be removed.
     </li>
+    <li>
+      Starting from v7.0, <code>GeneralizedPrior</code> will be renamed to <code>GeneralizedPenalty</code>.<br>
+      <code>QuadraticPrior</code>, <code>RelativeDifferencePrior</code>, and <code>LogcoshPrior</code> will be removed in favor of the new <code>GibbsPenalty</code> framework.
+      <a href=https://github.com/UCL/STIR/issues/1426>Issue #1426</a>
+    </li>
   </ul>
 
   <h3>Test changes</h3>

recon_test_pack/OSMAPOSL_test_PMFromFile_QPweights.par

@@ -31,11 +31,11 @@ data_filename:=my_PMRT.pm
 End Projection Matrix By Bin From File Parameters:=
   End Projector Pair Using Matrix Parameters :=
 
-prior type := quadratic
-  Quadratic Prior Parameters:=
+prior type := Gibbs Quadratic
+  Gibbs Quadratic Penalty Parameters:=
   penalisation factor := .9
   weights:={{{0,1,0},{1,0,1},{0,1,0}}}
-  END Quadratic Prior Parameters:=
+  END Gibbs Quadratic Penalty Parameters:=
 
 
 end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=

recon_test_pack/OSMAPOSL_test_PM_QP.par

@@ -22,12 +22,12 @@ projector pair type := Matrix
   End Ray tracing matrix parameters :=
   End Projector Pair Using Matrix Parameters :=
 
-prior type := quadratic
-  Quadratic Prior Parameters:=
+prior type := Gibbs Quadratic
+  Gibbs Quadratic Penalty Parameters:=
   penalisation factor := 0.5
   ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D 
   only 2D:= 0
-  END Quadratic Prior Parameters:=
+  END Gibbs Quadratic Penalty Parameters:=
 
 end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
 

recon_test_pack/OSMAPOSL_test_PM_QPweights.par

@@ -20,11 +20,11 @@ projector pair type := Matrix
   End Ray tracing matrix parameters :=
   End Projector Pair Using Matrix Parameters :=
 
-prior type := quadratic
-  Quadratic Prior Parameters:=
+prior type := Gibbs Quadratic
+  Gibbs Quadratic Penalty Parameters:=
   penalisation factor := .9
   weights:={{{0,1,0},{1,0,1},{0,1,0}}}
-  END Quadratic Prior Parameters:=
+  END Gibbs Quadratic Penalty Parameters:=
 
 
 end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=

src/include/stir/cuda_utilities.h

@@ -1,5 +1,6 @@
 /*
     Copyright (C) 2024, University College London
+    Copyright (C) 2025, University of Milano-Bicocca
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -16,13 +17,40 @@
   \brief some utilities for STIR and CUDA
 
   \author Kris Thielemans
+  \author Matteo Neel Colombo
 */
 #include "stir/Array.h"
 #include "stir/info.h"
+#include "stir/error.h"
 #include <vector>
 
 START_NAMESPACE_STIR
 
+#ifndef __CUDACC__
+#  ifndef __host__
+#    define __host__
+#  endif
+#  ifndef __device__
+#    define __device__
+#  endif
+#endif
+
+#ifndef __CUDACC__
+struct cuda_dim3
+{
+  unsigned int x = 1, y = 1, z = 1;
+};
+struct cuda_int3
+{
+  int x = 0, y = 0, z = 0;
+};
+#else
+#  include <cuda_runtime.h>
+typedef dim3 cuda_dim3;
+typedef int3 cuda_int3;
+#endif
+
+#ifdef __CUDACC__
 template <int num_dimensions, typename elemT>
 inline void
 array_to_device(elemT* dev_data, const Array<num_dimensions, elemT>& stir_array)
@@ -64,6 +92,62 @@ array_to_host(Array<num_dimensions, elemT>& stir_array, const elemT* dev_data)
     }
 }
 
-END_NAMESPACE_STIR
+//! \brief Performs a parallel reduction sum on shared memory within a CUDA thread block, final value stored in shared_mem[0].
+template <typename elemT>
+__device__ inline void
+blockReduction(elemT* shared_mem, int thread_in_block, int block_threads)
+{
+  for (int stride = block_threads / 2; stride > 0; stride /= 2)
+    {
+      if (thread_in_block < stride)
+        shared_mem[thread_in_block] += shared_mem[thread_in_block + stride];
+      __syncthreads();
+    }
+}
+
+//! \brief Provides atomic addition for double values with fallback for pre-Pascal GPU architectures.
+template <typename elemT>
+__device__ inline double
+atomicAddGeneric(double* address, elemT val)
+{
+#  if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 600
+  return atomicAdd(address, static_cast<double>(val));
+#  else
+  if (threadIdx.x == 0 && threadIdx.y == 0 && threadIdx.z == 0 && blockIdx.x == 0 && blockIdx.y == 0 && blockIdx.z == 0)
+    {
+      printf("CudaGibbsPenalty: atomicAdd(double) unsupported on this GPU. "
+             "Upgrade to compute capability >= 6.0 or check code at "
+             "sources/STIR/src/include/stir/cuda_utilities.h:108.\n");
+      asm volatile("trap;");
+    }
+  return 0.0; // never reached
+              // Emulate atomicAdd for double precision on pre-Pascal architectures
+              // unsigned long long int* address_as_ull = reinterpret_cast<unsigned long long int*>(address);
+              // unsigned long long int old = *address_as_ull, assumed;
 
+  // do
+  //   {
+  //     assumed = old;
+  //     double updated = __longlong_as_double(assumed) + dval;
+  //     old = atomicCAS(address_as_ull, assumed, __double_as_longlong(updated));
+  // } while (assumed != old);
+
+  // return __longlong_as_double(old);
+#  endif
+}
+
+//! \brief Utility function to check for CUDA errors and report them with context information.
+inline void
+checkCudaError(const std::string& operation)
+{
+  cudaError_t cuda_error = cudaGetLastError();
+  if (cuda_error != cudaSuccess)
+    {
+      const char* err = cudaGetErrorString(cuda_error);
+      error(std::string("CudaGibbsPrior: CUDA error in ") + operation + ": " + err);
+    }
+}
+#endif
+
+END_NAMESPACE_STIR
 #endif