Logo ROOT   6.14/05
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 
19 /** Inverts a symmetric matrix. Matrix is first scaled to have all ones on
20  the diagonal (equivalent to change of units) but no pivoting is done
21  since matrix is positive-definite.
22  */
23 
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.) return 1;
34  s(i) = 1./std::sqrt(si);
35  }
36 
37  for(unsigned int i = 0; i < nrow; i++)
38  for(unsigned int j = i; j < nrow; j++)
39  a(i,j) *= (s(i)*s(j));
40 
41  for(unsigned i = 0; i < nrow; i++) {
42  unsigned int k = i;
43  if(a(k,k) == 0.) return 1;
44  q(k) = 1./a(k,k);
45  pp(k) = 1.;
46  a(k,k) = 0.;
47  unsigned int kp1 = k + 1;
48  if(k != 0) {
49  for(unsigned int j = 0; j < k; j++) {
50  pp(j) = a(j,k);
51  q(j) = a(j,k)*q(k);
52  a(j,k) = 0.;
53  }
54  }
55  if (k != nrow-1) {
56  for(unsigned int j = kp1; j < nrow; j++) {
57  pp(j) = a(k,j);
58  q(j) = -a(k,j)*q(k);
59  a(k,j) = 0.;
60  }
61  }
62  for(unsigned int j = 0; j < nrow; j++)
63  for(k = j; k < nrow; k++)
64  a(j,k) += (pp(j)*q(k));
65  }
66 
67  for(unsigned int j = 0; j < nrow; j++)
68  for(unsigned int k = j; k < nrow; k++)
69  a(j,k) *= (s(j)*s(k));
70 
71  return 0;
72 }
73 
74  } // namespace Minuit2
75 
76 } // namespace ROOT
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
int mnvert(LASymMatrix &t)
Inverts a symmetric matrix.
Definition: mnvert.cxx:24
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
double sqrt(double)
auto * a
Definition: textangle.C:12
unsigned int Nrow() const
Definition: LASymMatrix.h:239
static constexpr double s
float * q
Definition: THbookFile.cxx:87