forked from pwittich/mictest
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Matrix.h
120 lines (87 loc) · 3.46 KB
/
Matrix.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#ifndef _matrix_
#define _matrix_
#include "Math/SMatrix.h"
#include "Config.h"
typedef ROOT::Math::SMatrix<float,6,6,ROOT::Math::MatRepSym<float,6> > SMatrixSym66;
typedef ROOT::Math::SMatrix<float,6> SMatrix66;
typedef ROOT::Math::SVector<float,6> SVector6;
typedef ROOT::Math::SMatrix<float,3> SMatrix33;
typedef ROOT::Math::SMatrix<float,3,3,ROOT::Math::MatRepSym<float,3> > SMatrixSym33;
typedef ROOT::Math::SVector<float,3> SVector3;
typedef ROOT::Math::SMatrix<float,2> SMatrix22;
typedef ROOT::Math::SMatrix<float,2,2,ROOT::Math::MatRepSym<float,2> > SMatrixSym22;
typedef ROOT::Math::SVector<float,2> SVector2;
typedef ROOT::Math::SMatrix<float,3,6> SMatrix36;
typedef ROOT::Math::SMatrix<float,6,3> SMatrix63;
// should work with any SMatrix
template<typename Matrix>
void dumpMatrix(Matrix m)
{
for (int r=0;r<m.kRows;++r) {
for (int c=0;c<m.kCols;++c) {
std::cout << std::setw(12) << m.At(r,c) << " ";
}
std::cout << std::endl;
}
}
//==============================================================================
// This should go elsewhere, eventually.
#include <sys/time.h>
inline double dtime()
{
double tseconds = 0.0;
struct timeval mytime;
gettimeofday(&mytime,(struct timezone*)0);
tseconds = (double)(mytime.tv_sec + mytime.tv_usec*1.0e-6);
return( tseconds );
}
inline float hipo(float x, float y)
{
return sqrt(x*x + y*y);
}
inline void sincos4(float x, float& sin, float& cos)
{
// Had this writen with explicit division by factorial.
// The *whole* fitting test ran like 2.5% slower on MIC, sigh.
cos = 1;
sin = x; x *= x * 0.5f;
cos -= x; x *= x * 0.33333333f;
sin -= x; x *= x * 0.25f;
cos += x;
}
//==============================================================================
// This ifdef needs to be changed to something like "use matriplex" and/or
// "is icc" as we can only do vectorization with icc now.
#ifdef USE_MATRIPLEX
#ifdef __INTEL_COMPILER
#define ASSUME_ALIGNED(a, b) __assume_aligned(a, b)
#else
#define ASSUME_ALIGNED(a, b) __builtin_assume_aligned(a, b)
#endif
#include "Matriplex/MatriplexSym.h"
const Matriplex::idx_t NN = MPT_SIZE; // "Length" of MPlex.
const Matriplex::idx_t LL = 6; // Dimension of large/long MPlex entities
const Matriplex::idx_t HH = 3; // Dimension of small/short MPlex entities
typedef Matriplex::Matriplex<float, LL, LL, NN> MPlexLL;
typedef Matriplex::Matriplex<float, LL, 1, NN> MPlexLV;
typedef Matriplex::MatriplexSym<float, LL, NN> MPlexLS;
typedef Matriplex::Matriplex<float, HH, HH, NN> MPlexHH;
typedef Matriplex::Matriplex<float, HH, 1, NN> MPlexHV;
typedef Matriplex::MatriplexSym<float, HH, NN> MPlexHS;
typedef Matriplex::Matriplex<float, 2, 2, NN> MPlex22;
typedef Matriplex::Matriplex<float, 2, 1, NN> MPlex2V;
typedef Matriplex::MatriplexSym<float, 2, NN> MPlex2S;
typedef Matriplex::Matriplex<float, LL, HH, NN> MPlexLH;
typedef Matriplex::Matriplex<float, HH, LL, NN> MPlexHL;
typedef Matriplex::Matriplex<float, 1, 1, NN> MPlexQF;
typedef Matriplex::Matriplex<int, 1, 1, NN> MPlexQI;
#endif
//==============================================================================
#include <random>
extern std::default_random_engine g_gen;
extern std::normal_distribution<float> g_gaus;
extern std::uniform_real_distribution<float> g_unif;
// All debug printouts are ifdefed with DEBUG
// #define DEBUG
extern bool g_dump;
#endif