Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-03-29 04:37:53

0001 #include "DataFormats/Math/interface/ExtVec.h"
0002 
0003 #include <cmath>
0004 #include <vector>
0005 
0006 #include <iostream>
0007 
0008 #ifdef __AVX2__
0009 #define CMS_USE_AVX2
0010 #endif /* __AVX2__ */
0011 
0012 void addScaleddiff(Vec3F& res, float s, Vec3F const& a, Vec3F const& b) { res = res + s * (a - b); }
0013 
0014 float dotV(Vec3F const& a, Vec3F const& b) { return dot(a, b); }
0015 
0016 float dotSimple(Vec3F const& a, Vec3F const& b) {
0017   Vec3F res = a * b;
0018   return res[0] + res[1] + res[2];
0019 }
0020 
0021 double dotSimple(Vec3D const& a, Vec3D const& b) {
0022   Vec3D res = a * b;
0023   return res[0] + res[1] + res[2];
0024 }
0025 
0026 float norm(Vec3F const& a) { return std::sqrt(dot(a, a)); }
0027 
0028 Vec3F toLocal(Vec3F const& a, Rot3<float> const& r) { return r.rotate(a); }
0029 
0030 Vec3F toGlobal(Vec3F const& a, Rot3<float> const& r) { return r.rotateBack(a); }
0031 
0032 // fake basicVector to check constructs...
0033 template <typename T>
0034 struct BaVec {
0035   typedef BaVec<T> self;
0036 
0037   BaVec() : theX(0), theY(0), theZ(0), theW(0) {}
0038 
0039   BaVec(float f1, float f2, float f3) : theX(f1), theY(f2), theZ(f3), theW(0) {}
0040 
0041   self& operator+=(self const& rh) { return *this; }
0042 
0043   T theX;
0044   T theY;
0045   T theZ;
0046   T theW;
0047 } __attribute__((aligned(16)));
0048 
0049 typedef BaVec<float> BaVecF;
0050 
0051 struct makeVec3F {
0052   makeVec3F(BaVecF& bv) : v(reinterpret_cast<Vec3F&>(bv)) {}
0053   Vec3F& v;
0054 };
0055 struct makeVec3FC {
0056   makeVec3FC(BaVecF const& bv) : v(reinterpret_cast<Vec3F const&>(bv)) {}
0057   Vec3F const& v;
0058 };
0059 
0060 template <>
0061 inline BaVecF& BaVecF::operator+=(BaVecF const& rh) {
0062   makeVec3FC v(rh);
0063   makeVec3F s(*this);
0064   s.v = s.v + v.v;
0065   return *this;
0066 }
0067 
0068 void sum(BaVecF& lh, BaVecF const& rh) { lh += rh; }
0069 
0070 void testBa() {
0071   std::cout << " test BA" << std::endl;
0072   BaVecF vx(2.0, 4.0, 5.0);
0073   BaVecF vy(-3.0, 2.0, -5.0);
0074   vx += vy;
0075   std::cout << vx.theX << ", " << vx.theY << ", " << vx.theZ << std::endl;
0076 }
0077 
0078 template <typename T>
0079 void go2d() {
0080   typedef Vec2<T> Vec2d;
0081   typedef Vec4<T> Vec3d;
0082 
0083   std::cout << "\n2d" << std::endl;
0084   std::cout << sizeof(Vec2d) << ' ' << alignof(Vec2d) << std::endl;
0085 
0086   std::vector<Vec2d> vec1;
0087   vec1.reserve(50);
0088   std::vector<T> vect(23);
0089   std::vector<Vec2d> vec2(53);
0090   std::vector<Vec2d> vec3;
0091   vec3.reserve(50234);
0092 
0093   constexpr Vec2d k{-2.0, 3.14};
0094   std::cout << k << std::endl;
0095   std::cout << k + k << std::endl;
0096   std::cout << k * k << std::endl;
0097   Vec3d x{2.0, 4.0, 5.0};
0098   Vec3d y{-3.0, 2.0, -5.0};
0099   std::cout << x << std::endl;
0100   std::cout << y << std::endl;
0101 
0102   Vec2d x2 = xy(x);
0103   Vec2d y2 = xy(y);
0104   std::cout << x2 << std::endl;
0105   Vec2d xx2 = xy(x);
0106   std::cout << xx2 << std::endl;
0107   std::cout << y2 << std::endl;
0108 
0109   std::cout << Vec2d(T(3.) * x2) << std::endl;
0110   std::cout << Vec2d(y2 * T(0.1)) << std::endl;
0111   std::cout << Vec2d(T(0.5) * (x2 + y2)) << std::endl;
0112   std::cout << apply(x2, [](T x) { return std::sqrt(x); }) << std::endl;
0113 
0114   std::cout << dot(x2, y2) << " = 2?" << std::endl;
0115   std::cout << dot2(x2, y) << " = 2?" << std::endl;
0116   std::cout << dot2(x, y) << " = 2?" << std::endl;
0117 
0118   T z = cross2(x2, y2);
0119   std::cout << z << " = 16?" << std::endl;
0120   z = cross2(x, y);
0121   std::cout << z << " = 16?" << std::endl;
0122 
0123   std::cout << std::sqrt(z) << " = 4?" << std::endl;
0124 }
0125 
0126 template <typename T>
0127 void go(bool dovec = true) {
0128   typedef Vec4<T> Vec;
0129   typedef Vec2<T> Vec2D;
0130 
0131   std::cout << std::endl;
0132   std::cout << sizeof(Vec) << ' ' << alignof(Vec) << std::endl;
0133 
0134   if (dovec) {
0135     std::vector<Vec> vec1;
0136     vec1.reserve(50);
0137     std::vector<T> vect(23);
0138     std::vector<Vec> vec2(53);
0139     std::vector<Vec> vec3;
0140     vec3.reserve(50234);
0141   }
0142 
0143   constexpr Vec zero{0, 0, 0, 0};
0144   constexpr Vec x{2.0f, 4.0f, 5.0f};
0145   constexpr Vec y{-3.0f, 2.0f, -5.0f};
0146   Vec x0 = x[0] + zero;
0147   //constexpr Vec xx = T(3.3) + zero;  // clang 3.8 does not like it
0148   const Vec xx = T(3.3) + zero;
0149   std::cout << x << std::endl;
0150   std::cout << (Vec4<float>){float(x[0]), float(x[1]), float(x[2]), float(x[3])} << std::endl;
0151   std::cout << (Vec4<double>){x[0], x[1], x[2], x[3]} << std::endl;
0152   std::cout << -x << std::endl;
0153   std::cout << Vec{x[2]} << std::endl;
0154   //std::cout <<  Vec(x[2]) << std::endl;
0155   std::cout << x0 << std::endl;
0156   std::cout << xx << std::endl;
0157   std::cout << Vec{} + x[2] << std::endl;
0158   std::cout << y << std::endl;
0159   std::cout << T(3.) * x << std::endl;
0160   std::cout << y * T(0.1) << std::endl;
0161   std::cout << (T(1) - y * T(0.1)) << std::endl;
0162   std::cout << apply(x, [](T x) { return std::sqrt(x); }) << std::endl;
0163 
0164   std::cout << dot(x, y) << std::endl;
0165   std::cout << dotSimple(x, y) << std::endl;
0166 
0167   //  std::cout << "equal" << (x==x ? " " : " not ") << "ok" << std::endl;
0168   // std::cout << "not equal" << (x==y ? " not " : " ") << "ok" << std::endl;
0169 
0170   Vec z = cross3(x, y);
0171   std::cout << z << std::endl;
0172   std::cout << cross2(x, y) << std::endl;
0173 
0174   std::cout << "\nrotations" << std::endl;
0175 
0176   // constexpr T a = 0.01;
0177   constexpr T ca = 0.9999500004166653;    // std::cos(a);  clang does not support constepxr math functions....
0178   constexpr T sa = 0.009999833334166664;  // std::sin(a);
0179 
0180   constexpr Rot3<T> r1(ca, sa, 0, -sa, ca, 0, 0, 0, 1);
0181 
0182   constexpr Rot2<T> r21(ca, sa, -sa, ca);
0183 
0184   constexpr Rot3<T> r2((Vec){0, 1, 0, 0}, (Vec){0, 0, 1, 0}, (Vec){1, 0, 0, 0});
0185   constexpr Rot2<T> r22((Vec2D){0, 1}, (Vec2D){1, 0});
0186 
0187   {
0188     std::cout << "\n3D rot" << std::endl;
0189     /* constexpr */ Vec xr = r1.rotate(x);
0190     std::cout << x << std::endl;
0191     std::cout << xr << std::endl;
0192     std::cout << r1.rotateBack(xr) << std::endl;
0193 
0194     Rot3<T> rt = r1.transpose();
0195     Vec xt = rt.rotate(xr);
0196     std::cout << x << std::endl;
0197     std::cout << xt << std::endl;
0198     std::cout << rt.rotateBack(xt) << std::endl;
0199 
0200     std::cout << r1 << std::endl;
0201     std::cout << rt << std::endl;
0202     std::cout << r1 * rt << std::endl;
0203     std::cout << r2 << std::endl;
0204     std::cout << r1 * r2 << std::endl;
0205     std::cout << r2 * r1 << std::endl;
0206     std::cout << r1 * r2.transpose() << std::endl;
0207     std::cout << r1.transpose() * r2 << std::endl;
0208   }
0209 
0210   {
0211     std::cout << "\n2D rot" << std::endl;
0212     Vec2D xr = r21.rotate(xy(x));
0213     std::cout << xy(x) << std::endl;
0214     std::cout << xr << std::endl;
0215     std::cout << r21.rotateBack(xr) << std::endl;
0216 
0217     Rot2<T> rt = r21.transpose();
0218     Vec2D xt = rt.rotate(xr);
0219     std::cout << xy(x) << std::endl;
0220     std::cout << xt << std::endl;
0221     std::cout << rt.rotateBack(xt) << std::endl;
0222 
0223     std::cout << r21 << std::endl;
0224     std::cout << rt << std::endl;
0225     std::cout << r21 * rt << std::endl;
0226     std::cout << r22 << std::endl;
0227     std::cout << r21 * r22 << std::endl;
0228     std::cout << r22 * r21 << std::endl;
0229     std::cout << r21 * r22.transpose() << std::endl;
0230     std::cout << r21.transpose() * r22 << std::endl;
0231   }
0232 }
0233 
0234 int main() {
0235 #ifdef CMS_USE_AVX2
0236   std::cout << "using AVX" << std::endl;
0237 #endif
0238   testBa();
0239   go<float>();
0240 #ifdef CMS_USE_AVX2
0241   go<double>(false);
0242 #else
0243   go<double>();
0244 #endif
0245   go2d<float>();
0246   go2d<double>();
0247 
0248   return 0;
0249 }