Logo ROOT  
Reference Guide
Dinv.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_Dinv
5 #define ROOT_Math_Dinv
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 03. Apr 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: square Matrix inversion
23 // Code was taken from CERNLIB::kernlib dfinv function, translated
24 // from FORTRAN to C++ and optimized.
25 // n: Order of the square matrix
26 // idim: First dimension of array A
27 //
28 // changes:
29 // 03 Apr 2001 (TG) creation
30 //
31 // ********************************************************************
32 #ifdef OLD_IMPL
33 #include "Math/Dfactir.h"
34 #include "Math/Dfinv.h"
35 #include "Math/Dsinv.h"
36 #endif
37 
38 #include "Math/CholeskyDecomp.h"
39 
41 
42 // #ifndef ROOT_Math_QRDecomposition
43 // #include "Math/QRDecomposition.h"
44 // #endif
45 
46 #include "TError.h"
47 
48 namespace ROOT {
49 
50  namespace Math {
51 
52 
53 
54 /**
55  Matrix Inverter class
56  Class to specialize calls to Dinv. Dinv computes the inverse of a square
57  matrix if dimension idim and order n. The content of the matrix will be
58  replaced by its inverse. In case the inversion fails, the matrix content is
59  destroyed. Invert specializes Dinv by the matrix order. E.g. if the order
60  of the matrix is two, the routine Inverter<2> is called which implements
61  Cramers rule.
62 
63  @author T. Glebe
64 */
65 //==============================================================================
66 // Inverter class
67 //==============================================================================
68 template <unsigned int idim, unsigned int n = idim>
69 class Inverter {
70 public:
71  /// matrix inversion for a generic square matrix using LU factorization
72  /// (code originally from CERNLIB and then ported in C++ for CLHEP)
73  /// implementation is in file Math/MatrixInversion.icc
74  template <class MatrixRep>
75  static bool Dinv(MatrixRep& rhs) {
76 
77 
78  /* Initialized data */
79  unsigned int work[n+1] = {0};
80 
81  typename MatrixRep::value_type det(0.0);
82 
83  if (DfactMatrix(rhs,det,work) != 0) {
84  Error("Inverter::Dinv","Dfact_matrix failed!!");
85  return false;
86  }
87 
88  int ifail = DfinvMatrix(rhs,work);
89  if (ifail == 0) return true;
90  return false;
91  } // Dinv
92 
93 
94  /// symmetric matrix inversion using
95  /// Bunch-kaufman pivoting method
96  /// implementation in Math/MatrixInversion.icc
97  template <class T>
98  static bool Dinv(MatRepSym<T,idim> & rhs) {
99  int ifail = 0;
100  InvertBunchKaufman(rhs,ifail);
101  if (ifail == 0) return true;
102  return false;
103  }
104 
105 
106  /**
107  LU Factorization method for inversion of general square matrices
108  (see implementation in Math/MatrixInversion.icc)
109  */
110  template <class T>
111  static int DfactMatrix(MatRepStd<T,idim,n> & rhs, T & det, unsigned int * work);
112  /**
113  LU inversion of general square matrices. To be called after DFactMatrix
114  (see implementation in Math/MatrixInversion.icc)
115  */
116  template <class T>
117  static int DfinvMatrix(MatRepStd<T,idim,n> & rhs, unsigned int * work);
118 
119  /**
120  Bunch-Kaufman method for inversion of symmetric matrices
121  */
122  template <class T>
123  static void InvertBunchKaufman(MatRepSym<T,idim> & rhs, int &ifail);
124 
125 
126 
127 }; // class Inverter
128 
129 // fast inverter class using Cramer inversion
130 // by default use other default inversion
131 /**
132  Fast Matrix Inverter class
133  Class to specialize calls to Dinv. Dinv computes the inverse of a square
134  matrix if dimension idim and order n. The content of the matrix will be
135  replaced by its inverse. In case the inversion fails, the matrix content is
136  destroyed. Invert specializes Dinv by the matrix order. E.g. if the order
137  of the matrix is less than 5 , the class implements
138  Cramers rule.
139  Be careful that for matrix with high condition the accuracy of the Cramer rule is much poorer
140 
141  @author L. Moneta
142 */
143 template <unsigned int idim, unsigned int n = idim>
145 public:
146  ///
147  template <class MatrixRep>
148  static bool Dinv(MatrixRep& rhs) {
149  return Inverter<idim,n>::Dinv(rhs);
150  }
151  template <class T>
152  static bool Dinv(MatRepSym<T,idim> & rhs) {
153  return Inverter<idim,n>::Dinv(rhs);
154  }
155 };
156 
157 
158 /** Inverter<0>.
159  In case of zero order, do nothing.
160 
161  @author T. Glebe
162 */
163 //==============================================================================
164 // Inverter<0>
165 //==============================================================================
166 template <>
167 class Inverter<0> {
168 public:
169  ///
170  template <class MatrixRep>
171  inline static bool Dinv(MatrixRep&) { return true; }
172 };
173 
174 
175 /**
176  1x1 matrix inversion \f$a_{11} \to 1/a_{11}\f$
177 
178  @author T. Glebe
179 */
180 //==============================================================================
181 // Inverter<1>
182 //==============================================================================
183 template <>
184 class Inverter<1> {
185 public:
186  ///
187  template <class MatrixRep>
188  static bool Dinv(MatrixRep& rhs) {
189 
190  if (rhs[0] == 0.) {
191  return false;
192  }
193  rhs[0] = 1. / rhs[0];
194  return true;
195  }
196 };
197 
198 
199 /**
200  2x2 matrix inversion using Cramers rule.
201 
202  @author T. Glebe
203 */
204 //==============================================================================
205 // Inverter<2>: Cramers rule
206 //==============================================================================
207 
208 template <>
209 class Inverter<2> {
210 public:
211  ///
212  template <class MatrixRep>
213  static bool Dinv(MatrixRep& rhs) {
214 
215  typedef typename MatrixRep::value_type T;
216  T det = rhs[0] * rhs[3] - rhs[2] * rhs[1];
217 
218  if (det == T(0.) ) { return false; }
219 
220  T s = T(1.0) / det;
221 
222  T c11 = s * rhs[3];
223 
224 
225  rhs[2] = -s * rhs[2];
226  rhs[1] = -s * rhs[1];
227  rhs[3] = s * rhs[0];
228  rhs[0] = c11;
229 
230 
231  return true;
232  }
233 
234  // specialization for the symmetric matrices
235  template <class T>
236  static bool Dinv(MatRepSym<T,2> & rep) {
237 
238  T * rhs = rep.Array();
239 
240  T det = rhs[0] * rhs[2] - rhs[1] * rhs[1];
241 
242 
243  if (det == T(0.)) { return false; }
244 
245  T s = T(1.0) / det;
246  T c11 = s * rhs[2];
247 
248  rhs[1] = -s * rhs[1];
249  rhs[2] = s * rhs[0];
250  rhs[0] = c11;
251  return true;
252  }
253 
254 };
255 
256 
257 /**
258  3x3 direct matrix inversion using Cramer Rule
259  use only for FastInverter
260 */
261 //==============================================================================
262 // FastInverter<3>
263 //==============================================================================
264 
265 template <>
266 class FastInverter<3> {
267 public:
268  ///
269  // use Cramer Rule
270  template <class MatrixRep>
271  static bool Dinv(MatrixRep& rhs);
272 
273  template <class T>
274  static bool Dinv(MatRepSym<T,3> & rhs);
275 
276 };
277 
278 /**
279  4x4 matrix inversion using Cramers rule.
280 */
281 template <>
282 class FastInverter<4> {
283 public:
284  ///
285  template <class MatrixRep>
286  static bool Dinv(MatrixRep& rhs);
287 
288  template <class T>
289  static bool Dinv(MatRepSym<T,4> & rhs);
290 
291 };
292 
293 /**
294  5x5 Matrix inversion using Cramers rule.
295 */
296 template <>
297 class FastInverter<5> {
298 public:
299  ///
300  template <class MatrixRep>
301  static bool Dinv(MatrixRep& rhs);
302 
303  template <class T>
304  static bool Dinv(MatRepSym<T,5> & rhs);
305 
306 };
307 
308 // inverter for Cholesky
309 // works only for symmetric matrices and will produce a
310 // compilation error otherwise
311 
312 template <unsigned int idim>
314 public:
315  ///
316  template <class MatrixRep>
317  static bool Dinv(MatrixRep&) {
318  STATIC_CHECK( false, Error_cholesky_SMatrix_type_is_not_symmetric );
319  return false;
320  }
321  template <class T>
322  inline static bool Dinv(MatRepSym<T,idim> & rhs) {
323  CholeskyDecomp<T, idim> decomp(rhs);
324  return decomp.Invert(rhs);
325  }
326 };
327 
328 
329  } // namespace Math
330 
331 } // namespace ROOT
332 
333 #include "CramerInversion.icc"
334 #include "CramerInversionSym.icc"
335 #include "MatrixInversion.icc"
336 
337 #endif /* ROOT_Math_Dinv */
ROOT::Math::Inverter::DfactMatrix
static int DfactMatrix(MatRepStd< T, idim, n > &rhs, T &det, unsigned int *work)
LU Factorization method for inversion of general square matrices (see implementation in Math/MatrixIn...
Definition: MatrixInversion.icc:447
n
const Int_t n
Definition: legend1.C:16
MatrixRepresentationsStatic.h
ROOT::Math::CholeskyDecomp
class to compute the Cholesky decomposition of a matrix
Definition: CholeskyDecomp.h:76
ROOT::Math::CholInverter
Definition: Dinv.h:313
CramerInversionSym.icc
ROOT::Math::MatRepSym
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
ROOT::Math::FastInverter::Dinv
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:148
ROOT::Math::Inverter< 2 >::Dinv
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:213
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:168
ROOT::Math::MatRepSym::Array
T * Array()
Definition: MatrixRepresentationsStatic.h:239
ROOT::Math::Inverter::Dinv
static bool Dinv(MatRepSym< T, idim > &rhs)
symmetric matrix inversion using Bunch-kaufman pivoting method implementation in Math/MatrixInversion...
Definition: Dinv.h:98
CholeskyDecomp.h
CramerInversion.icc
ROOT::Math::FastInverter::Dinv
static bool Dinv(MatRepSym< T, idim > &rhs)
Definition: Dinv.h:152
ROOT::Math::MatRepStd
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
ROOT::Math::Inverter::DfinvMatrix
static int DfinvMatrix(MatRepStd< T, idim, n > &rhs, unsigned int *work)
LU inversion of general square matrices.
Definition: MatrixInversion.icc:577
ROOT::Math::Inverter::Dinv
static bool Dinv(MatrixRep &rhs)
matrix inversion for a generic square matrix using LU factorization (code originally from CERNLIB and...
Definition: Dinv.h:75
ROOT::Math::FastInverter
Fast Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:144
ROOT::Math::Inverter
Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:69
MatrixInversion.icc
ROOT::Math::CholInverter::Dinv
static bool Dinv(MatRepSym< T, idim > &rhs)
Definition: Dinv.h:322
ROOT::Math::Inverter< 0 >::Dinv
static bool Dinv(MatrixRep &)
Definition: Dinv.h:171
ROOT::Math::Inverter< 1 >::Dinv
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:188
ROOT::Math::CholInverter::Dinv
static bool Dinv(MatrixRep &)
Definition: Dinv.h:317
ROOT::Math::Inverter< 2 >::Dinv
static bool Dinv(MatRepSym< T, 2 > &rep)
Definition: Dinv.h:236
ROOT::Math::Chebyshev::T
double T(double x)
Definition: ChebyshevPol.h:52
Dfactir.h
STATIC_CHECK
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
ROOT::Math::Inverter::InvertBunchKaufman
static void InvertBunchKaufman(MatRepSym< T, idim > &rhs, int &ifail)
Bunch-Kaufman method for inversion of symmetric matrices.
Definition: MatrixInversion.icc:40
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Dfinv.h
Math
Dsinv.h
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:188
TError.h
ROOT::Math::CholeskyDecomp::Invert
bool Invert(M &m) const
place the inverse into m
Definition: CholeskyDecomp.h:149