From 25069abe65cedcf47ffa8ec0731c1fac95cc4596 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 9 Oct 2025 12:29:55 +0200 Subject: [PATCH 1/3] Deprecate `vc` and `builtin_vc` build options --- README/ReleaseNotes/v638/index.md | 1 + cmake/modules/RootBuildOptions.cmake | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/README/ReleaseNotes/v638/index.md b/README/ReleaseNotes/v638/index.md index b61c524a7e5ed..d845c1e2c8c11 100644 --- a/README/ReleaseNotes/v638/index.md +++ b/README/ReleaseNotes/v638/index.md @@ -47,6 +47,7 @@ The following people have contributed to this new version: Relative RPATHs to the main ROOT libraries are unconditionally appended to all ROOT executables and libraries if the operating system supports it. If you want a ROOT build without RPATHs, use the canonical CMake variable `CMAKE_SKIP_INSTALL_RPATH=TRUE`. * The `TH1K` class is deprecated and will be removed in 6.40. It did not implement the `TH1` interface consistently, and limited the usability of the k-neighbors method it implemented by closely coupling the algorithm with the histogram class. Please use the new `TMath::KNNDensity` function that implements the same mathematical logic. +* The `vc` and `builtin_vc` build options are deprecated and will be removed in ROOT 6.40. They were not working anymore on newer platforms, and the vectorized math functions that they added to ROOT were not adopted by its users. ## Core Libraries * Behavior change: when selecting a template instantiation for a dictionary, all the template arguments have to be fully defined - the forward declarations are not enough any more. The error prompted by the dictionary generator will be `Warning: Unused class rule: MyTemplate`. diff --git a/cmake/modules/RootBuildOptions.cmake b/cmake/modules/RootBuildOptions.cmake index 9add71c9aa9fc..78d5695bdeaad 100644 --- a/cmake/modules/RootBuildOptions.cmake +++ b/cmake/modules/RootBuildOptions.cmake @@ -107,7 +107,7 @@ ROOT_BUILD_OPTION(builtin_pcre OFF "Build bundled copy of PCRE") ROOT_BUILD_OPTION(builtin_png OFF "Build bundled copy of libpng") ROOT_BUILD_OPTION(builtin_tbb OFF "Build TBB internally (requires network)") ROOT_BUILD_OPTION(builtin_unuran OFF "Build bundled copy of unuran") -ROOT_BUILD_OPTION(builtin_vc OFF "Build Vc internally (requires network)") +ROOT_BUILD_OPTION(builtin_vc OFF "Build Vc internally (requires network) [deprecated]") ROOT_BUILD_OPTION(builtin_vdt OFF "Build VDT internally (requires network)") ROOT_BUILD_OPTION(builtin_veccore OFF "Build VecCore internally (requires network)") ROOT_BUILD_OPTION(builtin_xrootd OFF "Build XRootD internally (requires network)") @@ -180,7 +180,7 @@ ROOT_BUILD_OPTION(unfold OFF "Enable the unfold package [GPL]") ROOT_BUILD_OPTION(unuran OFF "Enable support for UNURAN (package for generating non-uniform random numbers) [GPL]") ROOT_BUILD_OPTION(uring OFF "Enable support for io_uring (requires liburing and Linux kernel >= 5.1)") ROOT_BUILD_OPTION(use_gsl_cblas ON "Use the CBLAS library from GSL instead of finding a more optimized BLAS library automatically with FindBLAS (the GSL CBLAS is less performant but more portable)") -ROOT_BUILD_OPTION(vc OFF "Enable support for Vc (SIMD Vector Classes for C++)") +ROOT_BUILD_OPTION(vc OFF "Enable support for Vc (SIMD Vector Classes for C++) [deprecated]") ROOT_BUILD_OPTION(vdt ON "Enable support for VDT (fast and vectorisable mathematical functions)") ROOT_BUILD_OPTION(veccore OFF "Enable support for VecCore SIMD abstraction library") ROOT_BUILD_OPTION(vecgeom OFF "Enable support for VecGeom vectorized geometry library") @@ -402,7 +402,7 @@ foreach(opt afdsmgrd afs alien bonjour builtin_afterimage castor chirp cxx11 cxx endforeach() #---Deprecated options------------------------------------------------------------------------ -foreach(opt ) +foreach(opt builtin_vc vc) if(${opt}) message(DEPRECATION ">>> Option '${opt}' is deprecated and will be removed in the next release of ROOT. Please contact root-dev@cern.ch should you still need it.") endif() From fcb14ff4650125cf070e050ae3fd3daa6ebdc627 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 9 Oct 2025 12:30:15 +0200 Subject: [PATCH 2/3] [ci] Don't enabled deprecated build option `builtin_vc` on macOS 14 --- .github/workflows/root-ci-config/buildconfig/mac14.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/root-ci-config/buildconfig/mac14.txt b/.github/workflows/root-ci-config/buildconfig/mac14.txt index 56115d057e880..ea2b60a0ea6ec 100644 --- a/.github/workflows/root-ci-config/buildconfig/mac14.txt +++ b/.github/workflows/root-ci-config/buildconfig/mac14.txt @@ -18,7 +18,6 @@ builtin_openssl=ON builtin_pcre=ON builtin_tbb=ON builtin_unuran=ON -builtin_vc=ON builtin_vdt=ON builtin_veccore=ON builtin_xrootd=ON From 57c38f57d7453fbdf2dbd889210d753c2e6c225e Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 9 Oct 2025 12:30:57 +0200 Subject: [PATCH 3/3] [DO NOT MERGE] Remove code related to `vc` and `builtin_vc` This is to show how the removal would look like in the ROOT 6.40 development cycle. --- cmake/modules/SearchInstalledSoftware.cmake | 106 ----- config/RConfigure.in | 7 - core/clingutils/CMakeLists.txt | 3 - interpreter/cling/include/cling/vc.modulemap | 4 - math/mathcore/CMakeLists.txt | 10 - math/mathcore/inc/Math/Types.h | 23 -- math/mathcore/inc/VectorizedTMath.h | 33 -- math/mathcore/src/VectorizedTMath.cxx | 250 ------------ math/mathcore/test/CMakeLists.txt | 5 - math/mathcore/test/testVectorizedTMath.cxx | 106 ----- test/CMakeLists.txt | 14 - test/test-veccore.cxx | 3 - test/testGenVectorVc.cxx | 389 ------------------- test/testVc.cxx | 94 ----- 14 files changed, 1047 deletions(-) delete mode 100644 interpreter/cling/include/cling/vc.modulemap delete mode 100644 math/mathcore/inc/VectorizedTMath.h delete mode 100644 math/mathcore/src/VectorizedTMath.cxx delete mode 100644 math/mathcore/test/testVectorizedTMath.cxx delete mode 100644 test/testGenVectorVc.cxx delete mode 100644 test/testVc.cxx diff --git a/cmake/modules/SearchInstalledSoftware.cmake b/cmake/modules/SearchInstalledSoftware.cmake index 189f2fb70f1fd..ff838868529cf 100644 --- a/cmake/modules/SearchInstalledSoftware.cmake +++ b/cmake/modules/SearchInstalledSoftware.cmake @@ -1365,106 +1365,12 @@ if(builtin_tbb) set(TBB_TARGET TBB) endif() -#---Check for Vc--------------------------------------------------------------------- -if(builtin_vc) - unset(Vc_FOUND) - unset(Vc_FOUND CACHE) - set(vc ON CACHE BOOL "Enabled because builtin_vc requested (${vc_description})" FORCE) -elseif(vc) - if(fail-on-missing) - find_package(Vc 1.4.4 CONFIG QUIET REQUIRED) - else() - find_package(Vc 1.4.4 CONFIG QUIET) - if(NOT Vc_FOUND) - message(STATUS "Vc library not found, support for it disabled.") - message(STATUS "Please enable the option 'builtin_vc' to build Vc internally.") - set(vc OFF CACHE BOOL "Disabled because Vc not found (${vc_description})" FORCE) - endif() - endif() - if(Vc_FOUND) - set_property(DIRECTORY APPEND PROPERTY INCLUDE_DIRECTORIES ${Vc_INCLUDE_DIR}) - BUILD_ROOT_INCLUDE_PATH("${Vc_INCLUDE_DIR}") - endif() -endif() - -if(vc AND NOT Vc_FOUND) - ROOT_CHECK_CONNECTION_AND_DISABLE_OPTION("vc") -endif() - -if(vc AND NOT Vc_FOUND) - set(Vc_VERSION "1.4.4") - set(Vc_PROJECT "Vc-${Vc_VERSION}") - set(Vc_SRC_URI "${lcgpackages}/${Vc_PROJECT}.tar.gz") - set(Vc_DESTDIR "${CMAKE_BINARY_DIR}/externals") - set(Vc_ROOTDIR "${Vc_DESTDIR}/${CMAKE_INSTALL_PREFIX}") - set(Vc_LIBNAME "${CMAKE_STATIC_LIBRARY_PREFIX}Vc${CMAKE_STATIC_LIBRARY_SUFFIX}") - set(Vc_LIBRARY "${Vc_ROOTDIR}/lib/${Vc_LIBNAME}") - - ExternalProject_Add(VC - URL ${Vc_SRC_URI} - URL_HASH SHA256=5933108196be44c41613884cd56305df320263981fe6a49e648aebb3354d57f3 - BUILD_IN_SOURCE 0 - BUILD_BYPRODUCTS ${Vc_LIBRARY} - LOG_DOWNLOAD 1 LOG_CONFIGURE 1 LOG_BUILD 1 LOG_INSTALL 1 LOG_OUTPUT_ON_FAILURE 1 - CMAKE_ARGS -G ${CMAKE_GENERATOR} - -DCMAKE_POLICY_VERSION_MINIMUM=3.5 - -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} - -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} - -DCMAKE_C_FLAGS=${CMAKE_C_FLAGS} - -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} - -DCMAKE_CXX_FLAGS=${ROOT_EXTERNAL_CXX_FLAGS} - -DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX} - INSTALL_COMMAND env DESTDIR=${Vc_DESTDIR} ${CMAKE_COMMAND} --build . --target install - TIMEOUT 600 - ) - - set(VC_TARGET Vc) - set(Vc_LIBRARIES Vc) - set(Vc_INCLUDE_DIR ${Vc_ROOTDIR}/include) - set(Vc_CMAKE_MODULES_DIR ${Vc_ROOTDIR}/lib/cmake/Vc) - - add_library(VcExt STATIC IMPORTED) - set_property(TARGET VcExt PROPERTY IMPORTED_LOCATION ${Vc_LIBRARY}) - add_dependencies(VcExt VC) - - add_library(Vc INTERFACE) - target_include_directories(Vc SYSTEM BEFORE INTERFACE $) - target_link_libraries(Vc INTERFACE VcExt) - - find_package_handle_standard_args(Vc - FOUND_VAR Vc_FOUND - REQUIRED_VARS Vc_INCLUDE_DIR Vc_LIBRARIES Vc_CMAKE_MODULES_DIR - VERSION_VAR Vc_VERSION) - - # FIXME: This is a workaround to let ROOT find the headers at runtime if - # they are in the build directory. This is necessary until we decide how to - # treat externals with headers used by ROOT - if(NOT EXISTS ${CMAKE_BINARY_DIR}/include/Vc) - if (NOT EXISTS ${CMAKE_BINARY_DIR}/include) - execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_BINARY_DIR}/include) - endif() - execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink - ${Vc_INCLUDE_DIR}/Vc ${CMAKE_BINARY_DIR}/include/Vc) - endif() - # end of workaround - - install(DIRECTORY ${Vc_ROOTDIR}/ DESTINATION ".") -endif() - -if(Vc_FOUND) - # Missing from VcConfig.cmake - set(Vc_INCLUDE_DIRS ${Vc_INCLUDE_DIR}) -endif() - #---Check for VecCore-------------------------------------------------------------------- if(builtin_veccore) unset(VecCore_FOUND) unset(VecCore_FOUND CACHE) set(veccore ON CACHE BOOL "Enabled because builtin_veccore requested (${veccore_description})" FORCE) elseif(veccore) - if(vc) - set(VecCore_COMPONENTS Vc) - endif() if(fail-on-missing) find_package(VecCore 0.4.2 CONFIG QUIET REQUIRED COMPONENTS ${VecCore_COMPONENTS}) else() @@ -1521,18 +1427,6 @@ if(builtin_veccore) target_include_directories(VecCore SYSTEM INTERFACE $) add_dependencies(VecCore VECCORE) - if (Vc_FOUND) - set(VecCore_Vc_FOUND True) - set(VecCore_Vc_DEFINITIONS -DVECCORE_ENABLE_VC) - set(VecCore_Vc_INCLUDE_DIR ${Vc_INCLUDE_DIR}) - set(VecCore_Vc_LIBRARIES ${Vc_LIBRARIES}) - - set(VecCore_DEFINITIONS ${VecCore_Vc_DEFINITIONS}) - list(APPEND VecCore_INCLUDE_DIRS ${VecCore_Vc_INCLUDE_DIR}) - set(VecCore_LIBRARIES ${VecCore_LIBRARIES} ${Vc_LIBRARIES}) - target_link_libraries(VecCore INTERFACE ${Vc_LIBRARIES}) - endif() - find_package_handle_standard_args(VecCore FOUND_VAR VecCore_FOUND REQUIRED_VARS VecCore_INCLUDE_DIRS VecCore_LIBRARIES diff --git a/config/RConfigure.in b/config/RConfigure.in index a08b29d21fdd1..001d6e9d18ae3 100644 --- a/config/RConfigure.in +++ b/config/RConfigure.in @@ -36,7 +36,6 @@ #@haspthread@ R__HAS_PTHREAD /**/ #@hasxft@ R__HAS_XFT /**/ #@hascocoa@ R__HAS_COCOA /**/ -#@hasvc@ R__HAS_VC /**/ #@hasvdt@ R__HAS_VDT /**/ #@hasveccore@ R__HAS_VECCORE /**/ #@usecxxmodules@ R__USE_CXXMODULES /**/ @@ -54,12 +53,6 @@ #@hastbb@ R__HAS_TBB /**/ #define R__HARDWARE_INTERFERENCE_SIZE @hardwareinterferencesize@ /*Determined at CMake configure to be stable across all TUs*/ -#if defined(R__HAS_VECCORE) && defined(R__HAS_VC) -#ifndef VECCORE_ENABLE_VC -#define VECCORE_ENABLE_VC -#endif -#endif - #@uselz4@ R__HAS_DEFAULT_LZ4 /**/ #@usezlib@ R__HAS_DEFAULT_ZLIB /**/ #@uselzma@ R__HAS_DEFAULT_LZMA /**/ diff --git a/core/clingutils/CMakeLists.txt b/core/clingutils/CMakeLists.txt index 2a74a10eb0703..853aafb5b0865 100644 --- a/core/clingutils/CMakeLists.txt +++ b/core/clingutils/CMakeLists.txt @@ -117,9 +117,6 @@ set(clinginclude ${CMAKE_SOURCE_DIR}/interpreter/cling/include) set(custom_modulemaps) if (runtime_cxxmodules) set(custom_modulemaps boost.modulemap tinyxml2.modulemap cuda.modulemap module.modulemap.build) - if(vc) - set(custom_modulemaps ${custom_modulemaps} vc.modulemap) - endif() # We need to override the default modulemap because instead of producing a # single std.pcm, produces hundreds of pcms. This changed with MacOSX14.4.sdk diff --git a/interpreter/cling/include/cling/vc.modulemap b/interpreter/cling/include/cling/vc.modulemap deleted file mode 100644 index 20ee6cb240317..0000000000000 --- a/interpreter/cling/include/cling/vc.modulemap +++ /dev/null @@ -1,4 +0,0 @@ -module Vc { - header "Vc/Vc" - export * -} diff --git a/math/mathcore/CMakeLists.txt b/math/mathcore/CMakeLists.txt index f93c3557f9bd3..1c7892cd9bd9f 100644 --- a/math/mathcore/CMakeLists.txt +++ b/math/mathcore/CMakeLists.txt @@ -95,7 +95,6 @@ set(HEADERS TRandom3.h TRandomGen.h TStatistic.h - VectorizedTMath.h ) if(runtime_cxxmodules) @@ -104,14 +103,6 @@ if(runtime_cxxmodules) # 'abs' shadow declaration (from stl math.h) is broken. # FIXME: Revise after a llvm upgrade or reproduce it outside rootcling. list(REMOVE_ITEM HEADERS "Math/Math.h") - - if(vc) - # We do not link against libVc.a thus it makes no sense to check for - # version compatibility between libraries and header files. This fixes - # ROOT-11002 where upon building the modules.idx we run the static ctor - # runLibraryAbiCheck which fails to find the corresponding symbol. - set(dictoptions "-m" "Vc" "-mByproduct" "Vc" "-D" "Vc_NO_VERSION_CHECK") - endif(vc) endif() ROOT_ADD_C_FLAG(_flags -Wno-strict-overflow) # Avoid what it seems a compiler false positive warning @@ -189,7 +180,6 @@ ROOT_STANDARD_LIBRARY_PACKAGE(MathCore src/TStatistic.cxx src/UnBinData.cxx src/Util.cxx - src/VectorizedTMath.cxx LIBRARIES ${MATHCORE_LIBRARIES} DICTIONARY_OPTIONS diff --git a/math/mathcore/inc/Math/Types.h b/math/mathcore/inc/Math/Types.h index 95d71a78c7619..1a453031109b3 100644 --- a/math/mathcore/inc/Math/Types.h +++ b/math/mathcore/inc/Math/Types.h @@ -5,36 +5,13 @@ #ifdef R__HAS_VECCORE -#if defined(R__HAS_VC) - -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wall" -#pragma GCC diagnostic ignored "-Wunused-parameter" -#pragma GCC diagnostic ignored "-Wdeprecated-declarations" -#if (__cplusplus >= 202002L) // only for C++20 -#pragma GCC diagnostic ignored "-Wdeprecated-enum-enum-conversion" -#endif - -#ifdef __clang__ -#pragma clang diagnostic ignored "-Wconditional-uninitialized" -#pragma clang diagnostic ignored "-Wdeprecated-copy" -#endif - -#include -#pragma GCC diagnostic pop -#endif - #include namespace ROOT { namespace Internal { using ScalarBackend = vecCore::backend::Scalar; -#ifdef VECCORE_ENABLE_VC - using VectorBackend = vecCore::backend::VcVector; -#else using VectorBackend = vecCore::backend::Scalar; -#endif } using Float_v = typename Internal::VectorBackend::Float_v; using Double_v = typename Internal::VectorBackend::Double_v; diff --git a/math/mathcore/inc/VectorizedTMath.h b/math/mathcore/inc/VectorizedTMath.h deleted file mode 100644 index de95214c86194..0000000000000 --- a/math/mathcore/inc/VectorizedTMath.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef ROOT_VectorizedTMath -#define ROOT_VectorizedTMath - -#include "RtypesCore.h" -#include "Math/Types.h" -#include "TMath.h" - -#if defined(R__HAS_VECCORE) && defined(R__HAS_VC) - -namespace TMath { -::ROOT::Double_v Log2(::ROOT::Double_v &x); -::ROOT::Double_v BreitWigner(::ROOT::Double_v &x, Double_t mean = 0, Double_t gamma = 1); -::ROOT::Double_v Gaus(::ROOT::Double_v &x, Double_t mean = 0, Double_t sigma = 1, Bool_t norm = kFALSE); -::ROOT::Double_v LaplaceDist(::ROOT::Double_v &x, Double_t alpha = 0, Double_t beta = 1); -::ROOT::Double_v LaplaceDistI(::ROOT::Double_v &x, Double_t alpha = 0, Double_t beta = 1); -::ROOT::Double_v Freq(::ROOT::Double_v &x); -::ROOT::Double_v BesselI0_Split_More(::ROOT::Double_v &ax); -::ROOT::Double_v BesselI0_Split_Less(::ROOT::Double_v &x); -::ROOT::Double_v BesselI0(::ROOT::Double_v &x); -::ROOT::Double_v BesselI1_Split_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x); -::ROOT::Double_v BesselI1_Split_Less(::ROOT::Double_v &x); -::ROOT::Double_v BesselI1(::ROOT::Double_v &x); -::ROOT::Double_v BesselJ0_Split1_More(::ROOT::Double_v &ax); -::ROOT::Double_v BesselJ0_Split1_Less(::ROOT::Double_v &x); -::ROOT::Double_v BesselJ0(::ROOT::Double_v &x); -::ROOT::Double_v BesselJ1_Split1_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x); -::ROOT::Double_v BesselJ1_Split1_Less(::ROOT::Double_v &x); -::ROOT::Double_v BesselJ1(::ROOT::Double_v &x); -} // namespace TMath - -#endif // VECCORE and VC exist check - -#endif diff --git a/math/mathcore/src/VectorizedTMath.cxx b/math/mathcore/src/VectorizedTMath.cxx deleted file mode 100644 index bd3bc6b5fa8f4..0000000000000 --- a/math/mathcore/src/VectorizedTMath.cxx +++ /dev/null @@ -1,250 +0,0 @@ -#include "VectorizedTMath.h" - -#if defined(R__HAS_VECCORE) && defined(R__HAS_VC) - -namespace TMath { -//////////////////////////////////////////////////////////////////////////////// -::ROOT::Double_v Log2(::ROOT::Double_v &x) -{ - return vecCore::math::Log2(x); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Calculate a Breit Wigner function with mean and gamma. -::ROOT::Double_v BreitWigner(::ROOT::Double_v &x, Double_t mean, Double_t gamma) -{ - return 0.5 * M_1_PI * (gamma / (0.25 * gamma * gamma + (x - mean) * (x - mean))); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Calculate a gaussian function with mean and sigma. -/// If norm=kTRUE (default is kFALSE) the result is divided -/// by sqrt(2*Pi)*sigma. -::ROOT::Double_v Gaus(::ROOT::Double_v &x, Double_t mean, Double_t sigma, Bool_t norm) -{ - if (sigma == 0) - return ::ROOT::Double_v(1.e30); - - ::ROOT::Double_v inv_sigma = 1.0 / ::ROOT::Double_v(sigma); - ::ROOT::Double_v arg = (x - ::ROOT::Double_v(mean)) * inv_sigma; - - // For those entries of |arg| > 39 result is zero in double precision - ::ROOT::Double_v out = - vecCore::Blend<::ROOT::Double_v>(vecCore::math::Abs(arg) < ::ROOT::Double_v(39.0), - vecCore::math::Exp(::ROOT::Double_v(-0.5) * arg * arg), ::ROOT::Double_v(0.0)); - if (norm) - out *= 0.3989422804014327 * inv_sigma; // 1/sqrt(2*Pi)=0.3989422804014327 - return out; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Computes the probability density function of Laplace distribution -/// at point x, with location parameter alpha and shape parameter beta. -/// By default, alpha=0, beta=1 -/// This distribution is known under different names, most common is -/// double exponential distribution, but it also appears as -/// the two-tailed exponential or the bilateral exponential distribution -::ROOT::Double_v LaplaceDist(::ROOT::Double_v &x, Double_t alpha, Double_t beta) -{ - ::ROOT::Double_v beta_v_inv = ::ROOT::Double_v(1.0) / ::ROOT::Double_v(beta); - ::ROOT::Double_v out = vecCore::math::Exp(-vecCore::math::Abs((x - ::ROOT::Double_v(alpha)) * beta_v_inv)); - out *= ::ROOT::Double_v(0.5) * beta_v_inv; - return out; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Computes the distribution function of Laplace distribution -/// at point x, with location parameter alpha and shape parameter beta. -/// By default, alpha=0, beta=1 -/// This distribution is known under different names, most common is -/// double exponential distribution, but it also appears as -/// the two-tailed exponential or the bilateral exponential distribution -::ROOT::Double_v LaplaceDistI(::ROOT::Double_v &x, Double_t alpha, Double_t beta) -{ - ::ROOT::Double_v alpha_v = ::ROOT::Double_v(alpha); - ::ROOT::Double_v beta_v_inv = ::ROOT::Double_v(1.0) / ::ROOT::Double_v(beta); - return vecCore::Blend<::ROOT::Double_v>( - x <= alpha_v, 0.5 * vecCore::math::Exp(-vecCore::math::Abs((x - alpha_v) * beta_v_inv)), - 1 - 0.5 * vecCore::math::Exp(-vecCore::math::Abs((x - alpha_v) * beta_v_inv))); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Computation of the normal frequency function freq(x). -/// Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x. -/// -/// Translated from CERNLIB C300 by Rene Brun. -::ROOT::Double_v Freq(::ROOT::Double_v &x) -{ - Double_t c1 = 0.56418958354775629; - Double_t w2 = 1.41421356237309505; - - Double_t p10 = 2.4266795523053175e+2, q10 = 2.1505887586986120e+2, p11 = 2.1979261618294152e+1, - q11 = 9.1164905404514901e+1, p12 = 6.9963834886191355e+0, q12 = 1.5082797630407787e+1, - p13 = -3.5609843701815385e-2; - - Double_t p20 = 3.00459261020161601e+2, q20 = 3.00459260956983293e+2, p21 = 4.51918953711872942e+2, - q21 = 7.90950925327898027e+2, p22 = 3.39320816734343687e+2, q22 = 9.31354094850609621e+2, - p23 = 1.52989285046940404e+2, q23 = 6.38980264465631167e+2, p24 = 4.31622272220567353e+1, - q24 = 2.77585444743987643e+2, p25 = 7.21175825088309366e+0, q25 = 7.70001529352294730e+1, - p26 = 5.64195517478973971e-1, q26 = 1.27827273196294235e+1, p27 = -1.36864857382716707e-7; - - Double_t p30 = -2.99610707703542174e-3, q30 = 1.06209230528467918e-2, p31 = -4.94730910623250734e-2, - q31 = 1.91308926107829841e-1, p32 = -2.26956593539686930e-1, q32 = 1.05167510706793207e+0, - p33 = -2.78661308609647788e-1, q33 = 1.98733201817135256e+0, p34 = -2.23192459734184686e-2, q34 = 1; - - ::ROOT::Double_v v = vecCore::math::Abs(x) / w2; - - ::ROOT::Double_v result; - - vecCore::Mask<::ROOT::Double_v> mask1 = v < ::ROOT::Double_v(0.5); - vecCore::Mask<::ROOT::Double_v> mask2 = !mask1 && v < ::ROOT::Double_v(4.0); - vecCore::Mask<::ROOT::Double_v> mask3 = !(mask1 || mask2); - - ::ROOT::Double_v v2 = v * v; - ::ROOT::Double_v v3 = v2 * v; - ::ROOT::Double_v v4 = v3 * v; - ::ROOT::Double_v v5 = v4 * v; - ::ROOT::Double_v v6 = v5 * v; - ::ROOT::Double_v v7 = v6 * v; - ::ROOT::Double_v v8 = v7 * v; - - vecCore::MaskedAssign<::ROOT::Double_v>( - result, mask1, v * (p10 + p11 * v2 + p12 * v4 + p13 * v6) / (q10 + q11 * v2 + q12 * v4 + v6)); - vecCore::MaskedAssign<::ROOT::Double_v>( - result, mask2, - ::ROOT::Double_v(1.0) - - (p20 + p21 * v + p22 * v2 + p23 * v3 + p24 * v4 + p25 * v5 + p26 * v6 + p27 * v7) / - (vecCore::math::Exp(v2) * (q20 + q21 * v + q22 * v2 + q23 * v3 + q24 * v4 + q25 * v5 + q26 * v6 + v7))); - vecCore::MaskedAssign<::ROOT::Double_v>(result, mask3, - ::ROOT::Double_v(1.0) - - (c1 + (p30 * v8 + p31 * v6 + p32 * v4 + p33 * v2 + p34) / - ((q30 * v8 + q31 * v6 + q32 * v4 + q33 * v2 + q34) * v2)) / - (v * vecCore::math::Exp(v2))); - - return vecCore::Blend<::ROOT::Double_v>(x > 0, ::ROOT::Double_v(0.5) + ::ROOT::Double_v(0.5) * result, - ::ROOT::Double_v(0.5) * (::ROOT::Double_v(1) - result)); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Vectorized implementation of Bessel function I_0(x) for a vector x. -::ROOT::Double_v BesselI0_Split_More(::ROOT::Double_v &ax) -{ - ::ROOT::Double_v y = 3.75 / ax; - return (vecCore::math::Exp(ax) / vecCore::math::Sqrt(ax)) * - (0.39894228 + - y * (1.328592e-2 + - y * (2.25319e-3 + - y * (-1.57565e-3 + - y * (9.16281e-3 + - y * (-2.057706e-2 + y * (2.635537e-2 + y * (-1.647633e-2 + y * 3.92377e-3)))))))); -} - -::ROOT::Double_v BesselI0_Split_Less(::ROOT::Double_v &x) -{ - ::ROOT::Double_v y = x * x * 0.071111111; - - return 1.0 + - y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (3.60768e-2 + y * 4.5813e-3))))); -} - -::ROOT::Double_v BesselI0(::ROOT::Double_v &x) -{ - ::ROOT::Double_v ax = vecCore::math::Abs(x); - - return vecCore::Blend<::ROOT::Double_v>(ax <= 3.75, BesselI0_Split_Less(x), BesselI0_Split_More(ax)); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Vectorized implementation of modified Bessel function I_1(x) for a vector x. -::ROOT::Double_v BesselI1_Split_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x) -{ - ::ROOT::Double_v y = 3.75 / ax; - ::ROOT::Double_v result = - (vecCore::math::Exp(ax) / vecCore::math::Sqrt(ax)) * - (0.39894228 + y * (-3.988024e-2 + - y * (-3.62018e-3 + - y * (1.63801e-3 + y * (-1.031555e-2 + - y * (2.282967e-2 + y * (-2.895312e-2 + - y * (1.787654e-2 + y * -4.20059e-3)))))))); - return vecCore::Blend<::ROOT::Double_v>(x < 0, ::ROOT::Double_v(-1.0) * result, result); -} - -::ROOT::Double_v BesselI1_Split_Less(::ROOT::Double_v &x) -{ - ::ROOT::Double_v y = x * x * 0.071111111; - - return x * (0.5 + y * (0.87890594 + - y * (0.51498869 + y * (0.15084934 + y * (2.658733e-2 + y * (3.01532e-3 + y * 3.2411e-4)))))); -} - -::ROOT::Double_v BesselI1(::ROOT::Double_v &x) -{ - ::ROOT::Double_v ax = vecCore::math::Abs(x); - - return vecCore::Blend<::ROOT::Double_v>(ax <= 3.75, BesselI1_Split_Less(x), BesselI1_Split_More(ax, x)); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Vectorized implementation of Bessel function J0(x) for a vector x. -::ROOT::Double_v BesselJ0_Split1_More(::ROOT::Double_v &ax) -{ - ::ROOT::Double_v z = 8 / ax; - ::ROOT::Double_v y = z * z; - ::ROOT::Double_v xx = ax - 0.785398164; - ::ROOT::Double_v result1 = - 1 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6))); - ::ROOT::Double_v result2 = - -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7))); - return vecCore::math::Sqrt(0.636619772 / ax) * - (vecCore::math::Cos(xx) * result1 - z * vecCore::math::Sin(xx) * result2); -} - -::ROOT::Double_v BesselJ0_Split1_Less(::ROOT::Double_v &x) -{ - ::ROOT::Double_v y = x * x; - return (57568490574.0 + - y * (-13362590354.0 + y * (651619640.7 + y * (-11214424.18 + y * (77392.33017 + y * -184.9052456))))) / - (57568490411.0 + y * (1029532985.0 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y))))); -} - -::ROOT::Double_v BesselJ0(::ROOT::Double_v &x) -{ - ::ROOT::Double_v ax = vecCore::math::Abs(x); - return vecCore::Blend<::ROOT::Double_v>(ax < 8, BesselJ0_Split1_Less(x), BesselJ0_Split1_More(ax)); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Vectorized implementation of Bessel function J1(x) for a vector x. -::ROOT::Double_v BesselJ1_Split1_More(::ROOT::Double_v &ax, ::ROOT::Double_v &x) -{ - ::ROOT::Double_v z = 8 / ax; - ::ROOT::Double_v y = z * z; - ::ROOT::Double_v xx = ax - 2.356194491; - ::ROOT::Double_v result1 = - 1 + y * (0.183105e-2 + y * (-0.3516396496e-4 + y * (0.2457520174e-5 + y * -0.240337019e-6))); - ::ROOT::Double_v result2 = - 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5 + y * (-0.88228987e-6 - y * 0.105787412e-6))); - ::ROOT::Double_v result = - vecCore::math::Sqrt(0.636619772 / ax) * (vecCore::math::Cos(xx) * result1 - z * vecCore::math::Sin(xx) * result2); - vecCore::MaskedAssign<::ROOT::Double_v>(result, x < 0, -result); - return result; -} - -::ROOT::Double_v BesselJ1_Split1_Less(::ROOT::Double_v &x) -{ - ::ROOT::Double_v y = x * x; - return x * - (72362614232.0 + - y * (-7895059235.0 + y * (242396853.1 + y * (-2972611.439 + y * (15704.48260 + y * -30.16036606))))) / - (144725228442.0 + y * (2300535178.0 + y * (18583304.74 + y * (99447.43394 + y * (376.9991397 + y))))); -} - -::ROOT::Double_v BesselJ1(::ROOT::Double_v &x) -{ - ::ROOT::Double_v ax = vecCore::math::Abs(x); - return vecCore::Blend<::ROOT::Double_v>(ax < 8, BesselJ1_Split1_Less(x), BesselJ1_Split1_More(ax, x)); -} - -} // namespace TMath - -#endif diff --git a/math/mathcore/test/CMakeLists.txt b/math/mathcore/test/CMakeLists.txt index 94e72e706df49..8508a87c266fb 100644 --- a/math/mathcore/test/CMakeLists.txt +++ b/math/mathcore/test/CMakeLists.txt @@ -86,11 +86,6 @@ ROOT_ADD_GTEST(RanluxLCGUnit ranlux_lcg.cxx) ROOT_ADD_GTEST(RanluxppEngineTests RanluxppEngine.cxx LIBRARIES Core MathCore) -if(veccore AND vc) - ROOT_ADD_GTEST(VectorizedTMathUnit testVectorizedTMath.cxx - LIBRARIES Core MathCore) -endif() - ROOT_ADD_GTEST(testRootFinder testRootFinder.cxx LIBRARIES ${Libraries}) ROOT_ADD_GTEST(testKahan testKahan.cxx LIBRARIES Core MathCore) diff --git a/math/mathcore/test/testVectorizedTMath.cxx b/math/mathcore/test/testVectorizedTMath.cxx deleted file mode 100644 index c70591fa6e06e..0000000000000 --- a/math/mathcore/test/testVectorizedTMath.cxx +++ /dev/null @@ -1,106 +0,0 @@ -#include -#include -#include "VectorizedTMath.h" -#include - -#define N 16384 - -Double_t uniform_random(Double_t a, Double_t b) -{ - return a + (b - a) * drand48(); -} - -class VectorizedTMathTest : public ::testing::Test { -protected: - VectorizedTMathTest() {} - - size_t kVS = vecCore::VectorSize(); - Double_t input_array1[N] __attribute__((aligned(VECCORE_SIMD_ALIGN))); - Double_t input_array2[N] __attribute__((aligned(VECCORE_SIMD_ALIGN))); - - Double_t output_array[N] __attribute__((aligned(VECCORE_SIMD_ALIGN))); -}; - -#define TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, a, b) \ - TEST_F(VectorizedTMathTest, tmathfunc) \ - { \ - int trials = N; \ - for (int i = 0; i < trials; i++) \ - input_array1[i] = uniform_random(a, b); \ - \ - ROOT::Double_v x, y; \ - for (int j = 0; j < trials; j += kVS) { \ - vecCore::Load(x, &input_array1[j]); \ - y = TMath::tmathfunc(x); \ - vecCore::Store(y, &output_array[j]); \ - } \ - for (int j = 0; j < trials; j++) { \ - Double_t scalar_output = TMath::tmathfunc(input_array1[j]); \ - Double_t vec_output = output_array[j]; \ - Double_t re = \ - (scalar_output == vec_output && scalar_output == 0) ? 0 : (vec_output - scalar_output) / scalar_output; \ - EXPECT_NEAR(0, re, 1e9 * std::numeric_limits::epsilon()); \ - } \ - } - -#define TEST_VECTORIZED_TMATH_FUNCTION2(tmathfunc, a, b, c, d) \ - TEST_F(VectorizedTMathTest, tmathfunc) \ - { \ - int trials = N; \ - for (int i = 0; i < trials; i++) { \ - input_array1[i] = uniform_random(a, b); \ - input_array2[i] = uniform_random(c, d); \ - } \ - ROOT::Double_v x1, x2, y; \ - for (int j = 0; j < trials; j += kVS) { \ - vecCore::Load(x1, &input_array1[j]); \ - vecCore::Load(x2, &input_array2[j]); \ - y = TMath::tmathfunc(x1, x2); \ - vecCore::Store(y, &output_array[j]); \ - } \ - for (int j = 0; j < trials; j++) { \ - Double_t scalar_output = TMath::tmathfunc(input_array1[j], input_array2[j]); \ - Double_t vec_output = output_array[j]; \ - Double_t re = \ - (scalar_output == vec_output && scalar_output == 0) ? 0 : (vec_output - scalar_output) / scalar_output; \ - EXPECT_NEAR(0, re, 1e10 * std::numeric_limits::epsilon()); \ - } \ - } - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_RANGE(tmathfunc) TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, -FLT_MAX, FLT_MAX) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_POS_RANGE(tmathfunc) \ - TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, FLT_MIN, FLT_MAX) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_LOG(tmathfunc) TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, 1, FLT_MAX) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_EXP(tmathfunc) \ - TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, FLT_MIN_10_EXP, FLT_MAX_10_EXP) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_EXP_POS(tmathfunc) \ - TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, FLT_MIN, FLT_MAX_10_EXP) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_PROB(tmathfunc) TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, 0, 1) - -#define TEST_VECTORIZED_TMATH_FUNCTION_FLT_BESSEL(tmathfunc) TEST_VECTORIZED_TMATH_FUNCTION(tmathfunc, 0, 1e5) - -TEST_VECTORIZED_TMATH_FUNCTION_FLT_LOG(Log2); -TEST_VECTORIZED_TMATH_FUNCTION_FLT_RANGE(BreitWigner); -TEST_VECTORIZED_TMATH_FUNCTION_FLT_RANGE(Gaus); -TEST_VECTORIZED_TMATH_FUNCTION_FLT_EXP(LaplaceDist); -TEST_VECTORIZED_TMATH_FUNCTION_FLT_EXP(LaplaceDistI); - -// Freq function has exp(x^2). Hence range -6 to 6. -TEST_VECTORIZED_TMATH_FUNCTION(Freq, -6, 6); - -TEST_VECTORIZED_TMATH_FUNCTION_FLT_EXP(BesselI0); -TEST_VECTORIZED_TMATH_FUNCTION(BesselI1, FLT_MIN, FLT_MAX_10_EXP); - -TEST_VECTORIZED_TMATH_FUNCTION_FLT_BESSEL(BesselJ0); -TEST_VECTORIZED_TMATH_FUNCTION_FLT_BESSEL(BesselJ1); - -int main(int argc, char *argv[]) -{ - ::testing::InitGoogleTest(&argc, argv); - return RUN_ALL_TESTS(); -} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2c2242420599a..20174de89b787 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -332,13 +332,6 @@ if(MSVC) endif() ROOT_ADD_TEST(test-TFormulaTests COMMAND TFormulaTests FAILREGEX "FAILED|Error in") -#--Vc basic test----------------------------------------------------------------------------------- -if(ROOT_vc_FOUND) - ROOT_EXECUTABLE(testVc testVc.cxx LIBRARIES ${Vc_LIBRARIES} BUILTINS Vc) - target_include_directories(testVc SYSTEM BEFORE PRIVATE ${Vc_INCLUDE_DIRS}) - ROOT_ADD_TEST(test-Vc COMMAND testVc FAILREGEX "FAILED|Error in") -endif() - #--VecCore basic test------------------------------------------------------------------------------ if(ROOT_veccore_FOUND) ROOT_EXECUTABLE(test-veccore test-veccore.cxx LIBRARIES ${VecCore_LIBRARIES} BUILTINS VECCORE) @@ -351,13 +344,6 @@ if(ROOT_veccore_FOUND) endif() endif() -#--Vc GenVector test----------------------------------------------------------------------------------- -if(ROOT_vc_FOUND) - ROOT_EXECUTABLE(testGenVectorVc testGenVectorVc.cxx LIBRARIES Physics GenVector ${Vc_LIBRARIES} BUILTINS Vc) - target_include_directories(testGenVectorVc SYSTEM BEFORE PRIVATE ${Vc_INCLUDE_DIRS}) - ROOT_ADD_TEST(test-GenVector-Vc COMMAND testGenVectorVc FAILREGEX "FAILED|Error in") -endif() - #--g2root------------------------------------------------------------------------------------------ if(TARGET g2root) ROOT_ADD_TEST(test-g2root COMMAND g2root) diff --git a/test/test-veccore.cxx b/test/test-veccore.cxx index 6a8315b55880c..b7c9c1853af33 100644 --- a/test/test-veccore.cxx +++ b/test/test-veccore.cxx @@ -6,9 +6,6 @@ int main(int, char **) printf("VecCore available backends: Scalar "); #ifdef VECCORE_ENABLE_UMESIMD printf("UME::SIMD "); -#endif -#ifdef VECCORE_ENABLE_VC - printf("Vc"); #endif printf("\n"); return 0; diff --git a/test/testGenVectorVc.cxx b/test/testGenVectorVc.cxx deleted file mode 100644 index 4c09ee87bdc9a..0000000000000 --- a/test/testGenVectorVc.cxx +++ /dev/null @@ -1,389 +0,0 @@ -// ROOT -#include "Math/GenVector/PositionVector3D.h" -#include "Math/GenVector/DisplacementVector3D.h" -#include "Math/GenVector/Plane3D.h" -#include "Math/GenVector/Transform3D.h" -#include "TStopwatch.h" - -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wall" -#pragma GCC diagnostic ignored "-Wunused-parameter" -#include -#pragma GCC diagnostic pop - -// STL -#include -#include -#include -#include -#include -#include -#include - -template -T relativeError(const T &x, const T &y) -{ - if (x == y) - return 0; - - T diff = std::abs(x - y); - - if (x * y == T(0) || diff < std::numeric_limits::epsilon()) - return diff; - - return diff / (std::abs(x) + std::abs(y)); -} - -int compare(double x, double y, double tolerance = 1.0e-12) -{ - double error = relativeError(x, y); - - if (error > tolerance) { - int pr = std::cerr.precision(16); - std::cerr << "Error above tolerance:" << std::endl - << " expected = " << x << std::endl - << "true value = " << y << std::endl - << "abs. error = " << std::abs(x-y) << std::endl - << "rel. error = " << error << std::endl - << "tolerance = " << tolerance << std::endl; - std::cerr.precision(pr); - return 1; - } - - return 0; -} - -// randomn generator -static std::default_random_engine gen; -// Distributions for each member -static std::uniform_real_distribution p_x(-800, 800), p_y(-600, 600), p_z(10000, 10500); -static std::uniform_real_distribution d_x(-0.2, 0.2), d_y(-0.1, 0.1), d_z(0.95, 0.99); -static std::uniform_real_distribution c_x(3100, 3200), c_y(10, 15), c_z(3200, 3300); -static std::uniform_real_distribution r_rad(8500, 8600); -static std::uniform_real_distribution p0(-0.002, 0.002), p1(-0.2, 0.2), p2(0.97, 0.99), p3(-1300, 1300); - -template -class Data { -public: - typedef std::vector> Vector; - -public: - POINT position; - VECTOR direction; - POINT CoC; - PLANE plane; - FTYPE radius{0}; - -public: - template - Data(const INDATA &ind) - : position(ind.position.x(), ind.position.y(), ind.position.z()), - direction(ind.direction.x(), ind.direction.y(), ind.direction.z()), CoC(ind.CoC.x(), ind.CoC.y(), ind.CoC.z()), - plane(ind.plane.A(), ind.plane.B(), ind.plane.C(), ind.plane.D()), radius(ind.radius) - { - } - Data() - : position(p_x(gen), p_y(gen), p_z(gen)), direction(d_x(gen), d_y(gen), d_z(gen)), - CoC(c_x(gen), c_y(gen), c_z(gen)), plane(p0(gen), p1(gen), p2(gen), p3(gen)), radius(r_rad(gen)) - { - } -}; - -template -void fill(const INDATA &in, OUTDATA &out) -{ - out.clear(); - out.reserve(in.size()); - for (const auto &i : in) { - out.emplace_back(i); - } -} - -template ::value && - std::is_arithmetic::value && - std::is_arithmetic::value>::type> -inline bool reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, const FTYPE radius) -{ - constexpr FTYPE zero(0), two(2.0), four(4.0), half(0.5); - const FTYPE a = direction.Mag2(); - const VECTOR delta = position - CoC; - const FTYPE b = two * direction.Dot(delta); - const FTYPE c = delta.Mag2() - radius * radius; - const FTYPE discr = b * b - four * a * c; - const bool OK = discr > zero; - if (OK) { - const FTYPE dist = half * (std::sqrt(discr) - b) / a; - // change position to the intersection point - position += dist * direction; - // reflect the vector - // r = u - 2(u.n)n, r=reflection, u=incident, n=normal - const VECTOR normal = position - CoC; - direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal; - } - return OK; -} - -template ::value && - !std::is_arithmetic::value && - !std::is_arithmetic::value>::type> -inline typename FTYPE::mask_type reflectSpherical(POINT &position, VECTOR &direction, const POINT &CoC, - const FTYPE radius) -{ - const FTYPE two(2.0), four(4.0), half(0.5); - const FTYPE a = direction.Mag2(); - const VECTOR delta = position - CoC; - const FTYPE b = two * direction.Dot(delta); - const FTYPE c = delta.Mag2() - radius * radius; - FTYPE discr = b * b - four * a * c; - typename FTYPE::mask_type OK = discr > FTYPE::Zero(); - if (any_of(OK)) { - // Zero out the negative values in discr, to prevent sqrt(-ve) - discr(!OK) = FTYPE::Zero(); - // compute the distance - const FTYPE dist = half * (sqrt(discr) - b) / a; - // change position to the intersection point - position += dist * direction; - // reflect the vector - // r = u - 2(u.n)n, r=reflection, u=incident, n=normal - const VECTOR normal = position - CoC; - direction -= (two * normal.Dot(direction) / normal.Mag2()) * normal; - } - // return the mask indicating which results should be used - return OK; -} - -template ::value && - std::is_arithmetic::value>::type> -inline bool reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane) -{ - constexpr typename POINT::Scalar two(2.0); - const bool OK = true; - // Plane normal - const auto &normal = plane.Normal(); - // compute distance to the plane - const auto scalar = direction.Dot(normal); - const auto distance = -(plane.Distance(position)) / scalar; - // change position to reflection point and update direction - position += distance * direction; - direction -= two * scalar * normal; - return OK; -} - -template ::value && - !std::is_arithmetic::value>::type> -inline typename FTYPE::mask_type reflectPlane(POINT &position, VECTOR &direction, const PLANE &plane) -{ - const typename POINT::Scalar two(2.0); - const typename FTYPE::mask_type OK(true); - // Plane normal - const VECTOR normal = plane.Normal(); - // compute distance to the plane - const FTYPE scalar = direction.Dot(normal); - const FTYPE distance = -(plane.Distance(position)) / scalar; - // change position to reflection point and update direction - position += distance * direction; - direction -= two * scalar * normal; - return OK; -} - -template -using PositionVector = ROOT::Math::PositionVector3D, ROOT::Math::DefaultCoordinateSystemTag>; -template -using Vector = ROOT::Math::DisplacementVector3D, ROOT::Math::DefaultCoordinateSystemTag>; -template -using Plane = ROOT::Math::Impl::Plane3D; - -int main(int /*argc*/, char ** /*argv*/) -{ - int ret = 0; - - { - - const unsigned int nPhotons = 100; - std::cout << "Creating " << nPhotons << " random photons ..." << std::endl; - - // Scalar Types - Data, Vector, Plane, double>::Vector scalar_data(nPhotons); - - // Vc Types - Data, Vector, Plane, Vc::double_v>::Vector vc_data; - // Clone the exact random values from the Scalar vector - // Note we are making the same number of entries in the container, but each entry is a vector entry - // with Vc::double_t::Size entries. - fill(scalar_data, vc_data); - - // Loop over the two containers and compare - std::cout << "Ray Tracing :-" << std::endl; - - for (size_t i = 0; i < nPhotons; ++i) { - auto &sc = scalar_data[i]; - auto &vc = vc_data[i]; - - // ray tracing - reflectSpherical(sc.position, sc.direction, sc.CoC, sc.radius); - reflectPlane(sc.position, sc.direction, sc.plane); - reflectSpherical(vc.position, vc.direction, vc.CoC, vc.radius); - reflectPlane(vc.position, vc.direction, vc.plane); - - std::cout << "Position " << sc.position << " " << vc.position << std::endl; - std::cout << "Direction " << sc.direction << " " << vc.direction << std::endl; - - for (std::size_t j = 0; j < Vc::double_v::Size; ++j) { - ret |= compare(sc.position.x(), vc.position.x()[j]); - ret |= compare(sc.position.y(), vc.position.y()[j]); - ret |= compare(sc.position.z(), vc.position.z()[j]); - ret |= compare(sc.direction.x(), vc.direction.x()[j]); - ret |= compare(sc.direction.y(), vc.direction.y()[j]); - ret |= compare(sc.direction.z(), vc.direction.z()[j]); - } - } - - // Now test Transformation3D - std::cout << "Transforms :-" << std::endl; - for (size_t i = 0; i < nPhotons; ++i) { - auto &sc = scalar_data[i]; - auto &vc = vc_data[i]; - - // make 6 random scalar PositionVectors - PositionVector sp1(p_x(gen), p_y(gen), p_z(gen)); - PositionVector sp2(p_x(gen), p_y(gen), p_z(gen)); - PositionVector sp3(p_x(gen), p_y(gen), p_z(gen)); - PositionVector sp4(p_x(gen), p_y(gen), p_z(gen)); - PositionVector sp5(p_x(gen), p_y(gen), p_z(gen)); - PositionVector sp6(p_x(gen), p_y(gen), p_z(gen)); - // clone to Vc versions - PositionVector vp1(sp1.x(), sp1.y(), sp1.z()); - PositionVector vp2(sp2.x(), sp2.y(), sp2.z()); - PositionVector vp3(sp3.x(), sp3.y(), sp3.z()); - PositionVector vp4(sp4.x(), sp4.y(), sp4.z()); - PositionVector vp5(sp5.x(), sp5.y(), sp5.z()); - PositionVector vp6(sp6.x(), sp6.y(), sp6.z()); - - // Make transformations from points - // note warnings about axis not having the same angles expected here... - // point is to check scalar and vector versions do the same thing - const ROOT::Math::Impl::Transform3D st(sp1, sp2, sp3, sp4, sp5, sp6); - const ROOT::Math::Impl::Transform3D vt(vp1, vp2, vp3, vp4, vp5, vp6); - - // transform the vectors - const auto sv = st * sc.direction; - const auto vv = vt * vc.direction; - std::cout << "Transformed Direction " << sv << " " << vv << std::endl; - - // invert the transformations - const auto st_i = st.Inverse(); - const auto vt_i = vt.Inverse(); - - // Move the points back - const auto sv_i = st_i * sv; - const auto vv_i = vt_i * vv; - std::cout << "Transformed Back Direction " << sc.direction << " " << sv_i << " " << vv_i << std::endl; - - for (std::size_t j = 0; j < Vc::double_v::Size; ++j) { - ret |= compare(sv.x(), vv.x()[j]); - ret |= compare(sv.y(), vv.y()[j]); - ret |= compare(sv.z(), vv.z()[j]); - ret |= compare(sc.direction.x(), vv_i.x()[j]); - ret |= compare(sc.direction.y(), vv_i.y()[j]); - ret |= compare(sc.direction.z(), vv_i.z()[j]); - } - - ret |= compare(sc.direction.x(), sv_i.x()); - ret |= compare(sc.direction.y(), sv_i.y()); - ret |= compare(sc.direction.z(), sv_i.z()); - - // Make a scalar Plane - const double a(p0(gen)), b(p1(gen)), c(p2(gen)), d(p3(gen)); - Plane sc_plane(a, b, c, d); - // make a vector plane - Plane vc_plane(a, b, c, d); - - // transform the planes - const auto new_sc_plane = st * sc_plane; - const auto new_vc_plane = vt * vc_plane; - std::cout << "Transformed plane " << new_sc_plane << " " << new_vc_plane << std::endl; - - // now transform the planes back - const auto sc_plane_i = st_i * new_sc_plane; - const auto vc_plane_i = vt_i * new_vc_plane; - std::cout << "Transformed Back plane " << sc_plane_i << " " << vc_plane_i << std::endl; - - for (std::size_t j = 0; j < Vc::double_v::Size; ++j) { - ret |= compare(vc_plane.A()[j], vc_plane_i.A()[j]); - ret |= compare(vc_plane.B()[j], vc_plane_i.B()[j]); - ret |= compare(vc_plane.C()[j], vc_plane_i.C()[j]); - ret |= compare(vc_plane.D()[j], vc_plane_i.D()[j]); - ret |= compare(sc_plane_i.A(), vc_plane_i.A()[j]); - ret |= compare(sc_plane_i.B(), vc_plane_i.B()[j]); - ret |= compare(sc_plane_i.C(), vc_plane_i.C()[j]); - ret |= compare(sc_plane_i.D(), vc_plane_i.D()[j]); - } - } - } - - // now run some timing tests - { - const unsigned int nPhotons = 96000; // Must be multiple of 16 to avoid padding issues below... - - const unsigned int nTests = 1000; // number of tests to run - - // scalar data - Data, Vector, Plane, double>::Vector scalar_data(nPhotons); - // vector data with total equal number of photons (including vectorised size) - Data, Vector, Plane, Vc::double_v>::Vector vc_data( - nPhotons / Vc::double_v::Size); - - TStopwatch t; - - double best_time_scalar{9e30}, best_time_vector{9e30}; - - // time the scalar implementation - for (unsigned int i = 0; i < nTests; ++i) { - t.Start(); - for (auto &sc : scalar_data) { - reflectSpherical(sc.position, sc.direction, sc.CoC, sc.radius); - reflectPlane(sc.position, sc.direction, sc.plane); - } - t.Stop(); - const auto time = t.RealTime(); - if (time < best_time_scalar) { - best_time_scalar = time; - } - } - - // time the Vc implementation - for (unsigned int i = 0; i < nTests; ++i) { - t.Start(); - for (auto &vc : vc_data) { - reflectSpherical(vc.position, vc.direction, vc.CoC, vc.radius); - reflectPlane(vc.position, vc.direction, vc.plane); - } - t.Stop(); - const auto time = t.RealTime(); - if (time < best_time_vector) { - best_time_vector = time; - } - } - - std::cout << "Scalar best time = " << best_time_scalar << std::endl; - std::cout << "Vectorised Vc best time = " << best_time_vector << std::endl; - std::cout << "Vectorised Vc SIMD size = " << Vc::double_v::Size << std::endl; - std::cout << "Vectorised Vc speedup = " << best_time_scalar / best_time_vector << std::endl; - - // assert that the vector time is roughly Vc::double_v::Size times smaller than the scalar time - // allow 25% for 'safety' - // if (std::fabs((best_time_vector * Vc::double_v::Size) - best_time_scalar) > 0.25 * best_time_scalar) { - // ++ret; - // } - } - - if (ret) - std::cerr << "test FAILED !!! " << std::endl; - else - std::cout << "test OK " << std::endl; - return ret; -} diff --git a/test/testVc.cxx b/test/testVc.cxx deleted file mode 100644 index 3b8886a76408d..0000000000000 --- a/test/testVc.cxx +++ /dev/null @@ -1,94 +0,0 @@ -/* This file is part of the Vc project - Copyright (C) 2009-2010 Matthias Kretz - - Permission to use, copy, modify, and distribute this software - and its documentation for any purpose and without fee is hereby - granted, provided that the above copyright notice appear in all - copies and that both that the copyright notice and this - permission notice and warranty disclaimer appear in supporting - documentation, and that the name of the author not be used in - advertising or publicity pertaining to distribution of the - software without specific, written prior permission. - - The author disclaim all warranties with regard to this - software, including all implied warranties of merchantability - and fitness. In no event shall the author be liable for any - special, indirect or consequential damages or any damages - whatsoever resulting from loss of use, data or profits, whether - in an action of contract, negligence or other tortious action, - arising out of or in connection with the use or performance of - this software. - -*/ - -#ifdef __clang__ -#pragma clang diagnostic push -#pragma clang diagnostic ignored "-Wconditional-uninitialized" -#endif - -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wall" -#pragma GCC diagnostic ignored "-Wunused-parameter" -#include -#include -#pragma GCC diagnostic pop - -#ifdef __clang__ -#pragma clang diagnostic pop -#endif - -#include -#include - -template class Matrix; -template std::ostream &operator<<(std::ostream &, const Matrix &); - -template -class Matrix -{ - friend std::ostream &operator<< <>(std::ostream &, const Matrix &); - private: - typedef Vc::Vector V; - Vc::Memory m_mem; - public: - Matrix &operator=(const T &val) { - V vec(val); - for (unsigned int i = 0; i < m_mem.vectorsCount(); ++i) { - m_mem.vector(i) = vec; - } - return *this; - } - - Matrix &operator+=(const Matrix &rhs) { - for (unsigned int i = 0; i < m_mem.vectorsCount(); ++i) { - V v1(m_mem.vector(i)); - v1 += V(rhs.m_mem.vector(i)); - m_mem.vector(i) = v1; - } - return *this; - } -}; - -template -std::ostream &operator<<(std::ostream &out, const Matrix &m) -{ - for (unsigned int i = 0; i < Size; ++i) { - std::cout << "[" << std::setw(6) << m.m_mem[i * Size]; - for (unsigned int j = 1; j < Size; ++j) { - std::cout << std::setw(6) << m.m_mem[i * Size + j]; - } - std::cout << " ]\n"; - } - return out; -} - -int main() -{ - Matrix m1; - m1 = 1.f; - Matrix m2; - m2 = 2.f; - m1 += m2; - std::cout << m1 << std::endl; - return 0; -}