Logo ROOT   6.12/07
Reference Guide
Dsinv.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dsinv
5 #define ROOT_Math_Dsinv
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 22. Mar 2001
13 //
14 // author: Thorsten Glebe
15 // HERA-B Collaboration
16 // Max-Planck-Institut fuer Kernphysik
17 // Saupfercheckweg 1
18 // 69117 Heidelberg
19 // Germany
20 // E-mail: T.Glebe@mpi-hd.mpg.de
21 //
22 // Description: Inversion of a symmetric, positive definite matrix.
23 // Code was taken from CERNLIB::kernlib dsinv function, translated
24 // from FORTRAN to C++ and optimized.
25 //
26 // changes:
27 // 22 Mar 2001 (TG) creation
28 //
29 // ********************************************************************
30 
31 #include "Math/SMatrixDfwd.h"
32 
33 namespace ROOT {
34 
35  namespace Math {
36 
37 /** Dsinv.
38  Compute inverse of a symmetric, positive definite matrix of dimension
39  $idim$ and order $n$.
40 
41  @author T. Glebe
42 */
43 template <class T, int n, int idim>
44 class SInverter
45 {
46 
47 public:
48  template <class MatrixRep>
49  inline static bool Dsinv(MatrixRep& rhs) {
50 
51  /* Local variables */
52  int i, j, k, l;
53  T s31, s32;
54  int jm1, jp1;
55 
56  /* Parameter adjustments */
57  const int arrayOffset = -1*(idim + 1);
58 
59 
60  /* Function Body */
61  if (idim < n || n <= 1) {
62  return false;
63  }
64 
65  /* sfact.inc */
66  for (j = 1; j <= n; ++j) {
67  const int ja = j * idim;
68  const int jj = j + ja;
69  const int ja1 = ja + idim;
70 
71 
72  if (rhs[jj + arrayOffset] <= 0.) { return false; }
73  rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
74  if (j == n) { break; }
75 
76  for (l = j + 1; l <= n; ++l) {
77  rhs[j + (l * idim) + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ja + arrayOffset];
78  const int lj = l + ja1;
79  for (i = 1; i <= j; ++i) {
80  rhs[lj + arrayOffset] -= rhs[l + (i * idim) + arrayOffset] * rhs[i + ja1 + arrayOffset];
81  }
82  }
83  }
84 
85  /* sfinv.inc */
86  // idim << 1 is equal to idim * 2
87  // compiler will compute the arguments!
88  rhs[((idim << 1) + 1) + arrayOffset] = -rhs[((idim << 1) + 1) + arrayOffset];
89  rhs[idim + 2 + arrayOffset] = rhs[((idim << 1)) + 1 + arrayOffset] * rhs[((idim << 1)) + 2 + arrayOffset];
90 
91  if(n > 2) {
92 
93  for (j = 3; j <= n; ++j) {
94  const int jm2 = j - 2;
95  const int ja = j * idim;
96  const int jj = j + ja;
97  const int j1 = j - 1 + ja;
98 
99  for (k = 1; k <= jm2; ++k) {
100  s31 = rhs[k + ja + arrayOffset];
101 
102  for (i = k; i <= jm2; ++i) {
103  s31 += rhs[k + ((i + 1) * idim) + arrayOffset] * rhs[i + 1 + ja + arrayOffset];
104  } // for i
105  rhs[k + ja + arrayOffset] = -s31;
106  rhs[j + (k * idim) + arrayOffset] = -s31 * rhs[jj + arrayOffset];
107  } // for k
108  rhs[j1 + arrayOffset] *= -1;
109  // rhs[j1] = -rhs[j1];
110  rhs[jj - idim + arrayOffset] = rhs[j1 + arrayOffset] * rhs[jj + arrayOffset];
111  } // for j
112  } // if (n>2)
113 
114  j = 1;
115  do {
116  const int jad = j * idim;
117  const int jj = j + jad;
118 
119  jp1 = j + 1;
120  for (i = jp1; i <= n; ++i) {
121  rhs[jj + arrayOffset] += rhs[j + (i * idim) + arrayOffset] * rhs[i + jad + arrayOffset];
122  } // for i
123 
124  jm1 = j;
125  j = jp1;
126  const int ja = j * idim;
127 
128  for (k = 1; k <= jm1; ++k) {
129  s32 = 0.;
130  for (i = j; i <= n; ++i) {
131  s32 += rhs[k + (i * idim) + arrayOffset] * rhs[i + ja + arrayOffset];
132  } // for i
133  //rhs[k + ja + arrayOffset] = rhs[j + (k * idim) + arrayOffset] = s32;
134  rhs[k + ja + arrayOffset] = s32;
135  } // for k
136  } while(j < n);
137 
138  return true;
139  }
140 
141 
142  // for symmetric matrices
143 
144  static bool Dsinv(MatRepSym<T,n> & rhs) {
145  // not very efficient but need to re-do Dsinv for new storage of
146  // symmetric matrices
147  MatRepStd<T,n,n> tmp;
148  for (int i = 0; i< n*n; ++i)
149  tmp[i] = rhs[i];
150  // call dsinv
151  if (! SInverter<T,n,n>::Dsinv(tmp) ) return false;
152  //if (! Inverter<n>::Dinv(tmp) ) return false;
153  // recopy the data
154  for (int i = 0; i< n*n; ++i)
155  rhs[i] = tmp[i];
156 
157  return true;
158 
159  }
160 
161 }; // end of Dsinv
162 
163 
164 
165  } // namespace Math
166 
167 } // namespace ROOT
168 
169 
170 #endif /* ROOT_Math_Dsinv */
static bool Dsinv(MatRepSym< T, n > &rhs)
Definition: Dsinv.h:144
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double T(double x)
Definition: ChebyshevPol.h:34
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
static bool Dsinv(MatrixRep &rhs)
Definition: Dsinv.h:49
Namespace for new Math classes and functions.
auto * l
Definition: textangle.C:4
const Int_t n
Definition: legend1.C:16