Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 the 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
76TDecompQRH::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++) {
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
582{
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) {
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}
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:382
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 ApplyHouseHolder(const TVectorD &vc, Double_t up, Double_t b, Int_t lp, Int_t l, TMatrixDRow &cr)
Apply Householder-transformation.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float * q
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
TMatrixTColumn< Double_t > TMatrixDColumn
Decomposition Base class.
Definition TDecompBase.h:34
Int_t GetRowLwb() const
Definition TDecompBase.h:73
Double_t fDet1
Definition TDecompBase.h:37
Double_t fDet2
Definition TDecompBase.h:38
void ResetStatus()
Definition TDecompBase.h:43
Int_t GetColLwb() const
Definition TDecompBase.h:74
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
Int_t fRowLwb
Definition TDecompBase.h:40
Double_t fTol
Definition TDecompBase.h:36
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Double_t fCondition
Definition TDecompBase.h:39
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)
QR Decomposition class.
Definition TDecompQRH.h:26
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
TMatrixD Invert()
Definition TDecompQRH.h:77
Bool_t Solve(TVectorD &b) override
Solve Ax=b assuming the QR form of A is stored in fR,fQ and fW, but assume b has not been transformed...
Int_t GetNrows() const override
Definition TDecompQRH.h:50
void Print(Option_t *opt="") const override
Print the class members.
void Det(Double_t &d1, Double_t &d2) override
This routine calculates the absolute (!) value of the determinant |det| = d1*TMath::Power(2....
TDecompQRH & operator=(const TDecompQRH &source)
Assignment operator.
TMatrixD GetOrthogonalMatrix() const
For a matrix A(m,n), return the OtrhogonalMatrix Q such as A = Q * R.
TMatrixD fQ
Definition TDecompQRH.h:30
TMatrixD fR
Definition TDecompQRH.h:31
static Bool_t QRH(TMatrixD &q, TVectorD &diagR, TVectorD &up, TVectorD &w, Double_t tol)
Decomposition function .
Bool_t Decompose() override
QR decomposition of matrix a by Householder transformations, see Golub & Loan first edition p41 & Sec...
TVectorD fW
Definition TDecompQRH.h:33
Int_t GetNcols() const override
Definition TDecompQRH.h:51
TVectorD fUp
Definition TDecompQRH.h:32
Bool_t TransSolve(TVectorD &b) override
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...
Int_t GetNrows() const
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Int_t GetColLwb() const
Int_t GetNcols() const
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.
const TMatrixTBase< Element > * GetMatrix() const
Element * GetPtr() const
const Element * GetMatrixArray() const override
Definition TMatrixT.h:228
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:199
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:786
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:993
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:294
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition TVectorT.cxx:349
void Print(Option_t *option="") const override
Print the vector as a list of elements.
const Int_t n
Definition legend1.C:16
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition rsaaux.cxx:949