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