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
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
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
0264 int main() { return 0; }
0265 #endif
