Logo ROOT   6.08/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 #ifndef ROOT_Math_SMatrixDfwd
32 #include "Math/SMatrixDfwd.h"
33 #endif
34 
35 namespace ROOT {
36 
37  namespace Math {
38 
39 /** Dsinv.
40  Compute inverse of a symmetric, positive definite matrix of dimension
41  $idim$ and order $n$.
42 
43  @author T. Glebe
44 */
45 template <class T, int n, int idim>
46 class SInverter
47 {
48 
49 public:
50  template <class MatrixRep>
51  inline static bool Dsinv(MatrixRep& rhs) {
52 
53  /* Local variables */
54  int i, j, k, l;
55  T s31, s32;
56  int jm1, jp1;
57 
58  /* Parameter adjustments */
59  const int arrayOffset = -1*(idim + 1);
60 
61 
62  /* Function Body */
63  if (idim < n || n <= 1) {
64  return false;
65  }
66 
67  /* sfact.inc */
68  for (j = 1; j <= n; ++j) {
69  const int ja = j * idim;
70  const int jj = j + ja;
71  const int ja1 = ja + idim;
72 
73 
74  if (rhs[jj + arrayOffset] <= 0.) { return false; }
75  rhs[jj + arrayOffset] = 1. / rhs[jj + arrayOffset];
76  if (j == n) { break; }
77 
78  for (l = j + 1; l <= n; ++l) {
79  rhs[j + (l * idim) + arrayOffset] = rhs[jj + arrayOffset] * rhs[l + ja + arrayOffset];
80  const int lj = l + ja1;
81  for (i = 1; i <= j; ++i) {
82  rhs[lj + arrayOffset] -= rhs[l + (i * idim) + arrayOffset] * rhs[i + ja1 + arrayOffset];
83  }
84  }
85  }
86 
87  /* sfinv.inc */
88  // idim << 1 is equal to idim * 2
89  // compiler will compute the arguments!
90  rhs[((idim << 1) + 1) + arrayOffset] = -rhs[((idim << 1) + 1) + arrayOffset];
91  rhs[idim + 2 + arrayOffset] = rhs[((idim << 1)) + 1 + arrayOffset] * rhs[((idim << 1)) + 2 + arrayOffset];
92 
93  if(n > 2) {
94 
95  for (j = 3; j <= n; ++j) {
96  const int jm2 = j - 2;
97  const int ja = j * idim;
98  const int jj = j + ja;
99  const int j1 = j - 1 + ja;
100 
101  for (k = 1; k <= jm2; ++k) {
102  s31 = rhs[k + ja + arrayOffset];
103 
104  for (i = k; i <= jm2; ++i) {
105  s31 += rhs[k + ((i + 1) * idim) + arrayOffset] * rhs[i + 1 + ja + arrayOffset];
106  } // for i
107  rhs[k + ja + arrayOffset] = -s31;
108  rhs[j + (k * idim) + arrayOffset] = -s31 * rhs[jj + arrayOffset];
109  } // for k
110  rhs[j1 + arrayOffset] *= -1;
111  // rhs[j1] = -rhs[j1];
112  rhs[jj - idim + arrayOffset] = rhs[j1 + arrayOffset] * rhs[jj + arrayOffset];
113  } // for j
114  } // if (n>2)
115 
116  j = 1;
117  do {
118  const int jad = j * idim;
119  const int jj = j + jad;
120 
121  jp1 = j + 1;
122  for (i = jp1; i <= n; ++i) {
123  rhs[jj + arrayOffset] += rhs[j + (i * idim) + arrayOffset] * rhs[i + jad + arrayOffset];
124  } // for i
125 
126  jm1 = j;
127  j = jp1;
128  const int ja = j * idim;
129 
130  for (k = 1; k <= jm1; ++k) {
131  s32 = 0.;
132  for (i = j; i <= n; ++i) {
133  s32 += rhs[k + (i * idim) + arrayOffset] * rhs[i + ja + arrayOffset];
134  } // for i
135  //rhs[k + ja + arrayOffset] = rhs[j + (k * idim) + arrayOffset] = s32;
136  rhs[k + ja + arrayOffset] = s32;
137  } // for k
138  } while(j < n);
139 
140  return true;
141  }
142 
143 
144  // for symmetric matrices
145 
146  static bool Dsinv(MatRepSym<T,n> & rhs) {
147  // not very efficient but need to re-do Dsinv for new storage of
148  // symmetric matrices
149  MatRepStd<T,n,n> tmp;
150  for (int i = 0; i< n*n; ++i)
151  tmp[i] = rhs[i];
152  // call dsinv
153  if (! SInverter<T,n,n>::Dsinv(tmp) ) return false;
154  //if (! Inverter<n>::Dinv(tmp) ) return false;
155  // recopy the data
156  for (int i = 0; i< n*n; ++i)
157  rhs[i] = tmp[i];
158 
159  return true;
160 
161  }
162 
163 }; // end of Dsinv
164 
165 
166 
167  } // namespace Math
168 
169 } // namespace ROOT
170 
171 
172 #endif /* ROOT_Math_Dsinv */
static bool Dsinv(MatRepSym< T, n > &rhs)
Definition: Dsinv.h:146
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
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:51
TLine * l
Definition: textangle.C:4
Namespace for new Math classes and functions.
const Int_t n
Definition: legend1.C:16