Merge pull request #1449 from KrisThielemans/moreOpenMP

More openMP: find_min/max() and sum()

documentation/release_6.2.htm

@@ -100,15 +100,20 @@ <h3>Bug fixes</h3>
     but non-TOF instead. This did not happen in our normal reconstruction code, and would have thrown an error
     if it occured.
     <br>
-    Fixed in <a href=https://github.com/UCL/STIR/pull/1427>issue #1427</a>.
+    Fixed in <a href=https://github.com/UCL/STIR/pull/1427>PR #1427</a>.
   </li>
 </ul>
 
 
 <h3>Other code changes</h3>
 <ul>
   <li>
-    Fixes an incompatibility with C++20.
+    Fixed an incompatibility with C++20.
+  </li>
+  <li>
+    Enabled OpenMP for <code>Array</code> members <code>find_max(), find_min(), sum(), sum_positivie()</code>.
+    <br>
+    <a href=https://github.com/UCL/STIR/pull/1449>PR #1449</a>.
   </li>
 </ul>
 

src/CMakeLists.txt

@@ -136,7 +136,13 @@ if(STIR_OPENMP)
   find_package(OpenMP REQUIRED)  
   add_definitions(${OpenMP_CXX_FLAGS})
 
-  set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_CXX)
+  if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+    # work around https://gitlab.kitware.com/cmake/cmake/-/issues/26037
+    set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_CXX -latomic)
+    message(STATUS "OpenMP Linker flags for Clang: ${OpenMP_EXE_LINKER_FLAGS}")
+  else()
+    set (OpenMP_EXE_LINKER_FLAGS OpenMP::OpenMP_CXX)
+  endif()
 endif()
 
 #### Flags for compatibility between different systems

src/include/stir/Array.inl

@@ -240,6 +240,11 @@ Array<num_dimensions, elemT>::size_all() const
 {
   this->check_state();
   size_t acc = 0;
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(+ : acc)
+#  endif
+#endif
   for (int i = this->get_min_index(); i <= this->get_max_index(); i++)
     acc += this->num[i].size_all();
   return acc;
@@ -327,6 +332,11 @@ Array<num_dimensions, elemT>::sum() const
   this->check_state();
   typename HigherPrecision<elemT>::type acc;
   assign(acc, 0);
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(+ : acc)
+#  endif
+#endif
   for (int i = this->get_min_index(); i <= this->get_max_index(); i++)
     acc += this->num[i].sum();
   return static_cast<elemT>(acc);
@@ -339,6 +349,11 @@ Array<num_dimensions, elemT>::sum_positive() const
   this->check_state();
   typename HigherPrecision<elemT>::type acc;
   assign(acc, 0);
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(+ : acc)
+#  endif
+#endif
   for (int i = this->get_min_index(); i <= this->get_max_index(); i++)
     acc += this->num[i].sum_positive();
   return static_cast<elemT>(acc);
@@ -352,6 +367,11 @@ Array<num_dimensions, elemT>::find_max() const
   if (this->size() > 0)
     {
       elemT maxval = this->num[this->get_min_index()].find_max();
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(max : maxval)
+#  endif
+#endif
       for (int i = this->get_min_index() + 1; i <= this->get_max_index(); i++)
         {
           maxval = std::max(this->num[i].find_max(), maxval);
@@ -373,6 +393,11 @@ Array<num_dimensions, elemT>::find_min() const
   if (this->size() > 0)
     {
       elemT minval = this->num[this->get_min_index()].find_min();
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(min : minval)
+#  endif
+#endif
       for (int i = this->get_min_index() + 1; i <= this->get_max_index(); i++)
         {
           minval = std::min(this->num[i].find_min(), minval);
@@ -716,8 +741,13 @@ Array<1, elemT>::sum() const
   this->check_state();
   typename HigherPrecision<elemT>::type acc;
   assign(acc, 0);
-  for (int i = this->get_min_index(); i <= this->get_max_index(); acc += this->num[i++])
-    {}
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(+ : acc)
+#  endif
+#endif
+  for (int i = this->get_min_index(); i <= this->get_max_index(); ++i)
+    acc += this->num[i];
   return static_cast<elemT>(acc);
 };
 
@@ -728,6 +758,11 @@ Array<1, elemT>::sum_positive() const
   this->check_state();
   typename HigherPrecision<elemT>::type acc;
   assign(acc, 0);
+#ifdef STIR_OPENMP
+#  if _OPENMP >= 201107
+#    pragma omp parallel for reduction(+ : acc)
+#  endif
+#endif
   for (int i = this->get_min_index(); i <= this->get_max_index(); i++)
     {
       if (this->num[i] > 0)