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
48namespace 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//==============================================================================
68template <unsigned int idim, unsigned int n = idim>
69class Inverter {
70public:
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*/
143template <unsigned int idim, unsigned int n = idim>
145public:
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//==============================================================================
166template <>
167class Inverter<0> {
168public:
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//==============================================================================
183template <>
184class Inverter<1> {
185public:
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
208template <>
209class Inverter<2> {
210public:
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
265template <>
266class FastInverter<3> {
267public:
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*/
281template <>
282class FastInverter<4> {
283public:
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*/
296template <>
297class FastInverter<5> {
298public:
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
312template <unsigned int idim>
314public:
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 */
header file containing the templated implementation of matrix inversion routines for use with ROOT's ...
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
static bool Dinv(MatRepSym< T, idim > &rhs)
Definition: Dinv.h:322
static bool Dinv(MatrixRep &)
Definition: Dinv.h:317
class to compute the Cholesky decomposition of a matrix
bool Invert(M &m) const
place the inverse into m
Fast Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:144
static bool Dinv(MatRepSym< T, idim > &rhs)
Definition: Dinv.h:152
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:148
static bool Dinv(MatrixRep &)
Definition: Dinv.h:171
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:188
static bool Dinv(MatRepSym< T, 2 > &rep)
Definition: Dinv.h:236
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:213
Matrix Inverter class Class to specialize calls to Dinv.
Definition: Dinv.h:69
static bool Dinv(MatrixRep &rhs)
matrix inversion for a generic square matrix using LU factorization (code originally from CERNLIB and...
Definition: Dinv.h:75
static void InvertBunchKaufman(MatRepSym< T, idim > &rhs, int &ifail)
Bunch-Kaufman method for inversion of symmetric matrices.
static bool Dinv(MatRepSym< T, idim > &rhs)
symmetric matrix inversion using Bunch-kaufman pivoting method implementation in Math/MatrixInversion...
Definition: Dinv.h:98
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...
static int DfinvMatrix(MatRepStd< T, idim, n > &rhs, unsigned int *work)
LU inversion of general square matrices.
Expression wrapper class for Matrix objects.
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
double T(double x)
Definition: ChebyshevPol.h:34
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static constexpr double s