|
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
|
#ifndef DataFormat_Math_ProjectMatrix_H
#define DataFormat_Math_ProjectMatrix_H
#define SMATRIX_USE_CONSTEXPR
#include "Math/SMatrix.h"
template <typename T, unsigned int N, unsigned int D>
struct ProjectMatrix {
typedef ROOT::Math::SMatrix<T, D, D, ROOT::Math::MatRepSym<T, D> > SMatDD;
typedef ROOT::Math::SMatrix<T, N, N> SMatNN;
typedef ROOT::Math::SMatrix<T, N, D> SMatND;
typedef ROOT::Math::SMatrix<T, D, N> SMatDN;
// no constructor
// H
SMatDN matrix() const {
SMatDN r;
for (unsigned int i = 0; i < D; i++)
r(i, index[i]) = T(1.);
return r;
}
// constuct given H
void fromH(SMatDN const& H) {
for (unsigned int i = 0; i < D; i++) {
index[i] = N;
for (unsigned int j = 0; j < N; j++)
if (H(i, j) > 0) {
index[i] = j;
break;
}
}
}
// H*S
SMatND project(SMatDD const& s) const {
SMatND r;
for (unsigned int i = 0; i < D; i++)
for (unsigned int j = 0; j < D; j++)
r(index[i], j) = s(i, j);
return r;
}
// K*H
SMatNN project(SMatND const& k) const {
SMatNN s;
for (unsigned int i = 0; i < N; i++)
for (unsigned int j = 0; j < D; j++)
s(i, index[j]) = k(i, j);
return s;
}
// S-K*H
void projectAndSubtractFrom(SMatNN& __restrict__ s, SMatND const& __restrict__ k) const {
for (unsigned int i = 0; i < N; i++)
for (unsigned int j = 0; j < D; j++)
s(i, index[j]) -= k(i, j);
}
// only H(i,index[i])=1.
unsigned int index[D];
};
#endif
|