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
14namespace 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
float * q
Definition: THbookFile.cxx:87
double sqrt(double)
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
int mnvert(LASymMatrix &t)
Inverts a symmetric matrix.
Definition: mnvert.cxx:24
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
static constexpr double s
auto * a
Definition: textangle.C:12