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