Logo ROOT   6.21/01
Reference Guide
TDecompQRH.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 TDecompQRH
13  \ingroup Matrix
14 
15  QR Decomposition class
16 
17  Decompose a general (m x n) matrix A into A = fQ' fR H where
18 
19 ~~~
20  fQ : (m x n) - internal Q' matrix (not orthoginal)
21  fR : (n x n) - upper triangular matrix
22  H : HouseHolder matrix which is stored through
23  fUp: (n) - vector with Householder up's
24  fW : (n) - vector with Householder beta's
25 ~~~
26 
27  If row/column index of A starts at (rowLwb,colLwb) then
28  the decomposed matrices start from :
29 ~~~
30  fQ' : (rowLwb,0)
31  fR : (0,colLwb)
32  and the decomposed vectors start from :
33  fUp : (0)
34  fW : (0)
35 ~~~
36 
37  In order to get thw QR dcomposition of A (i.e. A = QR )
38  The orthoginal matrix Q needs to be computed from the internal Q' and
39  the up's and beta's vector defining the Householder transformation
40 
41  The orthogonal Q matrix is returned to the user by calling the
42  function TDecompQRH::GetOrthogonalMatrix()
43 
44  Errors arise from formation of reflectors i.e. singularity .
45  Note it attempts to handle the cases where the nRow <= nCol .
46 */
47 
48 #include "TDecompQRH.h"
49 #include "TError.h" // For R__ASSERT
51 
52 ////////////////////////////////////////////////////////////////////////////////
53 /// Constructor for (nrows x ncols) matrix
54 
56 {
57  if (nrows < ncols) {
58  Error("TDecompQRH(Int_t,Int_t","matrix rows should be >= columns");
59  return;
60  }
61 
62  fQ.ResizeTo(nrows,ncols);
63  fR.ResizeTo(ncols,ncols);
64  if (nrows <= ncols) {
65  fW.ResizeTo(nrows);
66  fUp.ResizeTo(nrows);
67  } else {
68  fW.ResizeTo(ncols);
69  fUp.ResizeTo(ncols);
70  }
71 }
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Constructor for ([row_lwb..row_upb] x [col_lwb..col_upb]) matrix
75 
76 TDecompQRH::TDecompQRH(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb)
77 {
78  const Int_t nrows = row_upb-row_lwb+1;
79  const Int_t ncols = col_upb-col_lwb+1;
80 
81  if (nrows < ncols) {
82  Error("TDecompQRH(Int_t,Int_t,Int_t,Int_t","matrix rows should be >= columns");
83  return;
84  }
85 
86  fRowLwb = row_lwb;
87  fColLwb = col_lwb;
88 
89  fQ.ResizeTo(nrows,ncols);
90  fR.ResizeTo(ncols,ncols);
91  if (nrows <= ncols) {
92  fW.ResizeTo(nrows);
93  fUp.ResizeTo(nrows);
94  } else {
95  fW.ResizeTo(ncols);
96  fUp.ResizeTo(ncols);
97  }
98 }
99 
100 ////////////////////////////////////////////////////////////////////////////////
101 /// Constructor for general matrix A .
102 
104 {
105  R__ASSERT(a.IsValid());
106  if (a.GetNrows() < a.GetNcols()) {
107  Error("TDecompQRH(const TMatrixD &","matrix rows should be >= columns");
108  return;
109  }
110 
112  fCondition = a.Norm1();
113  fTol = a.GetTol();
114  if (tol > 0.0)
115  fTol = tol;
116 
117  fRowLwb = a.GetRowLwb();
118  fColLwb = a.GetColLwb();
119  const Int_t nRow = a.GetNrows();
120  const Int_t nCol = a.GetNcols();
121 
122  fQ.ResizeTo(nRow,nCol);
123  memcpy(fQ.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
124  fR.ResizeTo(nCol,nCol);
125  if (nRow <= nCol) {
126  fW.ResizeTo(nRow);
127  fUp.ResizeTo(nRow);
128  } else {
129  fW.ResizeTo(nCol);
130  fUp.ResizeTo(nCol);
131  }
132 }
133 
134 ////////////////////////////////////////////////////////////////////////////////
135 /// Copy constructor
136 
138 {
139  *this = another;
140 }
141 
142 ////////////////////////////////////////////////////////////////////////////////
143 /// QR decomposition of matrix a by Householder transformations,
144 /// see Golub & Loan first edition p41 & Sec 6.2.
145 /// First fR is returned in upper triang of fQ and diagR. fQ returned in
146 /// 'u-form' in lower triang of fQ and fW, the latter containing the
147 /// "Householder betas".
148 /// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
149 
151 {
152  if (TestBit(kDecomposed)) return kTRUE;
153 
154  if ( !TestBit(kMatrixSet) ) {
155  Error("Decompose()","Matrix has not been set");
156  return kFALSE;
157  }
158 
159  const Int_t nRow = this->GetNrows();
160  const Int_t nCol = this->GetNcols();
161  const Int_t rowLwb = this->GetRowLwb();
162  const Int_t colLwb = this->GetColLwb();
163 
164  TVectorD diagR;
165  Double_t work[kWorkMax];
166  if (nCol > kWorkMax) diagR.ResizeTo(nCol);
167  else diagR.Use(nCol,work);
168 
169  if (QRH(fQ,diagR,fUp,fW,fTol)) {
170  for (Int_t i = 0; i < nRow; i++) {
171  const Int_t ic = (i < nCol) ? i : nCol;
172  for (Int_t j = ic ; j < nCol; j++)
173  fR(i,j) = fQ(i,j);
174  }
175  TMatrixDDiag diag(fR); diag = diagR;
176 
177  fQ.Shift(rowLwb,0);
178  fR.Shift(0,colLwb);
179 
181  }
182 
183  return kTRUE;
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// Decomposition function .
188 
190 {
191  const Int_t nRow = q.GetNrows();
192  const Int_t nCol = q.GetNcols();
193 
194  const Int_t n = (nRow <= nCol) ? nRow-1 : nCol;
195 
196  for (Int_t k = 0 ; k < n ; k++) {
197  const TVectorD qc_k = TMatrixDColumn_const(q,k);
198  if (!DefHouseHolder(qc_k,k,k+1,up(k),w(k),tol))
199  return kFALSE;
200  diagR(k) = qc_k(k)-up(k);
201  if (k < nCol-1) {
202  // Apply HouseHolder to sub-matrix
203  for (Int_t j = k+1; j < nCol; j++) {
204  TMatrixDColumn qc_j = TMatrixDColumn(q,j);
205  ApplyHouseHolder(qc_k,up(k),w(k),k,k+1,qc_j);
206  }
207  }
208  }
209 
210  if (nRow <= nCol) {
211  diagR(nRow-1) = q(nRow-1,nRow-1);
212  up(nRow-1) = 0;
213  w(nRow-1) = 0;
214  }
215 
216  return kTRUE;
217 }
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// Set matrix to be decomposed
221 
223 {
224  R__ASSERT(a.IsValid());
225 
226  ResetStatus();
227  if (a.GetNrows() < a.GetNcols()) {
228  Error("TDecompQRH(const TMatrixD &","matrix rows should be >= columns");
229  return;
230  }
231 
233  fCondition = a.Norm1();
234 
235  fRowLwb = a.GetRowLwb();
236  fColLwb = a.GetColLwb();
237  const Int_t nRow = a.GetNrows();
238  const Int_t nCol = a.GetNcols();
239 
240  fQ.ResizeTo(nRow,nCol);
241  memcpy(fQ.GetMatrixArray(),a.GetMatrixArray(),nRow*nCol*sizeof(Double_t));
242  fR.ResizeTo(nCol,nCol);
243  if (nRow <= nCol) {
244  fW.ResizeTo(nRow);
245  fUp.ResizeTo(nRow);
246  } else {
247  fW.ResizeTo(nCol);
248  fUp.ResizeTo(nCol);
249  }
250 }
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 /// Solve Ax=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
254 /// has *not* been transformed. Solution returned in b.
255 
257 {
258  R__ASSERT(b.IsValid());
259  if (TestBit(kSingular)) {
260  Error("Solve()","Matrix is singular");
261  return kFALSE;
262  }
263  if ( !TestBit(kDecomposed) ) {
264  if (!Decompose()) {
265  Error("Solve()","Decomposition failed");
266  return kFALSE;
267  }
268  }
269 
270  if (fQ.GetNrows() != b.GetNrows() || fQ.GetRowLwb() != b.GetLwb()) {
271  Error("Solve(TVectorD &","vector and matrix incompatible");
272  return kFALSE;
273  }
274 
275  const Int_t nQRow = fQ.GetNrows();
276  const Int_t nQCol = fQ.GetNcols();
277 
278  // Calculate Q^T.b
279  const Int_t nQ = (nQRow <= nQCol) ? nQRow-1 : nQCol;
280  for (Int_t k = 0; k < nQ; k++) {
281  const TVectorD qc_k = TMatrixDColumn_const(fQ,k);
282  ApplyHouseHolder(qc_k,fUp(k),fW(k),k,k+1,b);
283  }
284 
285  const Int_t nRCol = fR.GetNcols();
286 
287  const Double_t *pR = fR.GetMatrixArray();
288  Double_t *pb = b.GetMatrixArray();
289 
290  // Backward substitution
291  for (Int_t i = nRCol-1; i >= 0; i--) {
292  const Int_t off_i = i*nRCol;
293  Double_t r = pb[i];
294  for (Int_t j = i+1; j < nRCol; j++)
295  r -= pR[off_i+j]*pb[j];
296  if (TMath::Abs(pR[off_i+i]) < fTol)
297  {
298  Error("Solve(TVectorD &)","R[%d,%d]=%.4e < %.4e",i,i,pR[off_i+i],fTol);
299  return kFALSE;
300  }
301  pb[i] = r/pR[off_i+i];
302  }
303 
304  return kTRUE;
305 }
306 
307 ////////////////////////////////////////////////////////////////////////////////
308 /// Solve Ax=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
309 /// has *not* been transformed. Solution returned in b.
310 
312 {
313  TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
314  R__ASSERT(b->IsValid());
315  if (TestBit(kSingular)) {
316  Error("Solve()","Matrix is singular");
317  return kFALSE;
318  }
319  if ( !TestBit(kDecomposed) ) {
320  if (!Decompose()) {
321  Error("Solve()","Decomposition failed");
322  return kFALSE;
323  }
324  }
325 
326  if (fQ.GetNrows() != b->GetNrows() || fQ.GetRowLwb() != b->GetRowLwb())
327  {
328  Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
329  return kFALSE;
330  }
331 
332  const Int_t nQRow = fQ.GetNrows();
333  const Int_t nQCol = fQ.GetNcols();
334 
335  // Calculate Q^T.b
336  const Int_t nQ = (nQRow <= nQCol) ? nQRow-1 : nQCol;
337  for (Int_t k = 0; k < nQ; k++) {
338  const TVectorD qc_k = TMatrixDColumn_const(fQ,k);
339  ApplyHouseHolder(qc_k,fUp(k),fW(k),k,k+1,cb);
340  }
341 
342  const Int_t nRCol = fR.GetNcols();
343 
344  const Double_t *pR = fR.GetMatrixArray();
345  Double_t *pcb = cb.GetPtr();
346  const Int_t inc = cb.GetInc();
347 
348  // Backward substitution
349  for (Int_t i = nRCol-1; i >= 0; i--) {
350  const Int_t off_i = i*nRCol;
351  const Int_t off_i2 = i*inc;
352  Double_t r = pcb[off_i2];
353  for (Int_t j = i+1; j < nRCol; j++)
354  r -= pR[off_i+j]*pcb[j*inc];
355  if (TMath::Abs(pR[off_i+i]) < fTol)
356  {
357  Error("Solve(TMatrixDColumn &)","R[%d,%d]=%.4e < %.4e",i,i,pR[off_i+i],fTol);
358  return kFALSE;
359  }
360  pcb[off_i2] = r/pR[off_i+i];
361  }
362 
363  return kTRUE;
364 }
365 
366 ////////////////////////////////////////////////////////////////////////////////
367 /// Solve A^T x=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
368 /// has *not* been transformed. Solution returned in b.
369 
371 {
372  R__ASSERT(b.IsValid());
373  if (TestBit(kSingular)) {
374  Error("TransSolve()","Matrix is singular");
375  return kFALSE;
376  }
377  if ( !TestBit(kDecomposed) ) {
378  if (!Decompose()) {
379  Error("TransSolve()","Decomposition failed");
380  return kFALSE;
381  }
382  }
383 
384  if (fQ.GetNrows() != fQ.GetNcols() || fQ.GetRowLwb() != fQ.GetColLwb()) {
385  Error("TransSolve(TVectorD &","matrix should be square");
386  return kFALSE;
387  }
388 
389  if (fR.GetNrows() != b.GetNrows() || fR.GetRowLwb() != b.GetLwb()) {
390  Error("TransSolve(TVectorD &","vector and matrix incompatible");
391  return kFALSE;
392  }
393 
394  const Double_t *pR = fR.GetMatrixArray();
395  Double_t *pb = b.GetMatrixArray();
396 
397  const Int_t nRCol = fR.GetNcols();
398 
399  // Backward substitution
400  for (Int_t i = 0; i < nRCol; i++) {
401  const Int_t off_i = i*nRCol;
402  Double_t r = pb[i];
403  for (Int_t j = 0; j < i; j++) {
404  const Int_t off_j = j*nRCol;
405  r -= pR[off_j+i]*pb[j];
406  }
407  if (TMath::Abs(pR[off_i+i]) < fTol)
408  {
409  Error("TransSolve(TVectorD &)","R[%d,%d]=%.4e < %.4e",i,i,pR[off_i+i],fTol);
410  return kFALSE;
411  }
412  pb[i] = r/pR[off_i+i];
413  }
414 
415  const Int_t nQRow = fQ.GetNrows();
416 
417  // Calculate Q.b; it was checked nQRow == nQCol
418  for (Int_t k = nQRow-1; k >= 0; k--) {
419  const TVectorD qc_k = TMatrixDColumn_const(fQ,k);
420  ApplyHouseHolder(qc_k,fUp(k),fW(k),k,k+1,b);
421  }
422 
423  return kTRUE;
424 }
425 
426 ////////////////////////////////////////////////////////////////////////////////
427 /// Solve A^T x=b assuming the QR form of A is stored in fR,fQ and fW, but assume b
428 /// has *not* been transformed. Solution returned in b.
429 
431 {
432  TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
433  R__ASSERT(b->IsValid());
434  if (TestBit(kSingular)) {
435  Error("TransSolve()","Matrix is singular");
436  return kFALSE;
437  }
438  if ( !TestBit(kDecomposed) ) {
439  if (!Decompose()) {
440  Error("TransSolve()","Decomposition failed");
441  return kFALSE;
442  }
443  }
444 
445  if (fQ.GetNrows() != fQ.GetNcols() || fQ.GetRowLwb() != fQ.GetColLwb()) {
446  Error("TransSolve(TMatrixDColumn &","matrix should be square");
447  return kFALSE;
448  }
449 
450  if (fR.GetNrows() != b->GetNrows() || fR.GetRowLwb() != b->GetRowLwb()) {
451  Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
452  return kFALSE;
453  }
454 
455  const Double_t *pR = fR.GetMatrixArray();
456  Double_t *pcb = cb.GetPtr();
457  const Int_t inc = cb.GetInc();
458 
459  const Int_t nRCol = fR.GetNcols();
460 
461  // Backward substitution
462  for (Int_t i = 0; i < nRCol; i++) {
463  const Int_t off_i = i*nRCol;
464  const Int_t off_i2 = i*inc;
465  Double_t r = pcb[off_i2];
466  for (Int_t j = 0; j < i; j++) {
467  const Int_t off_j = j*nRCol;
468  r -= pR[off_j+i]*pcb[j*inc];
469  }
470  if (TMath::Abs(pR[off_i+i]) < fTol)
471  {
472  Error("TransSolve(TMatrixDColumn &)","R[%d,%d]=%.4e < %.4e",i,i,pR[off_i+i],fTol);
473  return kFALSE;
474  }
475  pcb[off_i2] = r/pR[off_i+i];
476  }
477 
478  const Int_t nQRow = fQ.GetNrows();
479 
480  // Calculate Q.b; it was checked nQRow == nQCol
481  for (Int_t k = nQRow-1; k >= 0; k--) {
482  const TVectorD qc_k = TMatrixDColumn_const(fQ,k);
483  ApplyHouseHolder(qc_k,fUp(k),fW(k),k,k+1,cb);
484  }
485 
486  return kTRUE;
487 }
488 
489 ////////////////////////////////////////////////////////////////////////////////
490 /// This routine calculates the absolute (!) value of the determinant
491 /// |det| = d1*TMath::Power(2.,d2)
492 
494 {
495  if ( !TestBit(kDetermined) ) {
496  if ( !TestBit(kDecomposed) )
497  Decompose();
498  if (TestBit(kSingular)) {
499  fDet1 = 0.0;
500  fDet2 = 0.0;
501  } else
502  TDecompBase::Det(d1,d2);
504  }
505  d1 = fDet1;
506  d2 = fDet2;
507 }
508 
509 ////////////////////////////////////////////////////////////////////////////////
510 /// For a matrix A(m,n), return the OtrhogonalMatrix Q such as
511 /// A = Q * R
512 ///
513 /// Note that this Q is not th einternal fQ matrix obtained in the QRH decomposition, but can be computed
514 /// from the fQ and the up and beta vector's defining the Householder transformation
515 
517 {
518  // apply HouseHolder transformation starting from the identity
519  // Calculate Q.b; it was checked nQRow == nQCol
520 
521  const Int_t nRow = this->GetNrows();
522  const Int_t nCol = this->GetNcols();
523  // remmber nCol <= nRow
524  TMatrixD orthogQ(nRow, nCol);
525  // start from identity matrix
526  for (int i = 0; i < nCol; ++i)
527  orthogQ(i, i) = 1;
528 
529 
530  // apply the HouseHolder transformations for each column of Q
531  for (int j = 0; j < nCol; ++j) {
532  TMatrixDColumn b = TMatrixDColumn(orthogQ, j);
533  int nQRow = fQ.GetNrows();
534  for (Int_t k = nQRow - 1; k >= 0; k--) {
535  const TVectorD qc_k = TMatrixDColumn_const(fQ, k);
536  ApplyHouseHolder(qc_k, fUp(k), fW(k), k, k + 1, b);
537  }
538  }
539  return orthogQ;
540 }
541 ////////////////////////////////////////////////////////////////////////////////
542 /// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
543 /// The user should always supply a matrix of size (m x m) !
544 /// If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
545 /// should be used .
546 
548 {
549  if (inv.GetNrows() != GetNrows() || inv.GetNcols() != GetNrows() ||
550  inv.GetRowLwb() != GetRowLwb() || inv.GetColLwb() != GetColLwb()) {
551  Error("Invert(TMatrixD &","Input matrix has wrong shape");
552  return kFALSE;
553  }
554 
555  inv.UnitMatrix();
556  const Bool_t status = MultiSolve(inv);
557 
558  return status;
559 }
560 
561 ////////////////////////////////////////////////////////////////////////////////
562 /// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
563 /// (n x m) Ainv is returned .
564 
566 {
567  const Int_t rowLwb = GetRowLwb();
568  const Int_t colLwb = GetColLwb();
569  const Int_t rowUpb = rowLwb+GetNrows()-1;
570  TMatrixD inv(rowLwb,rowUpb,colLwb,colLwb+GetNrows()-1);
571  inv.UnitMatrix();
572  status = MultiSolve(inv);
573  inv.ResizeTo(rowLwb,rowLwb+GetNcols()-1,colLwb,colLwb+GetNrows()-1);
574 
575  return inv;
576 }
577 
578 ////////////////////////////////////////////////////////////////////////////////
579 /// Print the class members
580 
581 void TDecompQRH::Print(Option_t *opt) const
582 {
583  TDecompBase::Print(opt);
584  fQ.Print("fQ");
585  fR.Print("fR");
586  fUp.Print("fUp");
587  fW.Print("fW");
588 }
589 
590 ////////////////////////////////////////////////////////////////////////////////
591 /// Assignment operator
592 
594 {
595  if (this != &source) {
596  TDecompBase::operator=(source);
597  fQ.ResizeTo(source.fQ);
598  fR.ResizeTo(source.fR);
599  fUp.ResizeTo(source.fUp);
600  fW.ResizeTo(source.fW);
601  fQ = source.fQ;
602  fR = source.fR;
603  fUp = source.fUp;
604  fW = source.fW;
605  }
606  return *this;
607 }
TMatrixD GetOrthogonalMatrix() const
For a matrix A(m,n), return the OtrhogonalMatrix Q such as A = Q * R.
Definition: TDecompQRH.cxx:516
QR Decomposition class.
Definition: TDecompQRH.h:25
virtual TMatrixTBase< Element > & Shift(Int_t row_shift, Int_t col_shift)
Shift the row index by adding row_shift and the column index by adding col_shift, respectively...
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
TMatrixD fQ
Definition: TDecompQRH.h:30
TMatrixD fR
Definition: TDecompQRH.h:31
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
Int_t GetInc() const
Double_t fCondition
Definition: TDecompBase.h:39
void Print(Option_t *opt="") const
Print the class members.
Definition: TDecompQRH.cxx:581
Int_t GetRowLwb() const
Definition: TDecompBase.h:73
Int_t fRowLwb
Definition: TDecompBase.h:40
const char Option_t
Definition: RtypesCore.h:62
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
TVectorD fUp
Definition: TDecompQRH.h:32
Decomposition Base class.
Definition: TDecompBase.h:33
#define R__ASSERT(e)
Definition: TError.h:96
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1361
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
Bool_t DefHouseHolder(const TVectorD &vc, Int_t lp, Int_t l, Double_t &up, Double_t &b, Double_t tol=0.0)
Define a Householder-transformation through the parameters up and b .
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual Int_t GetNrows() const
Definition: TDecompQRH.h:50
Double_t fTol
Definition: TDecompBase.h:36
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:347
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
Element * GetPtr() const
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
TDecompQRH & operator=(const TDecompQRH &source)
Assignment operator.
Definition: TDecompQRH.cxx:593
static Bool_t QRH(TMatrixD &q, TVectorD &diagR, TVectorD &up, TVectorD &w, Double_t tol)
Decomposition function .
Definition: TDecompQRH.cxx:189
Int_t fColLwb
Definition: TDecompBase.h:41
TVectorD fW
Definition: TDecompQRH.h:33
virtual Bool_t Decompose()
QR decomposition of matrix a by Householder transformations, see Golub & Loan first edition p41 & Sec...
Definition: TDecompQRH.cxx:150
virtual Bool_t Solve(TVectorD &b)
Solve Ax=b assuming the QR form of A is stored in fR,fQ and fW, but assume b has not been transformed...
Definition: TDecompQRH.cxx:256
ROOT::R::TRInterface & r
Definition: Object.C:4
TMatrixD Invert()
Definition: TDecompQRH.h:77
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
void ApplyHouseHolder(const TVectorD &vc, Double_t up, Double_t b, Int_t lp, Int_t l, TMatrixDRow &cr)
Apply Householder-transformation.
auto * a
Definition: textangle.C:12
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 ResetStatus()
Definition: TDecompBase.h:43
virtual void Det(Double_t &d1, Double_t &d2)
This routine calculates the absolute (!) value of the determinant |det| = d1*TMath::Power(2.,d2)
Definition: TDecompQRH.cxx:493
Linear Algebra Package.
Double_t fDet2
Definition: TDecompBase.h:38
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t fDet1
Definition: TDecompBase.h:37
#define ClassImp(name)
Definition: Rtypes.h:365
virtual Bool_t TransSolve(TVectorD &b)
Solve A^T x=b assuming the QR form of A is stored in fR,fQ and fW, but assume b has not been transfor...
Definition: TDecompQRH.cxx:370
double Double_t
Definition: RtypesCore.h:55
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual Int_t GetNcols() const
Definition: TDecompQRH.h:51
const TMatrixTBase< Element > * GetMatrix() const
TMatrixTColumn< Double_t > TMatrixDColumn
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
Definition: TDecompQRH.cxx:222
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
float * q
Definition: THbookFile.cxx:87
const Bool_t kTRUE
Definition: RtypesCore.h:87
Int_t GetColLwb() const
Definition: TDecompBase.h:74
const Int_t n
Definition: legend1.C:16