Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
125
126////////////////////////////////////////////////////////////////////////////////
127/// Default constructor
128
130{
131 fTol = std::numeric_limits<double>::epsilon();
132 fDet1 = 0;
133 fDet2 = 0;
134 fCondition = -1.0;
135 fRowLwb = 0;
136 fColLwb = 0;
137}
138
139////////////////////////////////////////////////////////////////////////////////
140/// Copy constructor
141
146
147////////////////////////////////////////////////////////////////////////////////
148
150{
151// Estimates lower bound for norm1 of inverse of A. Returns norm
152// estimate in est. iter sets the maximum number of iterations to be used.
153// The return value indicates the number of iterations remaining on exit from
154// loop, hence if this is non-zero the processed "converged".
155// This routine uses Hager's Convex Optimisation Algorithm.
156// See Applied Numerical Linear Algebra, p139 & SIAM J Sci Stat Comp 1984 pp 311-16
157
158 est = -1.0;
159
160 const TMatrixDBase &m = GetDecompMatrix();
161 if (!m.IsValid())
162 return iter;
163
164 const Int_t n = m.GetNrows();
165
166 TVectorD b(n); TVectorD y(n); TVectorD z(n);
167 b = Double_t(1.0/n);
168 Double_t inv_norm1 = 0.0;
169 Bool_t stop = kFALSE;
170 do {
171 y = b;
172 if (!Solve(y))
173 return iter;
174 const Double_t ynorm1 = y.Norm1();
175 if ( ynorm1 <= inv_norm1 ) {
176 stop = kTRUE;
177 } else {
179 Int_t i;
180 for (i = 0; i < n; i++)
181 z(i) = ( y(i) >= 0.0 ? 1.0 : -1.0 );
182 if (!TransSolve(z))
183 return iter;
184 Int_t imax = 0;
185 Double_t maxz = TMath::Abs(z(0));
186 for (i = 1; i < n; i++) {
187 const Double_t absz = TMath::Abs(z(i));
188 if ( absz > maxz ) {
189 maxz = absz;
190 imax = i;
191 }
192 }
193 stop = (maxz <= b*z);
194 if (!stop) {
195 b = 0.0;
196 b(imax) = 1.0;
197 }
198 }
199 iter--;
200 } while (!stop && iter);
201 est = inv_norm1;
202
203 return iter;
204}
205
206////////////////////////////////////////////////////////////////////////////////
207
209{
210// Returns product of matrix diagonal elements in d1 and d2. d1 is a mantissa and d2
211// an exponent for powers of 2. If matrix is in diagonal or triangular-matrix form this
212// will be the determinant.
213// Based on Bowler, Martin, Peters and Wilkinson in HACLA
214
215 const Double_t zero = 0.0;
216 const Double_t one = 1.0;
217 const Double_t four = 4.0;
218 const Double_t sixteen = 16.0;
219 const Double_t sixteenth = 0.0625;
220
221 const Int_t n = diag.GetNrows();
222
223 Double_t t1 = 1.0;
224 Double_t t2 = 0.0;
225 Int_t niter2 =0;
226 Int_t niter3 =0;
227 for (Int_t i = 0; (((i < n) && (t1 !=zero ))); i++) {
228 if (TMath::Abs(diag(i)) > tol) {
229 t1 *= (Double_t) diag(i);
230 while ( TMath::Abs(t1) >= one) {
231 t1 *= sixteenth;
232 t2 += four;
233 niter2++;
234 if ( niter2>100) break;
235 }
236 while ( TMath::Abs(t1) < sixteenth) {
237 t1 *= sixteen;
238 t2 -= four;
239 niter3++;
240 if (niter3>100) break;
241 }
242 } else {
243 t1 = zero;
244 t2 = zero;
245 }
246 }
247 d1 = t1;
248 d2 = t2;
249
250 return;
251}
252
253////////////////////////////////////////////////////////////////////////////////
254/// Matrix condition number
255
257{
258 if ( !TestBit(kCondition) ) {
259 fCondition = -1;
260 if (TestBit(kSingular))
261 return fCondition;
262 if ( !TestBit(kDecomposed) ) {
263 if (!Decompose())
264 return fCondition;
265 }
267 if (Hager(invNorm))
269 else // no convergence in Hager
270 Error("Condition()","Hager procedure did NOT converge");
272 }
273 return fCondition;
274}
275
276////////////////////////////////////////////////////////////////////////////////
277/// Solve set of equations with RHS in columns of B
278
280{
281 const TMatrixDBase &m = GetDecompMatrix();
282 R__ASSERT(m.IsValid() && B.IsValid());
283
284 const Int_t colLwb = B.GetColLwb();
285 const Int_t colUpb = B.GetColUpb();
286 Bool_t status = kTRUE;
287 for (Int_t icol = colLwb; icol <= colUpb && status; icol++) {
289 status &= Solve(b);
290 }
291
292 return status;
293}
294
295////////////////////////////////////////////////////////////////////////////////
296/// Matrix determinant det = d1*TMath::Power(2.,d2)
297
299{
300 if ( !TestBit(kDetermined) ) {
301 if ( !TestBit(kDecomposed) )
302 Decompose();
303 if (TestBit(kSingular) ) {
304 fDet1 = 0.0;
305 fDet2 = 0.0;
306 } else {
307 const TMatrixDBase &m = GetDecompMatrix();
308 R__ASSERT(m.IsValid());
309 TVectorD diagv(m.GetNrows());
310 for (Int_t i = 0; i < diagv.GetNrows(); i++)
311 diagv(i) = m(i,i);
313 }
315 }
316 d1 = fDet1;
317 d2 = fDet2;
318}
319
320////////////////////////////////////////////////////////////////////////////////
321/// Print class members
322
323void TDecompBase::Print(Option_t * /*opt*/) const
324{
325 printf("fTol = %.4e\n",fTol);
326 printf("fDet1 = %.4e\n",fDet1);
327 printf("fDet2 = %.4e\n",fDet2);
328 printf("fCondition = %.4e\n",fCondition);
329 printf("fRowLwb = %d\n",fRowLwb);
330 printf("fColLwb = %d\n",fColLwb);
331}
332
333////////////////////////////////////////////////////////////////////////////////
334/// Assignment operator
335
337{
338 if (this != &source) {
340 fTol = source.fTol;
341 fDet1 = source.fDet1;
342 fDet2 = source.fDet2;
343 fCondition = source.fCondition;
344 fRowLwb = source.fRowLwb;
345 fColLwb = source.fColLwb;
346 }
347 return *this;
348}
349
350////////////////////////////////////////////////////////////////////////////////
351/// Define a Householder-transformation through the parameters up and b .
352
355{
356 const Int_t n = vc.GetNrows();
357 const Double_t * const vp = vc.GetMatrixArray();
358
360 Int_t i;
361 for (i = l; i < n; i++)
362 c = TMath::Max(TMath::Abs(vp[i]),c);
363
364 up = 0.0;
365 beta = 0.0;
366 if (c <= tol) {
367// Warning("DefHouseHolder","max vector=%.4e < %.4e",c,tol);
368 return kFALSE;
369 }
370
371 Double_t sd = vp[lp]/c; sd *= sd;
372 for (i = l; i < n; i++) {
373 const Double_t tmp = vp[i]/c;
374 sd += tmp*tmp;
375 }
376
378 if (vp[lp] > 0.) vpprim = -vpprim;
379 up = vp[lp]-vpprim;
380 beta = 1./(vpprim*up);
381
382 return kTRUE;
383}
384
385////////////////////////////////////////////////////////////////////////////////
386/// Apply Householder-transformation.
387
390{
391 const Int_t nv = vc.GetNrows();
392 const Int_t nc = (cr.GetMatrix())->GetNcols();
393
394 if (nv > nc) {
395 Error("ApplyHouseHolder(const TVectorD &,..,TMatrixDRow &)","matrix row too short");
396 return;
397 }
398
399 const Int_t inc_c = cr.GetInc();
400 const Double_t * const vp = vc.GetMatrixArray();
401 Double_t * cp = cr.GetPtr();
402
403 Double_t s = cp[lp*inc_c]*up;
404 Int_t i;
405 for (i = l; i < nv; i++)
406 s += cp[i*inc_c]*vp[i];
407
408 s = s*beta;
409 cp[lp*inc_c] += s*up;
410 for (i = l; i < nv; i++)
411 cp[i*inc_c] += s*vp[i];
412}
413
414////////////////////////////////////////////////////////////////////////////////
415/// Apply Householder-transformation.
416
419{
420 const Int_t nv = vc.GetNrows();
421 const Int_t nc = (cc.GetMatrix())->GetNrows();
422
423 if (nv > nc) {
424 Error("ApplyHouseHolder(const TVectorD &,..,TMatrixDRow &)","matrix column too short");
425 return;
426 }
427
428 const Int_t inc_c = cc.GetInc();
429 const Double_t * const vp = vc.GetMatrixArray();
430 Double_t * cp = cc.GetPtr();
431
432 Double_t s = cp[lp*inc_c]*up;
433 Int_t i;
434 for (i = l; i < nv; i++)
435 s += cp[i*inc_c]*vp[i];
436
437 s = s*beta;
438 cp[lp*inc_c] += s*up;
439 for (i = l; i < nv; i++)
440 cp[i*inc_c] += s*vp[i];
441}
442
443////////////////////////////////////////////////////////////////////////////////
444/// Apply Householder-transformation.
445
448{
449 const Int_t nv = vc.GetNrows();
450 const Int_t nc = cv.GetNrows();
451
452 if (nv > nc) {
453 Error("ApplyHouseHolder(const TVectorD &,..,TVectorD &)","vector too short");
454 return;
455 }
456
457 const Double_t * const vp = vc.GetMatrixArray();
458 Double_t * cp = cv.GetMatrixArray();
459
460 Double_t s = cp[lp]*up;
461 Int_t i;
462 for (i = l; i < nv; i++)
463 s += cp[i]*vp[i];
464
465 s = s*beta;
466 cp[lp] += s*up;
467 for (i = l; i < nv; i++)
468 cp[i] += s*vp[i];
469}
470
471////////////////////////////////////////////////////////////////////////////////
472/// Defines a Givens-rotation by calculating 2 rotation parameters c and s.
473/// The rotation is defined with the vector components v1 and v2.
474
476{
477 const Double_t a1 = TMath::Abs(v1);
478 const Double_t a2 = TMath::Abs(v2);
479 if (a1 > a2) {
480 const Double_t w = v2/v1;
481 const Double_t q = TMath::Hypot(1.,w);
482 c = 1./q;
483 if (v1 < 0.) c = -c;
484 s = c*w;
485 } else {
486 if (v2 != 0) {
487 const Double_t w = v1/v2;
488 const Double_t q = TMath::Hypot(1.,w);
489 s = 1./q;
490 if (v2 < 0.) s = -s;
491 c = s*w;
492 } else {
493 c = 1.;
494 s = 0.;
495 }
496 }
497}
498
499////////////////////////////////////////////////////////////////////////////////
500/// Define and apply a Givens-rotation by calculating 2 rotation
501/// parameters c and s. The rotation is defined with and applied to the vector
502/// components v1 and v2.
503
505{
506 const Double_t a1 = TMath::Abs(v1);
507 const Double_t a2 = TMath::Abs(v2);
508 if (a1 > a2) {
509 const Double_t w = v2/v1;
510 const Double_t q = TMath::Hypot(1.,w);
511 c = 1./q;
512 if (v1 < 0.) c = -c;
513 s = c*w;
514 v1 = a1*q;
515 v2 = 0.;
516 } else {
517 if (v2 != 0) {
518 const Double_t w = v1/v2;
519 const Double_t q = TMath::Hypot(1.,w);
520 s = 1./q;
521 if (v2 < 0.) s = -s;
522 c = s*w;
523 v1 = a2*q;
524 v2 = 0.;
525 } else {
526 c = 1.;
527 s = 0.;
528 }
529 }
530}
531
532////////////////////////////////////////////////////////////////////////////////
533/// Apply a Givens transformation as defined by c and s to the vector components
534/// v1 and v2 .
535
537{
538 const Double_t w = z1*c+z2*s;
539 z2 = -z1*s+z2*c;
540 z1 = w;
541}
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
bool Bool_t
Boolean (0=false, 1=true) (bool)
Definition RtypesCore.h:77
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
float * q
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.
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
void Print(Option_t *opt="") const override
Print class members.
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
Int_t GetColLwb() const
Int_t GetColUpb() const
Bool_t IsValid() const
Mother of all ROOT objects.
Definition TObject.h:41
TObject & operator=(const TObject &rhs) noexcept
TObject assignment operator.
Definition TObject.h:299
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:202
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:864
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
Double_t y[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Double_t Hypot(Double_t x, Double_t y)
Returns sqrt(x*x + y*y)
Definition TMath.cxx:59
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
auto * t1
Definition textangle.C:20