Logo ROOT   6.12/07
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 
143 TDecompBase::TDecompBase(const TDecompBase &another) : TObject(another)
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 
324 void 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 
419  Int_t lp,Int_t l,TMatrixDColumn &cc)
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 }
static double B[]
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
auto * m
Definition: textangle.C:8
virtual Bool_t TransSolve(TVectorD &b)=0
Int_t GetInc() const
Int_t Hager(Double_t &est, Int_t iter=5)
Double_t fCondition
Definition: TDecompBase.h:39
Int_t fRowLwb
Definition: TDecompBase.h:40
const char Option_t
Definition: RtypesCore.h:62
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.
void ApplyHouseHolder(const TVectorD &vc, Double_t up, Double_t beta, Int_t lp, Int_t l, TMatrixDRow &cr)
Apply Householder-transformation.
virtual Bool_t Decompose()=0
Decomposition Base class.
Definition: TDecompBase.h:33
TDecompBase()
Default constructor.
#define R__ASSERT(e)
Definition: TError.h:96
virtual const TMatrixDBase & GetDecompMatrix() const =0
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
Int_t GetNrows() const
Definition: TVectorT.h:75
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
double beta(double x, double y)
Calculates the beta function.
Double_t fTol
Definition: TDecompBase.h:36
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 .
Int_t GetColLwb() const
Definition: TMatrixTBase.h:123
Element * GetPtr() const
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Int_t fColLwb
Definition: TDecompBase.h:41
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:271
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
Element * GetMatrixArray()
Definition: TVectorT.h:78
virtual Int_t GetNcols() const =0
virtual Bool_t Solve(TVectorD &b)=0
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
void Print(Option_t *opt="") const
Print class members.
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.
Element Norm1() const
Compute the 1-norm of the vector SUM{ |v[i]| }.
Definition: TVectorT.cxx:565
Int_t GetColUpb() const
Definition: TMatrixTBase.h:124
virtual Int_t GetNrows() const =0
Linear Algebra Package.
const TMatrixTBase< Element > * GetMatrix() const
REAL epsilon
Definition: triangle.c:617
Double_t fDet2
Definition: TDecompBase.h:38
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t fDet1
Definition: TDecompBase.h:37
Element * GetPtr() const
auto * t1
Definition: textangle.C:20
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
static constexpr double s
Double_t Hypot(Double_t x, Double_t y)
Definition: TMath.cxx:60
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
Mother of all ROOT objects.
Definition: TObject.h:37
const TMatrixTBase< Element > * GetMatrix() const
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
auto * l
Definition: textangle.C:4
virtual Double_t Condition()
Matrix condition number.
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
float * q
Definition: THbookFile.cxx:87
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
Bool_t IsValid() const
Definition: TMatrixTBase.h:145
Int_t GetInc() const