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