Fix segfault from constexpr_math.h and adapt testmisc.cc to higher tolerances when running through valgrind

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/madgraph/iolibs/template_files/gpu/constexpr_math.h

@@ -10,13 +10,20 @@
 
 #include <cassert>
 #include <cmath>
+#include <iomanip>
 #include <limits>
 
-// FOR DEBUGGING!
+// SWITCH ON/OFF DEBUGGING
 #undef CONSTEXPR_MATH_DEBUG // no-debug
 //#define CONSTEXPR_MATH_DEBUG 1 // debug
+
+// FOR DEBUGGING
 #ifdef CONSTEXPR_MATH_DEBUG
-#define constexpr const
+#define CONSTEXPRMATHFUN inline
+#define CONSTEXPRMATHVAR const
+#else
+#define CONSTEXPRMATHFUN constexpr
+#define CONSTEXPRMATHVAR constexpr
 #endif
 
 // NB: namespaces mg5amcGpu and mg5amcCpu includes types which are defined in different ways for CPU and GPU builds (see #318 and #725)
@@ -38,6 +45,9 @@ namespace mg5amcCpu
       : std::numeric_limits<long double>::quiet_NaN();
   }
 
+  // SQRT constants
+  constexpr long double constexpr_sqrt2 = constexpr_sqrt( 2 );
+
   // Constexpr implementation of floor (see https://stackoverflow.com/a/66146159)
   constexpr int constexpr_floor( const long double xx )
   {
@@ -76,10 +86,13 @@ namespace mg5amcCpu
 
   // Constexpr implementation of sin for 0<x<pi/4 (long double signature)
   // Taylor expansion : x - x**3/3! + x**5/5!
-  constexpr long double sinTaylor( const long double xx )
+  CONSTEXPRMATHFUN long double sinTaylor( const long double xx )
   {
-    assert( xx >= 0 && "The argument of sinTaylor is assumed to be in [0,pi/4)" );
-    assert( xx < constexpr_pi_by_4 && "The argument of sinTaylor is assumed to be in [0,pi/4)" );
+#ifdef CONSTEXPR_MATH_DEBUG
+    if( xx < 0 || xx >= constexpr_pi_by_4 ) std::cout << "sinTaylor xx=" << xx << std::endl;
+#endif
+    assert( xx >= 0 && "The argument of sinTaylor is lower than the expected range [0,pi/4)" );
+    assert( xx < constexpr_pi_by_4 && "The argument of sinTaylor is higher than the expected range [0,pi/4)" );
     long double sinx = 0;
     int ipow = 1;
     long double delta = xx;
@@ -88,84 +101,170 @@ namespace mg5amcCpu
       long double sinxlast = sinx;
       sinx += delta;
 #ifdef CONSTEXPR_MATH_DEBUG
-      std::cout << "ipow=" << ipow << ", delta=" << delta << ", sinx=" << sinx << std::endl; // for debugging (not constexpr)
+      //std::cout << "ipow=" << ipow << ", delta=" << delta << ", sinx=" << sinx << std::endl; // for debugging (not constexpr)
 #endif
       if( sinx == sinxlast ) break;
       // Next iteration
       ipow += 2;
       delta *= -xx * xx / ( ipow - 1 ) / ipow;
     }
+#ifdef CONSTEXPR_MATH_DEBUG
+    //std::cout << "ipow=" << ipow << ", delta=" << delta << ", sinx=" << sinx << std::endl; // for debugging (not constexpr)
+#endif
     return sinx;
   }
 
+#ifdef CONSTEXPR_MATH_DEBUG
+  // Debug printouts for trig functions
+  inline void debug_constexpr_trig( const long double xx, size_t call )
+  {
+    CONSTEXPRMATHVAR long double xxminuspi = xx - constexpr_pi;
+    CONSTEXPRMATHVAR long double twopiminusxx = 2 * constexpr_pi - xx;
+    std::cout << std::setprecision( 40 );
+    std::cout << "constexpr_sin_quad call=" << call << std::endl
+              << "  xx=" << xx << std::endl
+              << "  (xx-pi)=" << xxminuspi << std::endl
+              << "  (2pi-xx)=" << twopiminusxx << std::endl;
+    std::cout << std::setprecision( 6 );
+    if( xx < 0 ) // (-inf, 0)
+      std::cout << "  -- case 1 (xx < 0)" << std::endl;
+    else if( xx == 0 ) // [0] *** NEW
+      std::cout << "  -- case 2 (xx == 0)" << std::endl;
+    else if( xx < constexpr_pi_by_4 ) // (0, 1/4*pi)
+      std::cout << "  -- case 3 (xx < pi/4)" << std::endl;
+    else if( xx == constexpr_pi_by_4 ) // [1/4*pi] *** NEW (3rd fix #903 assert fails)
+      std::cout << "  -- case 4 (xx == pi/4)" << std::endl;
+    else if( xx < constexpr_pi_by_2 ) // (1/4*pi, 1/2*pi)
+      std::cout << "  -- case 5 (xx < pi/2)" << std::endl;
+    else if( xx == constexpr_pi_by_2 ) // [1/2*pi] ** NEW
+      std::cout << "  -- case 6 (xx == pi/2)" << std::endl;
+    else if( xx < 3 * constexpr_pi_by_4 ) // (1/2*pi, 3/4*pi)
+      std::cout << "  -- case 7 (xx < 3*pi/4)" << std::endl;
+    else if( xx == 3 * constexpr_pi_by_4 ) // [3/4*pi] ** NEW
+      std::cout << "  -- case 8 (xx == 3*pi/4)" << std::endl;
+    else if( xx < constexpr_pi ) // (3/4*pi, pi)
+      std::cout << "  -- case 9 (xx < pi)" << std::endl;
+    else if( xx == constexpr_pi ) // [pi] *** NEW (2nd fix #903 infinite recursion)
+      std::cout << "  -- case 10 (xx == pi)" << std::endl;
+    else if( xx < 2 * constexpr_pi ) // (pi, 2*pi)
+      std::cout << "  -- case 11 (xx < 2*pi)" << std::endl;
+    else if( xx == 2 * constexpr_pi ) // [2*pi] ** NEW
+      std::cout << "  -- case 12 (xx == 2*pi)" << std::endl;
+    else // (2*pi, +inf)
+      std::cout << "  -- case 13 (xx > 2*pi)" << std::endl;
+  }
+#endif
+
   // Mapping to [0,2*pi) range (long double signature)
   constexpr long double mapIn0to2Pi( const long double xx )
   {
     return xx - constexpr_floor( xx / 2 / constexpr_pi ) * 2 * constexpr_pi;
   }
 
   // Constexpr implementation of cos (long double signature)
-  constexpr long double constexpr_cos_quad( const long double xx, const bool assume0to2Pi = false )
+  CONSTEXPRMATHFUN long double constexpr_cos_quad( const long double xx, const bool assume0to2Pi = false )
   {
     if( assume0to2Pi )
     {
       assert( xx >= 0 && "The argument of constexpr_cos_quad is assumed to be in [0,2*pi)" );
       assert( xx < 2 * constexpr_pi && "The argument of constexpr_cos_quad is assumed to be in [0,2*pi)" );
     }
-    if( xx < 0 )
+#ifdef CONSTEXPR_MATH_DEBUG
+    static size_t call = 0;
+    if( !assume0to2Pi )
+      call = 0;
+    else
+      call++;
+    if( call > 10 ) debug_constexpr_trig( xx, call );
+    assert( call < 20 );
+#endif
+    if( xx < 0 ) // (-inf, 0)
       return constexpr_cos_quad( mapIn0to2Pi( xx ), true );
-    else if( xx < constexpr_pi_by_4 ) // [0/4*pi, 1/4*pi)
+    else if( xx == 0 ) // [0] *** NEW
+      return 1;
+    else if( xx < constexpr_pi_by_4 ) // (0, 1/4*pi)
       return constexpr_sqrt( 1 - constexpr_pow( sinTaylor( xx ), 2 ) );
-    else if( xx < constexpr_pi_by_2 ) // [1/4*pi, 2/4*pi)
+    else if( xx == constexpr_pi_by_4 ) // [1/4*pi] *** NEW (3rd fix #903 assert fails)
+      return 1 / constexpr_sqrt2;
+    else if( xx < constexpr_pi_by_2 ) // (1/4*pi, 1/2*pi)
       return sinTaylor( constexpr_pi_by_2 - xx );
-    else if( xx < 3 * constexpr_pi_by_4 ) // [2/4*pi, 3/4*pi)
+    else if( xx == constexpr_pi_by_2 ) // [1/2*pi] ** NEW
+      return 0;
+    else if( xx < 3 * constexpr_pi_by_4 ) // (1/2*pi, 3/4*pi)
       return -sinTaylor( xx - constexpr_pi_by_2 );
-    else if( xx < constexpr_pi ) // [3/4*pi, 4/4*pi)
+    else if( xx == 3 * constexpr_pi_by_4 ) // [3/4*pi] ** NEW
+      return -1 / constexpr_sqrt2;
+    else if( xx < constexpr_pi ) // (3/4*pi, pi)
       return -constexpr_sqrt( 1 - constexpr_pow( sinTaylor( constexpr_pi - xx ), 2 ) );
-    else if( xx < 2 * constexpr_pi ) // [4/4*pi, 8/4*pi)
+    else if( xx == constexpr_pi ) // [pi] *** NEW (2nd fix #903 infinite recursion)
+      return -1;
+    else if( xx < 2 * constexpr_pi ) // (pi, 2*pi)
       return constexpr_cos_quad( 2 * constexpr_pi - xx, true );
-    else // [8/4*pi, +inf)
+    else if( xx == 2 * constexpr_pi ) // [2*pi] ** NEW
+      return 1;
+    else // (2*pi, +inf)
       return constexpr_cos_quad( mapIn0to2Pi( xx ), true );
   }
 
   // Constexpr implementation of cos (double signature, internally implemented as long double)
-  constexpr double constexpr_cos( const double x )
+  CONSTEXPRMATHFUN double constexpr_cos( const double x )
   {
     return constexpr_cos_quad( x );
   }
 
   // Constexpr implementation of sin (long double signature)
-  constexpr long double constexpr_sin_quad( const long double xx, const bool assume0to2Pi = false )
+  CONSTEXPRMATHFUN long double constexpr_sin_quad( const long double xx, const bool assume0to2Pi = false )
   {
     if( assume0to2Pi )
     {
       assert( xx >= 0 && "The argument of constexpr_sin_quad is assumed to be in [0,2*pi)" );
       assert( xx < 2 * constexpr_pi && "The argument of constexpr_sin_quad is assumed to be in [0,2*pi)" );
     }
-    if( xx < 0 )
+#ifdef CONSTEXPR_MATH_DEBUG
+    static size_t call = 0;
+    if( !assume0to2Pi )
+      call = 0;
+    else
+      call++;
+    if( call > 10 ) debug_constexpr_trig( xx, call );
+    assert( call < 20 );
+#endif
+    if( xx < 0 ) // (-inf, 0)
       return constexpr_sin_quad( mapIn0to2Pi( xx ), true );
-    else if( xx < constexpr_pi_by_4 ) // [0/4*pi, 1/4*pi)
+    else if( xx == 0 ) // [0] *** NEW
+      return 0;
+    else if( xx < constexpr_pi_by_4 ) // (0, 1/4*pi)
       return sinTaylor( xx );
-    else if( xx < constexpr_pi_by_2 ) // [1/4*pi, 2/4*pi)
+    else if( xx == constexpr_pi_by_4 ) // [1/4*pi] *** NEW (3rd fix #903 assert fails)
+      return 1 / constexpr_sqrt2;
+    else if( xx < constexpr_pi_by_2 ) // (1/4*pi, 1/2*pi)
       return constexpr_sqrt( 1 - constexpr_pow( sinTaylor( constexpr_pi_by_2 - xx ), 2 ) );
-    else if( xx < 3 * constexpr_pi_by_4 ) // [2/4*pi, 3/4*pi)
+    else if( xx == constexpr_pi_by_2 ) // [1/2*pi] ** NEW
+      return 1;
+    else if( xx < 3 * constexpr_pi_by_4 ) // (1/2*pi, 3/4*pi)
       return constexpr_sqrt( 1 - constexpr_pow( sinTaylor( xx - constexpr_pi_by_2 ), 2 ) );
-    else if( xx < constexpr_pi ) // [3/4*pi, 4/4*pi)
+    else if( xx == 3 * constexpr_pi_by_4 ) // [3/4*pi] ** NEW
+      return 1 / constexpr_sqrt2;
+    else if( xx < constexpr_pi ) // (3/4*pi, pi)
       return sinTaylor( constexpr_pi - xx );
-    else if( xx < 2 * constexpr_pi ) // [4/4*pi, 8/4*pi)
+    else if( xx == constexpr_pi ) // [pi] *** NEW (1st fix #903 infinite recursion)
+      return 0;
+    else if( xx < 2 * constexpr_pi ) // (pi, 2*pi)
       return -constexpr_sin_quad( 2 * constexpr_pi - xx, true );
-    else // [8/4*pi, +inf)
+    else if( xx == 2 * constexpr_pi ) // [2*pi] ** NEW
+      return 0;
+    else // (2*pi, +inf)
       return constexpr_sin_quad( mapIn0to2Pi( xx ), true );
   }
 
   // Constexpr implementation of sin (double signature, internally implemented as long double)
-  constexpr double constexpr_sin( const double x )
+  CONSTEXPRMATHFUN double constexpr_sin( const double x )
   {
     return constexpr_sin_quad( x );
   }
 
   // Constexpr implementation of tan (long double signature)
-  constexpr long double constexpr_tan_quad( const long double xx, const bool assume0to2Pi = false )
+  CONSTEXPRMATHFUN long double constexpr_tan_quad( const long double xx, const bool assume0to2Pi = false )
   {
     if( assume0to2Pi )
     {
@@ -181,14 +280,14 @@ namespace mg5amcCpu
   }
 
   // Constexpr implementation of tan (double signature, internally implemented as long double)
-  constexpr double constexpr_tan( const double x )
+  CONSTEXPRMATHFUN double constexpr_tan( const double x )
   {
     return constexpr_tan_quad( x );
   }
 
   // Constexpr implementation of atan for -1<x<1 (long double signature)
   // Taylor expansion : x - x**3/3 + x**5/5...
-  constexpr long double atanTaylor( const long double xx )
+  CONSTEXPRMATHFUN long double atanTaylor( const long double xx )
   {
     assert( xx >= -1 && "The argument of atanTaylor is assumed to be in (-1,+1)" );
     assert( xx < 1 && "The argument of atanTaylor is assumed to be in (-1,+1)" );
@@ -200,7 +299,7 @@ namespace mg5amcCpu
       long double atanxlast = atanx;
       atanx += xpow / ipow;
 #ifdef CONSTEXPR_MATH_DEBUG
-      std::cout << "ipow=" << ipow << ", xpow=" << xpow << ", atanx=" << atanx << std::endl; // for debugging (not constexpr)
+      //std::cout << "ipow=" << ipow << ", xpow=" << xpow << ", atanx=" << atanx << std::endl; // for debugging (not constexpr)
 #endif
       if( atanx == atanxlast ) break;
       // Next iteration
@@ -211,7 +310,7 @@ namespace mg5amcCpu
   }
 
   // Constexpr implementation of atan (long double signature)
-  constexpr long double constexpr_atan_quad( const long double xx )
+  CONSTEXPRMATHFUN long double constexpr_atan_quad( const long double xx )
   {
     if( xx > 1 )
       return constexpr_pi_by_2 - atanTaylor( 1 / xx );
@@ -226,7 +325,7 @@ namespace mg5amcCpu
   }
 
   // Constexpr implementation of atan (double signature, internally implemented as long double)
-  constexpr double constexpr_atan( const double x )
+  CONSTEXPRMATHFUN double constexpr_atan( const double x )
   {
     return constexpr_atan_quad( x );
   }

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/madgraph/iolibs/template_files/gpu/cudacpp.mk

@@ -381,6 +381,12 @@ ifeq ($(USEOPENMP),1)
   else ifneq ($(shell $(CXX) --version | egrep '^Intel'),)
     override OMPFLAGS = -fopenmp
     ###override OMPFLAGS = # disable OpenMP MT on Intel (was ok without GPUCC but not ok with GPUCC before #578)
+  else ifneq ($(shell $(CXX) --version | egrep '^clang version 16'),)
+    ###override OMPFLAGS = # disable OpenMP on clang16 #904
+    $(error OpenMP is not supported by cudacpp on clang16 - issue #904)
+  else ifneq ($(shell $(CXX) --version | egrep '^clang version 17'),)
+    ###override OMPFLAGS = # disable OpenMP on clang17 #904
+    $(error OpenMP is not supported by cudacpp on clang17 - issue #904)
   else ifneq ($(shell $(CXX) --version | egrep '^(clang)'),)
     override OMPFLAGS = -fopenmp
     ###override OMPFLAGS = # disable OpenMP MT on clang (was not ok without or with nvcc before #578)
@@ -643,11 +649,21 @@ endif
 
 # Target (and build options): debug
 MAKEDEBUG=
-debug: OPTFLAGS   = -g -O0
+debug: OPTFLAGS = -g -O0
 debug: CUDA_OPTFLAGS = -G
 debug: MAKEDEBUG := debug
 debug: all.$(TAG)
 
+# Target (and build options): address sanitizer #207
+###CXXLIBFLAGSASAN =
+###GPULIBFLAGSASAN =
+###asan: OPTFLAGS = -g -O0 -fsanitize=address -fno-omit-frame-pointer
+###asan: CUDA_OPTFLAGS = -G $(XCOMPILERFLAG) -fsanitize=address $(XCOMPILERFLAG) -fno-omit-frame-pointer
+###asan: CXXLIBFLAGSASAN = -fsanitize=address
+###asan: GPULIBFLAGSASAN = -Xlinker -fsanitize=address -Xlinker -shared
+###asan: MAKEDEBUG := debug
+###asan: all.$(TAG)
+
 # Target: tag-specific build lockfiles
 override oldtagsb=`if [ -d $(BUILDDIR) ]; then find $(BUILDDIR) -maxdepth 1 -name '.build.*' ! -name '.build.$(TAG)' -exec echo $(shell pwd)/{} \; ; fi`
 $(BUILDDIR)/.build.$(TAG):
@@ -768,11 +784,13 @@ endif
 #-------------------------------------------------------------------------------
 
 # Target (and build rules): C++ and CUDA/HIP standalone executables
+###$(cxx_checkmain): LIBFLAGS += $(CXXLIBFLAGSASAN)
 $(cxx_checkmain): LIBFLAGS += $(CXXLIBFLAGSRPATH) # avoid the need for LD_LIBRARY_PATH
 $(cxx_checkmain): $(BUILDDIR)/check_sa_cpp.o $(LIBDIR)/lib$(MG5AMC_CXXLIB).so $(cxx_objects_exe) $(BUILDDIR)/CurandRandomNumberKernel_cpp.o $(BUILDDIR)/HiprandRandomNumberKernel_cpp.o
 	$(CXX) -o $@ $(BUILDDIR)/check_sa_cpp.o $(OMPFLAGS) -ldl -pthread $(LIBFLAGS) -L$(LIBDIR) -l$(MG5AMC_CXXLIB) $(cxx_objects_exe) $(BUILDDIR)/CurandRandomNumberKernel_cpp.o $(BUILDDIR)/HiprandRandomNumberKernel_cpp.o $(RNDLIBFLAGS)
 
 ifneq ($(GPUCC),)
+###$(gpu_checkmain): LIBFLAGS += $(GPULIBFLAGSASAN)
 ifneq ($(shell $(CXX) --version | grep ^Intel),)
 $(gpu_checkmain): LIBFLAGS += -lintlc # compile with icpx and link with GPUCC (undefined reference to `_intel_fast_memcpy')
 $(gpu_checkmain): LIBFLAGS += -lsvml # compile with icpx and link with GPUCC (undefined reference to `__svml_cos4_l9')
@@ -787,9 +805,11 @@ endif
 #-------------------------------------------------------------------------------
 
 # Generic target and build rules: objects from Fortran compilation
+# (NB In this makefile, this only applies to fcheck_sa_fortran.o)
+# (NB -fPIC was added to fix clang16 build #904, but this seems better for other cases too and is consistent to c++ and cuda builds)
 $(BUILDDIR)/%%_fortran.o : %%.f *.inc
 	@if [ ! -d $(BUILDDIR) ]; then echo "mkdir -p $(BUILDDIR)"; mkdir -p $(BUILDDIR); fi
-	$(FC) -I. -c $< -o $@
+	$(FC) -I. -fPIC -c $< -o $@
 
 # Generic target and build rules: objects from Fortran compilation
 ###$(BUILDDIR)/%%_fortran.o : %%.f *.inc
@@ -800,6 +820,7 @@ $(BUILDDIR)/%%_fortran.o : %%.f *.inc
 # Target (and build rules): Fortran standalone executables
 ###$(BUILDDIR)/fcheck_sa_fortran.o : $(INCDIR)/fbridge.inc
 
+###$(cxx_fcheckmain): LIBFLAGS += $(CXXLIBFLAGSASAN)
 ifeq ($(UNAME_S),Darwin)
 $(cxx_fcheckmain): LIBFLAGS += -L$(shell dirname $(shell $(FC) --print-file-name libgfortran.dylib)) # add path to libgfortran on Mac #375
 endif
@@ -812,6 +833,7 @@ else
 endif
 
 ifneq ($(GPUCC),)
+###$(gpu_fcheckmain): LIBFLAGS += $(GPULIBFLAGSASAN)
 ifneq ($(shell $(CXX) --version | grep ^Intel),)
 $(gpu_fcheckmain): LIBFLAGS += -lintlc # compile with icpx and link with GPUCC (undefined reference to `_intel_fast_memcpy')
 $(gpu_fcheckmain): LIBFLAGS += -lsvml # compile with icpx and link with GPUCC (undefined reference to `__svml_cos4_l9')
@@ -914,10 +936,12 @@ endif
 ###endif
 
 ifeq ($(GPUCC),) # link only runTest_cpp.o
+###$(cxx_testmain): LIBFLAGS += $(CXXLIBFLAGSASAN)
 $(cxx_testmain): LIBFLAGS += $(CXXLIBFLAGSRPATH) # avoid the need for LD_LIBRARY_PATH
 $(cxx_testmain): $(LIBDIR)/lib$(MG5AMC_COMMONLIB).so $(cxx_objects_lib) $(cxx_objects_exe) $(GTESTLIBS)
 	$(CXX) -o $@ $(cxx_objects_lib) $(cxx_objects_exe) -ldl -pthread $(LIBFLAGS)
 else # link only runTest_$(GPUSUFFIX).o (new: in the past, this was linking both runTest_cpp.o and runTest_$(GPUSUFFIX).o)
+###$(gpu_testmain): LIBFLAGS += $(GPULIBFLAGSASAN)
 $(gpu_testmain): LIBFLAGS += $(GPULIBFLAGSRPATH) # avoid the need for LD_LIBRARY_PATH
 $(gpu_testmain): $(LIBDIR)/lib$(MG5AMC_COMMONLIB).so $(gpu_objects_lib) $(gpu_objects_exe) $(GTESTLIBS)
 ifneq ($(findstring hipcc,$(GPUCC)),) # link fortran/c++/hip using $FC when hipcc is used #802