Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 10:50:38

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