File indexing completed on 2025-09-12 09:52:10
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 Vec2d nk = -k;
0103 std::cout << k << std::endl;
0104 std::cout << -k << std::endl;
0105 std::cout << nk << std::endl;
0106 std::cout << k + k << std::endl;
0107 std::cout << k * k << std::endl;
0108 Vec3d x(2.0, 4.0, 5.0);
0109 Vec3d y(-3.0, 2.0, -5.0);
0110 std::cout << x << std::endl;
0111 std::cout << y << std::endl;
0112
0113 Vec2d x2 = x.xy();
0114 Vec2d y2 = y.xy();
0115 std::cout << x2 << std::endl;
0116 Vec2d xx2 = x;
0117 std::cout << xx2 << std::endl;
0118 std::cout << y2 << std::endl;
0119
0120 std::cout << Vec2d(T(3.) * x2) << std::endl;
0121 std::cout << Vec2d(y2 * T(0.1)) << std::endl;
0122 std::cout << Vec2d(T(0.5) * (x2 + y2)) << std::endl;
0123 std::cout << mathSSE::sqrt(x2) << std::endl;
0124
0125 std::cout << dot(x2, y2) << " = 2?" << std::endl;
0126
0127 T z = cross(x2, y2);
0128 std::cout << z << " = 16?" << std::endl;
0129
0130 std::cout << mathSSE::sqrt(z) << " = 4?" << std::endl;
0131 }
0132
0133 template <typename T>
0134 void go(bool dovect = true) {
0135 typedef Vec4<T> Vec;
0136 typedef Vec2<T> Vec2D;
0137
0138 std::cout << std::endl;
0139 std::cout << sizeof(Vec) << ' ' << alignof(Vec) << std::endl;
0140 if (dovect) {
0141 std::vector<Vec> vec1;
0142 vec1.reserve(50);
0143 std::vector<T> vect(23);
0144 std::vector<Vec> vec2(53);
0145 std::vector<Vec> vec3;
0146 vec3.reserve(50234);
0147 }
0148
0149 Vec x(2.0, 4.0, 5.0);
0150 Vec y(-3.0, 2.0, -5.0);
0151 Vec nx = -x;
0152 std::cout << x << std::endl;
0153 std::cout << Vec4<float>(x) << std::endl;
0154 std::cout << Vec4<double>(x) << std::endl;
0155 std::cout << -x << std::endl;
0156 std::cout << nx << std::endl;
0157 std::cout << x.template get1<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 << (Vec(1) - y * T(0.1)) << std::endl;
0162 std::cout << mathSSE::sqrt(x) << std::endl;
0163
0164 std::cout << dot(x, y) << std::endl;
0165 std::cout << dot3(x, y) << std::endl;
0166 std::cout << dotSimple(x, y) << std::endl;
0167
0168 std::cout << "equal" << (x == x ? " " : " not ") << "ok" << std::endl;
0169 std::cout << "not equal" << (x == y ? " not " : " ") << "ok" << std::endl;
0170
0171 Vec z = cross(x, y);
0172 std::cout << z << std::endl;
0173
0174 std::cout << "rotations" << std::endl;
0175
0176 T a = 0.01;
0177 T ca = std::cos(a);
0178 T sa = std::sin(a);
0179
0180 Rot3<T> r1(ca, sa, 0, -sa, ca, 0, 0, 0, 1);
0181
0182 Rot2<T> r21(ca, sa, -sa, ca);
0183
0184 Rot3<T> r2(Vec(0, 1, 0), Vec(0, 0, 1), Vec(1, 0, 0));
0185 Rot2<T> r22(Vec2D(0, 1), Vec2D(1, 0));
0186
0187 {
0188 std::cout << "\n3D rot" << std::endl;
0189 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(x.xy());
0213 std::cout << x.xy() << 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 << x.xy() << 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 AVX2" << 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 }
0250
0251 #else
0252 int main() {
0253 typedef float T;
0254 typedef Vec4<T> Vec;
0255
0256 std::cout << std::endl;
0257 std::cout << sizeof(Vec) << ' ' << alignof(Vec) << std::endl;
0258 std::vector<Vec> vec1;
0259 vec1.reserve(50);
0260 std::vector<T> vect(23);
0261 std::vector<Vec> vec2(53);
0262 std::vector<Vec> vec3;
0263 vec3.reserve(50234);
0264
0265 return 0;
0266 }
0267 #endif
0268
0269 #else
0270 int main() { return 0; }
0271 #endif
