ROOT   Reference Guide
Searching...
No Matches
mnvert.cxx
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
10#include "Minuit2/MnMatrix.h"
11
12#include <cmath>
13
14namespace ROOT {
15
16namespace Minuit2 {
17
18/** Inverts a symmetric matrix. Matrix is first scaled to have all ones on
19 the diagonal (equivalent to change of units) but no pivoting is done
20 since matrix is positive-definite.
21 */
22
24{
25
26 unsigned int nrow = a.Nrow();
27 MnAlgebraicVector s(nrow);
28 MnAlgebraicVector q(nrow);
29 MnAlgebraicVector pp(nrow);
30
31 for (unsigned int i = 0; i < nrow; i++) {
32 double si = a(i, i);
33 if (si < 0.)
34 return 1;
35 s(i) = 1. / std::sqrt(si);
36 }
37
38 for (unsigned int i = 0; i < nrow; i++)
39 for (unsigned int j = i; j < nrow; j++)
40 a(i, j) *= (s(i) * s(j));
41
42 for (unsigned i = 0; i < nrow; i++) {
43 unsigned int k = i;
44 if (a(k, k) == 0.)
45 return 1;
46 q(k) = 1. / a(k, k);
47 pp(k) = 1.;
48 a(k, k) = 0.;
49 unsigned int kp1 = k + 1;
50 if (k != 0) {
51 for (unsigned int j = 0; j < k; j++) {
52 pp(j) = a(j, k);
53 q(j) = a(j, k) * q(k);
54 a(j, k) = 0.;
55 }
56 }
57 if (k != nrow - 1) {
58 for (unsigned int j = kp1; j < nrow; j++) {
59 pp(j) = a(k, j);
60 q(j) = -a(k, j) * q(k);
61 a(k, j) = 0.;
62 }
63 }
64 for (unsigned int j = 0; j < nrow; j++)
65 for (k = j; k < nrow; k++)
66 a(j, k) += (pp(j) * q(k));
67 }
68
69 for (unsigned int j = 0; j < nrow; j++)
70 for (unsigned int k = j; k < nrow; k++)
71 a(j, k) *= (s(j) * s(k));
72
73 return 0;
74}
75
76} // namespace Minuit2
77
78} // namespace ROOT
#define a(i)
Definition RSha256.hxx:99
float * q
Class describing a symmetric matrix of size n.
Definition LASymMatrix.h:45
int mnvert(LASymMatrix &t)
Inverts a symmetric matrix.
Definition mnvert.cxx:23
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...