Logo ROOT   6.18/05
Reference Guide
TDecompBase.cxx
Go to the documentation of this file.
1// @(#)root/matrix:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Dec 2003
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TDecompBase
13 \ingroup Matrix
14
15 Decomposition Base class
16
17 This class forms the base for all the decompositions methods in the
18 linear algebra package .
19 It or its derived classes have installed the methods to solve
20 equations,invert matrices and calculate determinants while monitoring
21 the accuracy.
22
23 Each derived class has always the following methods available:
24
25 #### Condition() :
26 In an iterative scheme the condition number for matrix inversion is
27 calculated . This number is of interest for estimating the accuracy
28 of x in the equation Ax=b
29 For example:
30 A is a (10x10) Hilbert matrix which looks deceivingly innocent
31 and simple, A(i,j) = 1/(i+j+1)
32 b(i) = Sum_j A(i,j), so a sum of a row in A
33
34 the solution is x(i) = 1. i=0,.,9
35
36 However,
37~~~
38 TMatrixD m....; TVectorD b.....
39 TDecompLU lu(m); lu.SetTol(1.0e-12); lu.Solve(b); b.Print()
40~~~
41 gives,
42
43~~~
44 {1.000,1.000,1.000,1.000,0.998,1.000,0.993,1.001,0.996,1.000}
45~~~
46
47 Looking at the condition number, this is in line with expected the
48 accuracy . The condition number is 3.957e+12 . As a simple rule of
49 thumb, a condition number of 1.0e+n means that you lose up to n
50 digits of accuracy in a solution . Since doubles are stored with 15
51 digits, we can expect the accuracy to be as small as 3 digits .
52
53 #### Det(Double_t &d1,Double_t &d2)
54 The determinant is d1*TMath::Power(2.,d2)
55 Expressing the determinant this way makes under/over-flow very
56 unlikely .
57
58 #### Decompose()
59 Here the actually decomposition is performed . One can change the
60 matrix A after the decomposition constructor has been called
61 without effecting the decomposition result
62
63 #### Solve(TVectorD &b)
64 Solve A x = b . x is supplied through the argument and replaced with
65 the solution .
66
67 #### TransSolve(TVectorD &b)
68 Solve A^T x = b . x is supplied through the argument and replaced
69 with the solution .
70
71 #### MultiSolve(TMatrixD &B)
72 Solve A X = B . where X and are now matrices . X is supplied through
73 the argument and replaced with the solution .
74
75 #### Invert(TMatrixD &inv)
76 This is of course just a call to MultiSolve with as input argument
77 the unit matrix . Note that for a matrix a(m,n) with m > n a
78 pseudo-inverse is calculated .
79
80 ### Tolerances and Scaling
81
82 The tolerance parameter (which is a member of this base class) plays
83 a crucial role in all operations of the decomposition classes . It
84 gives the user a powerful tool to monitor and steer the operations
85 Its default value is sqrt(epsilon) where 1+epsilon = 1
86
87 If you do not want to be bothered by the following considerations,
88 like in most other linear algebra packages, just set the tolerance
89 with SetTol to an arbitrary small number .
90
91 The tolerance number is used by each decomposition method to decide
92 whether the matrix is near singular, except of course SVD which can
93 handle singular matrices .
94 For each decomposition this will be checked in a different way; in LU
95 the matrix is considered singular when, at some point in the
96 decomposition, a diagonal element < fTol . Therefore, we had to set in
97 the example above of the (10x10) Hilbert, which is near singular, the
98 tolerance on 10e-12 . (The fact that we have to set the tolerance <
99 sqrt(epsilon) is a clear indication that we are losing precision .)
100
101 If the matrix is flagged as being singular, operations with the
102 decomposition will fail and will return matrices/vectors that are
103 invalid .
104
105 The observant reader will notice that by scaling the complete matrix
106 by some small number the decomposition will detect a singular matrix .
107 In this case the user will have to reduce the tolerance number by this
108 factor . (For CPU time saving we decided not to make this an automatic
109 procedure) .
110
111 Code for this could look as follows:
112~~~
113 const Double_t max_abs = Abs(a).Max();
114 const Double_t scale = TMath::Min(max_abs,1.);
115 a.SetTol(a.GetTol()*scale);
116~~~
117
118 For usage examples see $ROOTSYS/test/stressLinear.cxx
119*/
120
121#include "TDecompBase.h"
122#include "TMath.h"
123#include "TError.h"
124
126
127////////////////////////////////////////////////////////////////////////////////
128/// Default constructor
129
131{
133 fDet1 = 0;
134 fDet2 = 0;
135 fCondition = -1.0;
136 fRowLwb = 0;
137 fColLwb = 0;
138}
139
140////////////////////////////////////////////////////////////////////////////////
141/// Copy constructor
142
144{
145 *this = another;
146}
147
148////////////////////////////////////////////////////////////////////////////////
149
151{
152// Estimates lower bound for norm1 of inverse of A. Returns norm
153// estimate in est. iter sets the maximum number of iterations to be used.
154// The return value indicates the number of iterations remaining on exit from
155// loop, hence if this is non-zero the processed "converged".
156// This routine uses Hager's Convex Optimisation Algorithm.
157// See Applied Numerical Linear Algebra, p139 & SIAM J Sci Stat Comp 1984 pp 311-16
158
159 est = -1.0;
160
161 const TMatrixDBase &m = GetDecompMatrix();
162 if (!m.IsValid())
163 return iter;
164
165 const Int_t n = m.GetNrows();
166
167 TVectorD b(n); TVectorD y(n); TVectorD z(n);
168 b = Double_t(1.0/n);
169 Double_t inv_norm1 = 0.0;
170 Bool_t stop = kFALSE;
171 do {
172 y = b;
173 if (!Solve(y))
174 return iter;
175 const Double_t ynorm1 = y.Norm1();
176 if ( ynorm1 <= inv_norm1 ) {
177 stop = kTRUE;
178 } else {
179 inv_norm1 = ynorm1;
180 Int_t i;
181 for (i = 0; i < n; i++)
182 z(i) = ( y(i) >= 0.0 ? 1.0 : -1.0 );
183 if (!TransSolve(z))
184 return iter;
185 Int_t imax = 0;
186 Double_t maxz = TMath::Abs(z(0));
187 for (i = 1; i < n; i++) {
188 const Double_t absz = TMath::Abs(z(i));
189 if ( absz > maxz ) {
190 maxz = absz;
191 imax = i;
192 }
193 }
194 stop = (maxz <= b*z);
195 if (!stop) {
196 b = 0.0;
197 b(imax) = 1.0;
198 }
199 }
200 iter--;
201 } while (!stop && iter);
202 est = inv_norm1;
203
204 return iter;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208
210{
211// Returns product of matrix diagonal elements in d1 and d2. d1 is a mantissa and d2
212// an exponent for powers of 2. If matrix is in diagonal or triangular-matrix form this
213// will be the determinant.
214// Based on Bowler, Martin, Peters and Wilkinson in HACLA
215
216 const Double_t zero = 0.0;
217 const Double_t one = 1.0;
218 const Double_t four = 4.0;
219 const Double_t sixteen = 16.0;
220 const Double_t sixteenth = 0.0625;
221
222 const Int_t n = diag.GetNrows();
223
224 Double_t t1 = 1.0;
225 Double_t t2 = 0.0;
226 Int_t niter2 =0;
227 Int_t niter3 =0;
228 for (Int_t i = 0; (((i < n) && (t1 !=zero ))); i++) {
229 if (TMath::Abs(diag(i)) > tol) {
230 t1 *= (Double_t) diag(i);
231 while ( TMath::Abs(t1) < one) {
232 t1 *= sixteenth;
233 t2 += four;
234 niter2++;
235 if ( niter2>100) break;
236 }
237 while ( TMath::Abs(t1) < sixteenth) {
238 t1 *= sixteen;
239 t2 -= four;
240 niter3++;
241 if (niter3>100) break;
242 }
243 } else {
244 t1 = zero;
245 t2 = zero;
246 }
247 }
248 d1 = t1;
249 d2 = t2;
250
251 return;
252}
253
254////////////////////////////////////////////////////////////////////////////////
255/// Matrix condition number
256
258{
259 if ( !TestBit(kCondition) ) {
260 fCondition = -1;
261 if (TestBit(kSingular))
262 return fCondition;
263 if ( !TestBit(kDecomposed) ) {
264 if (!Decompose())
265 return fCondition;
266 }
267 Double_t invNorm;
268 if (Hager(invNorm))
269 fCondition *= invNorm;
270 else // no convergence in Hager
271 Error("Condition()","Hager procedure did NOT converge");
273 }
274 return fCondition;
275}
276
277////////////////////////////////////////////////////////////////////////////////
278/// Solve set of equations with RHS in columns of B
279
281{
282 const TMatrixDBase &m = GetDecompMatrix();
283 R__ASSERT(m.IsValid() && B.IsValid());
284
285 const Int_t colLwb = B.GetColLwb();
286 const Int_t colUpb = B.GetColUpb();
287 Bool_t status = kTRUE;
288 for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
289 TMatrixDColumn b(B,icol);
290 status &= Solve(b);
291 }
292
293 return status;
294}
295
296////////////////////////////////////////////////////////////////////////////////
297/// Matrix determinant det = d1*TMath::Power(2.,d2)
298
300{
301 if ( !TestBit(kDetermined) ) {
302 if ( !TestBit(kDecomposed) )
303 Decompose();
304 if (TestBit(kSingular) ) {
305 fDet1 = 0.0;
306 fDet2 = 0.0;
307 } else {
308 const TMatrixDBase &m = GetDecompMatrix();
309 R__ASSERT(m.IsValid());
310 TVectorD diagv(m.GetNrows());
311 for (Int_t i = 0; i < diagv.GetNrows(); i++)
312 diagv(i) = m(i,i);
313 DiagProd(diagv,fTol,fDet1,fDet2);
314 }
316 }
317 d1 = fDet1;
318 d2 = fDet2;
319}
320
321////////////////////////////////////////////////////////////////////////////////
322/// Print class members
323
324void TDecompBase::Print(Option_t * /*opt*/) const
325{
326 printf("fTol = %.4e\n",fTol);
327 printf("fDet1 = %.4e\n",fDet1);
328 printf("fDet2 = %.4e\n",fDet2);
329 printf("fCondition = %.4e\n",fCondition);
330 printf("fRowLwb = %d\n",fRowLwb);
331 printf("fColLwb = %d\n",fColLwb);
332}
333
334////////////////////////////////////////////////////////////////////////////////
335/// Assignment operator
336
338{
339 if (this != &source) {
340 TObject::operator=(source);
341 fTol = source.fTol;
342 fDet1 = source.fDet1;
343 fDet2 = source.fDet2;
344 fCondition = source.fCondition;
345 fRowLwb = source.fRowLwb;
346 fColLwb = source.fColLwb;
347 }
348 return *this;
349}
350
351////////////////////////////////////////////////////////////////////////////////
352/// Define a Householder-transformation through the parameters up and b .
353
355 Double_t tol)
356{
357 const Int_t n = vc.GetNrows();
358 const Double_t * const vp = vc.GetMatrixArray();
359
360 Double_t c = TMath::Abs(vp[lp]);
361 Int_t i;
362 for (i = l; i < n; i++)
363 c = TMath::Max(TMath::Abs(vp[i]),c);
364
365 up = 0.0;
366 beta = 0.0;
367 if (c <= tol) {
368// Warning("DefHouseHolder","max vector=%.4e < %.4e",c,tol);
369 return kFALSE;
370 }
371
372 Double_t sd = vp[lp]/c; sd *= sd;
373 for (i = l; i < n; i++) {
374 const Double_t tmp = vp[i]/c;
375 sd += tmp*tmp;
376 }
377
378 Double_t vpprim = c*TMath::Sqrt(sd);
379 if (vp[lp] > 0.) vpprim = -vpprim;
380 up = vp[lp]-vpprim;
381 beta = 1./(vpprim*up);
382
383 return kTRUE;
384}
385
386////////////////////////////////////////////////////////////////////////////////
387/// Apply Householder-transformation.
388
390 Int_t lp,Int_t l,TMatrixDRow &cr)
391{
392 const Int_t nv = vc.GetNrows();
393 const Int_t nc = (cr.GetMatrix())->GetNcols();
394
395 if (nv > nc) {
396 Error("ApplyHouseHolder(const TVectorD &,..,TMatrixDRow &)","matrix row too short");
397 return;
398 }
399
400 const Int_t inc_c = cr.GetInc();
401 const Double_t * const vp = vc.GetMatrixArray();
402 Double_t * cp = cr.GetPtr();
403
404 Double_t s = cp[lp*inc_c]*up;
405 Int_t i;
406 for (i = l; i < nv; i++)
407 s += cp[i*inc_c]*vp[i];
408
409 s = s*beta;
410 cp[lp*inc_c] += s*up;
411 for (i = l; i < nv; i++)
412 cp[i*inc_c] += s*vp[i];
413}
414
415////////////////////////////////////////////////////////////////////////////////
416/// Apply Householder-transformation.
417
420{
421 const Int_t nv = vc.GetNrows();
422 const Int_t nc = (cc.GetMatrix())->GetNrows();
423
424 if (nv > nc) {
425 Error("ApplyHouseHolder(const TVectorD &,..,TMatrixDRow &)","matrix column too short");
426 return;
427 }
428
429 const Int_t inc_c = cc.GetInc();
430 const Double_t * const vp = vc.GetMatrixArray();
431 Double_t * cp = cc.GetPtr();
432
433 Double_t s = cp[lp*inc_c]*up;
434 Int_t i;
435 for (i = l; i < nv; i++)
436 s += cp[i*inc_c]*vp[i];
437
438 s = s*beta;
439 cp[lp*inc_c] += s*up;
440 for (i = l; i < nv; i++)
441 cp[i*inc_c] += s*vp[i];
442}
443
444////////////////////////////////////////////////////////////////////////////////
445/// Apply Householder-transformation.
446
448 Int_t lp,Int_t l,TVectorD &cv)
449{
450 const Int_t nv = vc.GetNrows();
451 const Int_t nc = cv.GetNrows();
452
453 if (nv > nc) {
454 Error("ApplyHouseHolder(const TVectorD &,..,TVectorD &)","vector too short");
455 return;
456 }
457
458 const Double_t * const vp = vc.GetMatrixArray();
459 Double_t * cp = cv.GetMatrixArray();
460
461 Double_t s = cp[lp]*up;
462 Int_t i;
463 for (i = l; i < nv; i++)
464 s += cp[i]*vp[i];
465
466 s = s*beta;
467 cp[lp] += s*up;
468 for (i = l; i < nv; i++)
469 cp[i] += s*vp[i];
470}
471
472////////////////////////////////////////////////////////////////////////////////
473/// Defines a Givens-rotation by calculating 2 rotation parameters c and s.
474/// The rotation is defined with the vector components v1 and v2.
475
477{
478 const Double_t a1 = TMath::Abs(v1);
479 const Double_t a2 = TMath::Abs(v2);
480 if (a1 > a2) {
481 const Double_t w = v2/v1;
482 const Double_t q = TMath::Hypot(1.,w);
483 c = 1./q;
484 if (v1 < 0.) c = -c;
485 s = c*w;
486 } else {
487 if (v2 != 0) {
488 const Double_t w = v1/v2;
489 const Double_t q = TMath::Hypot(1.,w);
490 s = 1./q;
491 if (v2 < 0.) s = -s;
492 c = s*w;
493 } else {
494 c = 1.;
495 s = 0.;
496 }
497 }
498}
499
500////////////////////////////////////////////////////////////////////////////////
501/// Define and apply a Givens-rotation by calculating 2 rotation
502/// parameters c and s. The rotation is defined with and applied to the vector
503/// components v1 and v2.
504
506{
507 const Double_t a1 = TMath::Abs(v1);
508 const Double_t a2 = TMath::Abs(v2);
509 if (a1 > a2) {
510 const Double_t w = v2/v1;
511 const Double_t q = TMath::Hypot(1.,w);
512 c = 1./q;
513 if (v1 < 0.) c = -c;
514 s = c*w;
515 v1 = a1*q;
516 v2 = 0.;
517 } else {
518 if (v2 != 0) {
519 const Double_t w = v1/v2;
520 const Double_t q = TMath::Hypot(1.,w);
521 s = 1./q;
522 if (v2 < 0.) s = -s;
523 c = s*w;
524 v1 = a2*q;
525 v2 = 0.;
526 } else {
527 c = 1.;
528 s = 0.;
529 }
530 }
531}
532
533////////////////////////////////////////////////////////////////////////////////
534/// Apply a Givens transformation as defined by c and s to the vector components
535/// v1 and v2 .
536
538{
539 const Double_t w = z1*c+z2*s;
540 z2 = -z1*s+z2*c;
541 z1 = w;
542}
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
Bool_t DefHouseHolder(const TVectorD &vc, Int_t lp, Int_t l, Double_t &up, Double_t &beta, Double_t tol)
Define a Householder-transformation through the parameters up and b .
void ApplyHouseHolder(const TVectorD &vc, Double_t up, Double_t beta, Int_t lp, Int_t l, TMatrixDRow &cr)
Apply Householder-transformation.
void ApplyGivens(Double_t &z1, Double_t &z2, Double_t c, Double_t s)
Apply a Givens transformation as defined by c and s to the vector components v1 and v2 .
void DefGivens(Double_t v1, Double_t v2, Double_t &c, Double_t &s)
Defines a Givens-rotation by calculating 2 rotation parameters c and s.
void DefAplGivens(Double_t &v1, Double_t &v2, Double_t &c, Double_t &s)
Define and apply a Givens-rotation by calculating 2 rotation parameters c and s.
#define R__ASSERT(e)
Definition: TError.h:96
void Error(const char *location, const char *msgfmt,...)
float * q
Definition: THbookFile.cxx:87
Decomposition Base class.
Definition: TDecompBase.h:34
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
Double_t fDet1
Definition: TDecompBase.h:37
virtual Bool_t Decompose()=0
Double_t fDet2
Definition: TDecompBase.h:38
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
void Print(Option_t *opt="") const
Print class members.
virtual const TMatrixDBase & GetDecompMatrix() const =0
virtual Bool_t Solve(TVectorD &b)=0
Int_t fRowLwb
Definition: TDecompBase.h:40
Double_t fTol
Definition: TDecompBase.h:36
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Int_t Hager(Double_t &est, Int_t iter=5)
virtual Double_t Condition()
Matrix condition number.
TDecompBase()
Default constructor.
Double_t fCondition
Definition: TDecompBase.h:39
virtual Bool_t TransSolve(TVectorD &b)=0
Int_t fColLwb
Definition: TDecompBase.h:41
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetInc() const
Element * GetPtr() const
Int_t GetInc() const
const TMatrixTBase< Element > * GetMatrix() const
Element * GetPtr() const
Mother of all ROOT objects.
Definition: TObject.h:37
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:268
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
Int_t GetNrows() const
Definition: TVectorT.h:75
Element * GetMatrixArray()
Definition: TVectorT.h:78
double beta(double x, double y)
Calculates the beta function.
Double_t y[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
static double B[]
static constexpr double s
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:57
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * t1
Definition: textangle.C:20
REAL epsilon
Definition: triangle.c:617