Logo ROOT  

Reference Guide
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 
14 namespace ROOT {
15 
16 namespace 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 
23 int mnvert(MnAlgebraicSymMatrix &a)
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
ROOT::Minuit2::LAVector
Definition: LAVector.h:32
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
MnMatrix.h
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
q
float * q
Definition: THbookFile.cxx:89
a
auto * a
Definition: textangle.C:12
sqrt
double sqrt(double)
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
ROOT::Minuit2::mnvert
int mnvert(LASymMatrix &t)
Inverts a symmetric matrix.
Definition: mnvert.cxx:23