From f0481ff8c7bf0dc6a8dc630e891657bafbad3da7 Mon Sep 17 00:00:00 2001
From: Miha Muskinja <miha.muskinja@cern.ch>
Date: Thu, 16 Sep 2021 11:43:32 -0700
Subject: [PATCH] directly call the vector sincos function

---
 .../GeoSpecialShapes/LArWheelCalculator.h     |  7 --
 .../vec_parametrized_sincos.h                 | 95 ++++++++++---------
 .../src/LArWheelCalculator.cxx                |  3 -
 .../DistanceCalculatorSaggingOff.cxx          | 10 +-
 .../LArWheelCalculator_Impl/sincos_poly.cxx   | 10 +-
 5 files changed, 62 insertions(+), 63 deletions(-)

diff --git a/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/LArWheelCalculator.h b/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/LArWheelCalculator.h
index 0da3b1906837..345de3971416 100644
--- a/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/LArWheelCalculator.h
+++ b/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/LArWheelCalculator.h
@@ -126,9 +126,6 @@ class LArWheelCalculator
     std::pair<int, int> GetPhiGapAndSide(const CLHEP::Hep3Vector &p) const;
     double AmplitudeOfSurface(const CLHEP::Hep3Vector& P, int side, int fan_number) const;
 
-    /// Set the sincos calculator
-    // void SetVectorizedSincos() {m_sincos_calculator = &m_vsincos_par.eval;};
-
     /// @}
 
   private:
@@ -220,10 +217,6 @@ class LArWheelCalculator
     vsincos_par m_vsincos_par{};
 #endif
 
-  private:
-    typedef void (LArWheelCalculator::*SincosCalculator)(const double, double &, double &) const;
-    SincosCalculator m_sincos_calculator;
-    void EvalSincos(const double P, double &x, double &y) const {(this->*m_sincos_calculator)(P, x, y);};
 };
 
 #ifndef XAOD_STANDALONE
diff --git a/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/vec_parametrized_sincos.h b/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/vec_parametrized_sincos.h
index be996d7507b3..f662fbf17993 100644
--- a/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/vec_parametrized_sincos.h
+++ b/DetectorDescription/GeoModel/GeoSpecialShapes/GeoSpecialShapes/vec_parametrized_sincos.h
@@ -3,9 +3,10 @@
 */
 
 /**
- * @brief vectorized version of parametrized sincos
- * @ author Miha Muskinja, Chistos
- */
+  * @brief vectorized version of parametrized sincos
+  * see ATLASSIM-4753 for details
+  * @author Miha Muskinja, Chistos Anastopoulos
+*/
 #ifndef VEC_PARAMETRIZED_SINCOS_H
 #define VEC_PARAMETRIZED_SINCOS_H
 
@@ -14,61 +15,61 @@
 struct vsincos_par
 {
 
-  CxxUtils::vec<double, 4> param_0 = {};
-  CxxUtils::vec<double, 4> param_1 = {};
-  CxxUtils::vec<double, 4> param_2 = {};
+    CxxUtils::vec<double, 4> param_0 = {};
+    CxxUtils::vec<double, 4> param_1 = {};
+    CxxUtils::vec<double, 4> param_2 = {};
 
 #if HAVE_FUNCTION_MULTIVERSIONING
 #if defined(__x86_64__)
 
-  __attribute__((target("avx2,fma")))
+    __attribute__((target("avx2,fma")))
+
+    void
+    evaluate(const double r,
+             double &ATH_RESTRICT sin_a,
+             double &ATH_RESTRICT cos_a) const ATH_RESTRICT
+    {
+        const double r2 = r * r;
+        CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
+        P = r2 * P + param_2;
+        CxxUtils::vec<double, 4> P2 = {P[1], P[0], P[3], P[2]};
+        CxxUtils::vec<double, 4> res = r * P2 + P;
+        sin_a = res[0];
+        cos_a = res[2];
+    }
 
-  void
-  evaluate(const double r,
-           double& ATH_RESTRICT sin_a,
-           double& ATH_RESTRICT cos_a) const ATH_RESTRICT
-  {
-    const double r2 = r * r;
-    CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
-    P = r2 * P + param_2;
-    CxxUtils::vec<double, 4> P2 = { P[1], P[0], P[3], P[2] };
-    CxxUtils::vec<double, 4> res = r * P2 + P;
-    sin_a = res[0];
-    cos_a = res[2];
-  }
+    __attribute__((target("avx2")))
 
-  __attribute__((target("avx2")))
+    void
+    eval(const double r,
+         double &ATH_RESTRICT sin_a,
+         double &ATH_RESTRICT cos_a) const ATH_RESTRICT
+    {
+        const double r2 = r * r;
+        CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
+        P = r2 * P + param_2;
+        CxxUtils::vec<double, 4> P2 = {P[1], P[0], P[3], P[2]};
+        CxxUtils::vec<double, 4> res = r * P2 + P;
+        sin_a = res[0];
+        cos_a = res[2];
+    }
 
-  void
-  eval(const double r,
-       double& ATH_RESTRICT sin_a,
-       double& ATH_RESTRICT cos_a) const ATH_RESTRICT
-  {
-    const double r2 = r * r;
-    CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
-    P = r2 * P + param_2;
-    CxxUtils::vec<double, 4> P2 = { P[1], P[0], P[3], P[2] };
-    CxxUtils::vec<double, 4> res = r * P2 + P;
-    sin_a = res[0];
-    cos_a = res[2];
-  }
+    __attribute__((target("default")))
 
-  __attribute__((target("default")))
 #endif // x86
 #endif // FMV
 
-  void
-  eval(const double r,
-       double& ATH_RESTRICT sin_a,
-       double& ATH_RESTRICT cos_a) const ATH_RESTRICT
-  {
-    const double r2 = r * r;
-    CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
-    P = r2 * P + param_2;
-    sin_a = r * P[1] + P[0];
-    cos_a = r * P[3] + P[2];
-  }
+    void
+    eval(const double r,
+         double &ATH_RESTRICT sin_a,
+         double &ATH_RESTRICT cos_a) const ATH_RESTRICT
+    {
+        const double r2 = r * r;
+        CxxUtils::vec<double, 4> P = r2 * param_0 + param_1;
+        P = r2 * P + param_2;
+        sin_a = r * P[1] + P[0];
+        cos_a = r * P[3] + P[2];
+    }
 };
 
 #endif
-
diff --git a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator.cxx b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator.cxx
index cec1e276735f..11a481c860f3 100644
--- a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator.cxx
+++ b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator.cxx
@@ -443,9 +443,6 @@ LArWheelCalculator::LArWheelCalculator(LArG4::LArWheelCalculator_t a_wheelType,
      fclose(O);
      exit(0);
   */
-
-   // set the sincos calculator
-   m_sincos_calculator = &LArWheelCalculator::parameterized_sincos;
 }
 
 /* converts module gap number into wheel gap number */
diff --git a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/DistanceCalculatorSaggingOff.cxx b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/DistanceCalculatorSaggingOff.cxx
index 0163f1ba9442..06239cc8c1f4 100644
--- a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/DistanceCalculatorSaggingOff.cxx
+++ b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/DistanceCalculatorSaggingOff.cxx
@@ -48,7 +48,7 @@ namespace LArWheelCalculator_Impl
     const double cos_a = scalpha.cs, sin_a = scalpha.sn;
 #else // parameterized sine
     double cos_a, sin_a;
-    lwc()->EvalSincos(P.y(), sin_a, cos_a);
+    lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
 #endif
     // determination of the nearest quarter-wave number
     int nqwave = (z < 0.) ? 0 : int(z / lwc()->m_QuarterWaveLength);
@@ -160,7 +160,7 @@ namespace LArWheelCalculator_Impl
     double cos_a = scalpha.cs, sin_a = scalpha.sn;
 #else // parameterized sine
     double cos_a, sin_a;
-    lwc()->EvalSincos(P.y(), sin_a, cos_a);
+    lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
 #endif
 
     bool sqw = false;
@@ -251,7 +251,7 @@ namespace LArWheelCalculator_Impl
     const double cos_a = scalpha.cs, sin_a = scalpha.sn;
 #else // parameterized sine
     double cos_a, sin_a;
-    lwc()->EvalSincos(P.y(), sin_a, cos_a);
+    lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
 #endif
 
     int nqwave;
@@ -339,7 +339,7 @@ namespace LArWheelCalculator_Impl
     double cos_a = scalpha.cs, sin_a = scalpha.sn;
 #else // parameterized sine
     double cos_a, sin_a;
-    lwc()->EvalSincos(P.y(), sin_a, cos_a);
+    lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
 #endif
 
     bool sqw = false;
@@ -434,7 +434,7 @@ namespace LArWheelCalculator_Impl
     // parameterized sine
 #else
     double cos_a, sin_a;
-    lwc()->EvalSincos(P.y(), sin_a, cos_a);
+    lwc()->m_vsincos_par.eval(P.y(), sin_a, cos_a);
 #endif
 
     // determination of the nearest quarter-wave number
diff --git a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/sincos_poly.cxx b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/sincos_poly.cxx
index 917badd4cd08..deebabc0c293 100644
--- a/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/sincos_poly.cxx
+++ b/DetectorDescription/GeoModel/GeoSpecialShapes/src/LArWheelCalculator_Impl/sincos_poly.cxx
@@ -99,6 +99,8 @@ void LArWheelCalculator::fill_sincos_parameterization()
       m_cos_parametrization[i] = cos_parametrization[S][i];
     }
 
+    // Parameterization for the vectorized sincos calculation
+    // see ATLASSIM-4753 for details
     // s4, s5, c4, c5
     // s2, s3, c2, c3
     // s0, s1, c0, c1
@@ -155,6 +157,8 @@ void LArWheelCalculator::fill_sincos_parameterization()
     cos_parametrization[S][i] = params_cos[i];
   }
 
+  // Parameterization for the vectorized sincos calculation
+  // see ATLASSIM-4753 for details
   // s4, s5, c4, c5
   // s2, s3, c2, c3
   // s0, s1, c0, c1
@@ -191,8 +195,12 @@ void LArWheelCalculator::fill_sincos_parameterization()
     parameterized_sincos(r, sin_a, cos_a);
     m_vsincos_par.eval(r, sin_a_v, cos_a_v);
 #if DEBUGPRINT
-    std::cout << "default: " << r << " " << sin_a << " " << cos_a << std::endl;
+    std::streamsize ss = std::cout.precision();
+    std::cout.precision(16);
+    std::cout << "def: " << r << " " << sin_a << " " << cos_a << std::endl;
     std::cout << "vec: " << r << " " << sin_a_v << " " << cos_a_v << std::endl;
+    std::cout << "dif: " << r << " " << (sin_a - sin_a_v) / sin_a << " " << (cos_a - cos_a_v) / cos_a << std::endl;
+    std::cout.precision(ss);
 #endif
     double ds = fabs(scalpha.sn - sin_a);
     if(ds > dsin){
-- 
GitLab