File indexing completed on 2024-10-04 22:54:37
0001 #ifdef NEWSMATRIX
0002
0003 #define SMATRIX_USE_CONSTEXPR
0004 #include "Math/SMatrix.h"
0005
0006 typedef unsigned int IndexType;
0007
0008
0009 namespace ROOT {
0010
0011 namespace Math {
0012
0013 constexpr int fOff15x15[] = {
0014 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 1, 2, 4, 7, 11, 16,
0015 22, 29, 37, 46, 56, 67, 79, 92, 106, 3, 4, 5, 8, 12, 17, 23, 30, 38, 47, 57, 68,
0016 80, 93, 107, 6, 7, 8, 9, 13, 18, 24, 31, 39, 48, 58, 69, 81, 94, 108, 10, 11, 12,
0017 13, 14, 19, 25, 32, 40, 49, 59, 70, 82, 95, 109, 15, 16, 17, 18, 19, 20, 26, 33, 41,
0018 50, 60, 71, 83, 96, 110, 21, 22, 23, 24, 25, 26, 27, 34, 42, 51, 61, 72, 84, 97, 111,
0019 28, 29, 30, 31, 32, 33, 34, 35, 43, 52, 62, 73, 85, 98, 112, 36, 37, 38, 39, 40, 41,
0020 42, 43, 44, 53, 63, 74, 86, 99, 113, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 64, 75,
0021 87, 100, 114, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 76, 88, 101, 115, 66, 67, 68,
0022 69, 70, 71, 72, 73, 74, 75, 76, 77, 89, 102, 116, 78, 79, 80, 81, 82, 83, 84, 85, 86,
0023 87, 88, 89, 90, 103, 117, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 118,
0024 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119};
0025
0026 constexpr int const* pfOff15x15[] = {fOff15x15 + 0,
0027 fOff15x15 + 15,
0028 fOff15x15 + 2 * 15,
0029 fOff15x15 + 3 * 15,
0030 fOff15x15 + 4 * 15,
0031 fOff15x15 + 5 * 15,
0032 fOff15x15 + 6 * 15,
0033 fOff15x15 + 7 * 15,
0034 fOff15x15 + 8 * 15,
0035 fOff15x15 + 9 * 15,
0036 fOff15x15 + 10 * 15,
0037 fOff15x15 + 11 * 15,
0038 fOff15x15 + 12 * 15,
0039 fOff15x15 + 13 * 15,
0040 fOff15x15 + 14 * 15};
0041 template <>
0042 struct RowOffsets<15> {
0043 RowOffsets() {}
0044 int operator()(unsigned int i, unsigned int j) const { return pfOff15x15[i][j]; }
0045 int apply(unsigned int i) const { return fOff15x15[i]; }
0046 };
0047
0048
0049
0050
0051
0052
0053
0054
0055 struct AssignAsSym {
0056
0057 template <class T, IndexType D, class A, class R>
0058 static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const Expr<A, T, D, D, R>& rhs) {
0059
0060 for (IndexType i = 0; i < D; ++i)
0061
0062 for (IndexType j = 0; j <= i; ++j) {
0063 lhs(i, j) = rhs(i, j);
0064 }
0065 }
0066
0067 template <class T, IndexType D, class R>
0068 static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const SMatrix<T, D, D, R>& rhs) {
0069
0070 for (IndexType i = 0; i < D; ++i)
0071
0072 for (IndexType j = 0; j <= i; ++j) {
0073 lhs(i, j) = rhs(i, j);
0074 }
0075 }
0076 };
0077
0078 template <class T, IndexType D1, IndexType D2, class R>
0079 inline SMatrix<T, D1, D1, MatRepSym<T, D1> > Similarity(const SMatrix<T, D1, D2, R>& lhs,
0080 const SMatrix<T, D2, D2, MatRepStd<T, D2> >& rhs) {
0081 SMatrix<T, D1, D2, MatRepStd<T, D1, D2> > tmp = lhs * rhs;
0082 typedef SMatrix<T, D1, D1, MatRepSym<T, D1> > SMatrixSym;
0083 SMatrixSym mret;
0084 AssignSym::Evaluate(mret, tmp * Transpose(lhs));
0085 return mret;
0086 }
0087
0088 }
0089 }
0090
0091 template <typename T, IndexType N>
0092 inline void mult(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> >& a,
0093 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& rh,
0094 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& lh) {
0095
0096 for (IndexType i = 0; i != N; ++i) {
0097 IndexType off_i = a.fRep.offset(i, 0);
0098 for (IndexType k = 0; k != N; ++k) {
0099 IndexType off_k = a.fRep.offset(k, 0);
0100 if (k < i) {
0101 for (IndexType j = 0; j != (k + 1); ++j)
0102 a.Array()[off_i + j] += rh(i, k) * lh.Array()[off_k + j];
0103 for (IndexType j = k + 1; j != (i + 1); ++j)
0104 a.Array()[off_i + j] += rh(i, k) * lh(k, j);
0105 } else
0106 for (IndexType j = 0; j != (i + 1); ++j)
0107 a.Array()[off_i + j] += rh(i, k) * lh.Array()[off_k + j];
0108 }
0109 }
0110 }
0111
0112 template <typename T, IndexType N>
0113 inline void mult(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& a,
0114 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& rh,
0115 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& lh) {
0116
0117 for (IndexType i = 0; i != N; ++i) {
0118 for (IndexType j = 0; j <= i; ++j) {
0119
0120
0121 for (IndexType k = 0; k != N; ++k)
0122 a(i, j) += rh(i, k) * lh(j, k);
0123 }
0124 }
0125
0126 for (IndexType i = 0; i != N - 1; ++i)
0127 for (IndexType j = i + 1; j != N; ++j)
0128 a(i, j) = a(j, i);
0129 }
0130
0131
0132 template <typename T, IndexType N>
0133 inline void similarity(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& b,
0134 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& u,
0135 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& a) {
0136 for (IndexType i = 0; i != N; ++i)
0137 for (IndexType j = 0; j <= i; ++j)
0138 for (IndexType k = 0; k != N; ++k)
0139 for (IndexType l = 0; l != N; ++l)
0140 b(i, j) += u(i, k) * a(k, l) * u(j, l);
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155 for (IndexType i = 0; i != N - 1; ++i)
0156 for (IndexType j = i + 1; j != N; ++j)
0157 b(i, j) = b(j, i);
0158 }
0159
0160
0161 template <typename T, IndexType N>
0162 inline void similarityT(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& b,
0163 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& u,
0164 ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& a) {
0165 for (IndexType i = 0; i != N; ++i)
0166 for (IndexType j = 0; j <= i; ++j)
0167 for (IndexType k = 0; k != N; ++k)
0168 for (IndexType l = 0; l != N; ++l)
0169 b(i, j) += u(k, i) * a(k, l) * u(l, j);
0170
0171 for (IndexType i = 0; i != N - 1; ++i)
0172 for (IndexType j = i + 1; j != N; ++j)
0173 b(i, j) = b(j, i);
0174 }
0175
0176 template <typename M1, typename M2>
0177 double eps(M1 const& m1, M2 const& m2) {
0178 IndexType N = M1::kRows;
0179 double ret = 0.;
0180 for (IndexType i = 0; i != N; ++i)
0181 for (IndexType j = 0; j != N; ++j)
0182 ret = std::max(ret, std::abs(m1(i, j) - m2(i, j)));
0183 return ret;
0184 }
0185
0186 template <typename M>
0187 bool isSym(M const& m) {
0188 IndexType N = M::kRows;
0189 for (IndexType i = 0; i != N; ++i)
0190 for (IndexType j = 0; j <= i; ++j)
0191 if (m(i, j) != m(j, i))
0192 return false;
0193 return true;
0194 }
0195
0196 #include <iostream>
0197 #include "FWCore/Utilities/interface/HRRealTime.h"
0198
0199 #include <random>
0200
0201 namespace {
0202 std::mt19937 eng;
0203 std::uniform_real_distribution<double> rgen(-5., 5.);
0204 template <typename T, IndexType N>
0205 inline void fillRandom(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& a) {
0206 for (IndexType k = 0; k != N; ++k) {
0207 a(k, k) = std::abs(rgen(eng));
0208 for (IndexType j = 0; j < k; ++j)
0209 a(k, j) = rgen(eng);
0210 }
0211 for (IndexType i = 0; i != N - 1; ++i)
0212 for (IndexType j = i + 1; j != N; ++j)
0213 a(i, j) = -a(j, i);
0214 }
0215 }
0216
0217 template <typename T, IndexType N>
0218 void go(edm::HRTimeType& s1,
0219 edm::HRTimeType& s2,
0220 edm::HRTimeType& s3,
0221 edm::HRTimeType& s4,
0222 edm::HRTimeType& s5,
0223 edm::HRTimeType& s6,
0224 edm::HRTimeType& s7,
0225 edm::HRTimeType& s8,
0226 bool print) {
0227 typedef ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > Matrix;
0228 typedef ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > SymMatrix;
0229
0230 ROOT::Math::SMatrixIdentity id;
0231
0232 Matrix im1(id);
0233 Matrix im2(id);
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244 SymMatrix is1;
0245 ROOT::Math::AssignSym::Evaluate(is1, im1);
0246 SymMatrix is2;
0247 ROOT::Math::AssignSym::Evaluate(is2, im2);
0248
0249 Matrix j1;
0250 fillRandom(j1);
0251 Matrix j2;
0252 fillRandom(j2);
0253
0254 SymMatrix rh;
0255 s6 = edm::hrRealTime();
0256 rh = ROOT::Math::Similarity(j1, is1);
0257 s6 = edm::hrRealTime() - s6;
0258
0259 SymMatrix rh2;
0260 s7 = edm::hrRealTime();
0261 rh2 = ROOT::Math::Similarity(j1, im1);
0262 s7 = edm::hrRealTime() - s7;
0263
0264 Matrix mrh;
0265 s8 = edm::hrRealTime();
0266 similarity(mrh, j1, im1);
0267 s8 = edm::hrRealTime() - s8;
0268
0269 SymMatrix lh = ROOT::Math::SimilarityT(j2, is1);
0270
0271 Matrix mlh;
0272 similarityT(mlh, j2, im1);
0273
0274 SymMatrix a;
0275 s1 = edm::hrRealTime();
0276 mult(a, lh, rh);
0277 s1 = edm::hrRealTime() - s1;
0278
0279 SymMatrix b;
0280 s2 = edm::hrRealTime();
0281 ROOT::Math::AssignSym::Evaluate(b, lh * rh);
0282 s2 = edm::hrRealTime() - s2;
0283
0284 Matrix m0 = b;
0285
0286 Matrix m;
0287 s3 = edm::hrRealTime();
0288 m = mlh * mrh;
0289 s3 = edm::hrRealTime() - s3;
0290
0291 SymMatrix sm;
0292 s4 = edm::hrRealTime();
0293 ROOT::Math::AssignSym::Evaluate(sm, mlh * mrh);
0294 s4 = edm::hrRealTime() - s4;
0295
0296 Matrix m2;
0297 s5 = edm::hrRealTime();
0298 mult(m2, mlh, mrh);
0299 s5 = edm::hrRealTime() - s5;
0300
0301 SymMatrix smm;
0302 ROOT::Math::AssignSym::Evaluate(smm, m);
0303 SymMatrix smm2;
0304 ROOT::Math::AssignSym::Evaluate(smm2, m2);
0305
0306 if (print) {
0307 if (!isSym(im1))
0308 std::cout << " im is not sym" << std::endl;
0309 if (!isSym(mrh))
0310 std::cout << " rh is not sym" << std::endl;
0311 if (!isSym(mlh))
0312 std::cout << " lh is not sym" << std::endl;
0313 if (!isSym(m))
0314 std::cout << " m is not sym" << std::endl;
0315 if (!isSym(m2))
0316 std::cout << "m2 is not sym" << std::endl;
0317
0318 std::cout << "eps sim " << eps(rh, mrh) << std::endl;
0319 std::cout << "eps sim " << eps(rh, rh2) << std::endl;
0320 std::cout << "eps simT " << eps(lh, mlh) << std::endl;
0321 std::cout << "eps s m " << eps(m, b) << std::endl;
0322 std::cout << "eps s sm " << eps(m2, b) << std::endl;
0323
0324
0325
0326
0327 if (m != m0)
0328 std::cout << "problem with SMatrix Assign " << eps(m, m0) << std::endl;
0329 if (smm != b)
0330 std::cout << "problem with SMatrix * " << eps(smm, b) << std::endl;
0331 if (sm != b)
0332 std::cout << "problem with SMatrix evaluate " << eps(sm, b) << std::endl;
0333 if (a != b)
0334 std::cout << "problem with MulSymMatrix " << eps(a, b) << std::endl;
0335 if (a != smm)
0336 std::cout << "problem with MulSymMatrix twice " << eps(a, smm) << std::endl;
0337 if (m != m2)
0338 std::cout << "problem with MulMatrix " << eps(m, m2) << std::endl;
0339 if (smm != smm2)
0340 std::cout << "problem with MulMatrix twice " << eps(smm, smm2) << std::endl;
0341
0342 std::cout << "sym mult " << s1 << std::endl;
0343 std::cout << "sym * " << s2 << std::endl;
0344 std::cout << "mat * " << s3 << std::endl;
0345 std::cout << "mat as sym " << s4 << std::endl;
0346 std::cout << "mat mult " << s5 << std::endl;
0347 std::cout << "sym sim " << s6 << std::endl;
0348 std::cout << "std sim " << s7 << std::endl;
0349 std::cout << "loop sim " << s8 << std::endl;
0350 }
0351 }
0352
0353 int main() {
0354 edm::HRTimeType s1 = 0, s2 = 0, s3 = 0, s4 = 0, s5 = 0, s6 = 0, s7 = 0, s8 = 0;
0355 go<double, 3>(s1, s2, s3, s4, s5, s6, s7, s8, true);
0356
0357 go<double, 15>(s1, s2, s3, s4, s5, s6, s7, s8, true);
0358
0359 edm::HRTimeType t1 = 0;
0360 edm::HRTimeType t2 = 0;
0361 edm::HRTimeType t3 = 0;
0362 edm::HRTimeType t4 = 0;
0363 edm::HRTimeType t5 = 0;
0364 edm::HRTimeType t6 = 0;
0365 edm::HRTimeType t7 = 0;
0366 edm::HRTimeType t8 = 0;
0367
0368 for (int i = 0; i != 50000; ++i) {
0369 go<double, 15>(s1, s2, s3, s4, s5, s6, s7, s8, false);
0370 t1 += s1;
0371 t2 += s2;
0372 t3 += s3;
0373 t4 += s4;
0374 t5 += s5;
0375 t6 += s6;
0376 t7 += s7;
0377 t8 += s8;
0378 }
0379 std::cout << "sym mult " << t1 / 50000 << std::endl;
0380 std::cout << "sym * " << t2 / 50000 << std::endl;
0381 std::cout << "mat * " << t3 / 50000 << std::endl;
0382 std::cout << "mat as sym " << t4 / 50000 << std::endl;
0383 std::cout << "mat mult " << t5 / 50000 << std::endl;
0384 std::cout << "sym sim " << t6 / 50000 << std::endl;
0385 std::cout << "std sim " << t7 / 50000 << std::endl;
0386 std::cout << "loop sim " << t8 / 50000 << std::endl;
0387
0388 return 0;
0389 }
0390 #else
0391 int main() { return 0; }
0392 #endif