ROOT  6.06/09
Reference Guide
TDecompLU.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 #include "TDecompLU.h"
13 #include "TMath.h"
14 
16 
17 ///////////////////////////////////////////////////////////////////////////
18 // //
19 // LU Decomposition class //
20 // //
21 // Decompose a general n x n matrix A into P A = L U //
22 // //
23 // where P is a permutation matrix, L is unit lower triangular and U //
24 // is upper triangular. //
25 // L is stored in the strict lower triangular part of the matrix fLU. //
26 // The diagonal elements of L are unity and are not stored. //
27 // U is stored in the diagonal and upper triangular part of the matrix //
28 // fU. //
29 // P is stored in the index array fIndex : j = fIndex[i] indicates that //
30 // row j and row i should be swapped . //
31 // //
32 // fSign gives the sign of the permutation, (-1)^n, where n is the //
33 // number of interchanges in the permutation. //
34 // //
35 // fLU has the same indexing range as matrix A . //
36 // //
37 // The decomposition fails if a diagonal element of abs(fLU) is == 0, //
38 // The matrix fUL is made invalid . //
39 // //
40 ///////////////////////////////////////////////////////////////////////////
41 
42 ////////////////////////////////////////////////////////////////////////////////
43 /// Default constructor
44 
46 {
47  fSign = 0.0;
48  fNIndex = 0;
49  fIndex = 0;
50  fImplicitPivot = 0;
51 }
52 
53 ////////////////////////////////////////////////////////////////////////////////
54 /// Constructor for (nrows x nrows) matrix
55 
57 {
58  fSign = 1.0;
59  fNIndex = nrows;
60  fIndex = new Int_t[fNIndex];
61  memset(fIndex,0,fNIndex*sizeof(Int_t));
62  fImplicitPivot = 0;
63  fLU.ResizeTo(nrows,nrows);
64 }
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// Constructor for ([row_lwb..row_upb] x [row_lwb..row_upb]) matrix
68 
70 {
71  const Int_t nrows = row_upb-row_lwb+1;
72  fSign = 1.0;
73  fNIndex = nrows;
74  fIndex = new Int_t[fNIndex];
75  memset(fIndex,0,fNIndex*sizeof(Int_t));
76  fRowLwb = row_lwb;
77  fColLwb = row_lwb;
78  fImplicitPivot = 0;
79  fLU.ResizeTo(row_lwb,row_lwb+nrows-1,row_lwb,row_lwb+nrows-1);
80 }
81 
82 ////////////////////////////////////////////////////////////////////////////////
83 /// Constructor for matrix a
84 
86 {
87  R__ASSERT(a.IsValid());
88 
89  if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
90  Error("TDecompLU(const TMatrixD &","matrix should be square");
91  return;
92  }
93 
95  fCondition = a.Norm1();
96  fImplicitPivot = implicit;
97  fTol = a.GetTol();
98  if (tol > 0)
99  fTol = tol;
100 
101  fSign = 1.0;
102  fNIndex = a.GetNcols();
103  fIndex = new Int_t[fNIndex];
104  memset(fIndex,0,fNIndex*sizeof(Int_t));
105 
106  fRowLwb = a.GetRowLwb();
107  fColLwb = a.GetColLwb();
108  fLU.ResizeTo(a);
109  fLU = a;
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Copy constructor
114 
115 TDecompLU::TDecompLU(const TDecompLU &another) : TDecompBase(another)
116 {
117  fNIndex = 0;
118  fIndex = 0;
119  *this = another;
120 }
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 /// Matrix A is decomposed in components U and L so that P * A = U * L
124 /// If the decomposition succeeds, bit kDecomposed is set , otherwise kSingular
125 
127 {
128  if (TestBit(kDecomposed)) return kTRUE;
129 
130  if ( !TestBit(kMatrixSet) ) {
131  Error("Decompose()","Matrix has not been set");
132  return kFALSE;
133  }
134 
135  Int_t nrZeros = 0;
136  Bool_t ok;
137  if (fImplicitPivot)
138  ok = DecomposeLUCrout(fLU,fIndex,fSign,fTol,nrZeros);
139  else
140  ok = DecomposeLUGauss(fLU,fIndex,fSign,fTol,nrZeros);
141 
142  if (!ok) SetBit(kSingular);
143  else SetBit(kDecomposed);
144 
145  return ok;
146 }
147 
148 ////////////////////////////////////////////////////////////////////////////////
149 /// Reconstruct the original matrix using the decomposition parts
150 
152 {
153  if (TestBit(kSingular)) {
154  Error("GetMatrix()","Matrix is singular");
155  return TMatrixD();
156  }
157  if ( !TestBit(kDecomposed) ) {
158  if (!Decompose()) {
159  Error("GetMatrix()","Decomposition failed");
160  return TMatrixD();
161  }
162  }
163 
164  TMatrixD mL = fLU;
165  TMatrixD mU = fLU;
166  Double_t * const pU = mU.GetMatrixArray();
167  Double_t * const pL = mL.GetMatrixArray();
168  const Int_t n = fLU.GetNcols();
169  for (Int_t irow = 0; irow < n; irow++) {
170  const Int_t off_row = irow*n;
171  for (Int_t icol = 0; icol < n; icol++) {
172  if (icol > irow) pL[off_row+icol] = 0.;
173  else if (icol < irow) pU[off_row+icol] = 0.;
174  else pL[off_row+icol] = 1.;
175  }
176  }
177 
178  TMatrixD a = mL*mU;
179 
180  // swap rows
181 
182  Double_t * const pA = a.GetMatrixArray();
183  for (Int_t i = n-1; i >= 0; i--) {
184  const Int_t j = fIndex[i];
185  if (j != i) {
186  const Int_t off_j = j*n;
187  const Int_t off_i = i*n;
188  for (Int_t k = 0; k < n; k++) {
189  const Double_t tmp = pA[off_j+k];
190  pA[off_j+k] = pA[off_i+k];
191  pA[off_i+k] = tmp;
192  }
193  }
194  }
195 
196  return a;
197 }
198 
199 ////////////////////////////////////////////////////////////////////////////////
200 /// Set matrix to be decomposed
201 
203 {
204  R__ASSERT(a.IsValid());
205 
206  ResetStatus();
207  if (a.GetNrows() != a.GetNcols() || a.GetRowLwb() != a.GetColLwb()) {
208  Error("TDecompLU(const TMatrixD &","matrix should be square");
209  return;
210  }
211 
213  fCondition = a.Norm1();
214 
215  fSign = 1.0;
216  if (fNIndex != a.GetNcols()) {
217  fNIndex = a.GetNcols();
218  delete [] fIndex;
219  fIndex = new Int_t[fNIndex];
220  memset(fIndex,0,fNIndex*sizeof(Int_t));
221  }
222 
223  fRowLwb = a.GetRowLwb();
224  fColLwb = a.GetColLwb();
225  fLU.ResizeTo(a);
226  fLU = a;
227 }
228 
229 ////////////////////////////////////////////////////////////////////////////////
230 /// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
231 /// been transformed. Solution returned in b.
232 
234 {
235  R__ASSERT(b.IsValid());
236  if (TestBit(kSingular)) {
237  Error("Solve()","Matrix is singular");
238  return kFALSE;
239  }
240  if ( !TestBit(kDecomposed) ) {
241  if (!Decompose()) {
242  Error("Solve()","Decomposition failed");
243  return kFALSE;
244  }
245  }
246 
247  if (fLU.GetNrows() != b.GetNrows() || fLU.GetRowLwb() != b.GetLwb()) {
248  Error("Solve(TVectorD &","vector and matrix incompatible");
249  return kFALSE;
250  }
251 
252  const Int_t n = fLU.GetNrows();
253 
254  const Double_t *pLU = fLU.GetMatrixArray();
255  Double_t *pb = b.GetMatrixArray();
256 
257  Int_t i;
258 
259  // Check for zero diagonals
260  for (i = 0; i < n ; i++) {
261  const Int_t off_i = i*n;
262  if (TMath::Abs(pLU[off_i+i]) < fTol) {
263  Error("Solve(TVectorD &b)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
264  return kFALSE;
265  }
266  }
267 
268  // Transform b allowing for leading zeros
269  Int_t nonzero = -1;
270  for (i = 0; i < n; i++) {
271  const Int_t off_i = i*n;
272  const Int_t iperm = fIndex[i];
273  Double_t r = pb[iperm];
274  pb[iperm] = pb[i];
275  if (nonzero >= 0)
276  for (Int_t j = nonzero; j < i; j++)
277  r -= pLU[off_i+j]*pb[j];
278  else if (r != 0.0)
279  nonzero = i;
280  pb[i] = r;
281  }
282 
283  // Backward substitution
284  for (i = n-1; i >= 0; i--) {
285  const Int_t off_i = i*n;
286  Double_t r = pb[i];
287  for (Int_t j = i+1; j < n; j++)
288  r -= pLU[off_i+j]*pb[j];
289  pb[i] = r/pLU[off_i+i];
290  }
291 
292  return kTRUE;
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has *not*
297 /// been transformed. Solution returned in b.
298 
300 {
301  TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
302  R__ASSERT(b->IsValid());
303  if (TestBit(kSingular)) {
304  Error("Solve()","Matrix is singular");
305  return kFALSE;
306  }
307  if ( !TestBit(kDecomposed) ) {
308  if (!Decompose()) {
309  Error("Solve()","Decomposition failed");
310  return kFALSE;
311  }
312  }
313 
314  if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
315  Error("Solve(TMatrixDColumn &","vector and matrix incompatible");
316  return kFALSE;
317  }
318 
319  const Int_t n = fLU.GetNrows();
320  const Double_t *pLU = fLU.GetMatrixArray();
321 
322  Int_t i;
323 
324  // Check for zero diagonals
325  for (i = 0; i < n ; i++) {
326  const Int_t off_i = i*n;
327  if (TMath::Abs(pLU[off_i+i]) < fTol) {
328  Error("Solve(TMatrixDColumn &cb)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
329  return kFALSE;
330  }
331  }
332 
333  Double_t *pcb = cb.GetPtr();
334  const Int_t inc = cb.GetInc();
335 
336  // Transform b allowing for leading zeros
337  Int_t nonzero = -1;
338  for (i = 0; i < n; i++) {
339  const Int_t off_i = i*n;
340  const Int_t off_i2 = i*inc;
341  const Int_t iperm = fIndex[i];
342  const Int_t off_iperm = iperm*inc;
343  Double_t r = pcb[off_iperm];
344  pcb[off_iperm] = pcb[off_i2];
345  if (nonzero >=0)
346  for (Int_t j = nonzero; j <= i-1; j++)
347  r -= pLU[off_i+j]*pcb[j*inc];
348  else if (r != 0.0)
349  nonzero = i;
350  pcb[off_i2] = r;
351  }
352 
353  // Backward substitution
354  for (i = n-1; i >= 0; i--) {
355  const Int_t off_i = i*n;
356  const Int_t off_i2 = i*inc;
357  Double_t r = pcb[off_i2];
358  for (Int_t j = i+1; j < n; j++)
359  r -= pLU[off_i+j]*pcb[j*inc];
360  pcb[off_i2] = r/pLU[off_i+i];
361  }
362 
363  return kTRUE;
364 }
365 
366 
367 ////////////////////////////////////////////////////////////////////////////////
368 /// Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has *not*
369 /// been transformed. Solution returned in b.
370 
372 {
373  R__ASSERT(b.IsValid());
374  if (TestBit(kSingular)) {
375  Error("TransSolve()","Matrix is singular");
376  return kFALSE;
377  }
378  if ( !TestBit(kDecomposed) ) {
379  if (!Decompose()) {
380  Error("TransSolve()","Decomposition failed");
381  return kFALSE;
382  }
383  }
384 
385  if (fLU.GetNrows() != b.GetNrows() || fLU.GetRowLwb() != b.GetLwb()) {
386  Error("TransSolve(TVectorD &","vector and matrix incompatible");
387  return kFALSE;
388  }
389 
390  const Int_t n = fLU.GetNrows();
391 
392  const Double_t *pLU = fLU.GetMatrixArray();
393  Double_t *pb = b.GetMatrixArray();
394 
395  Int_t i;
396 
397  // Check for zero diagonals
398  for (i = 0; i < n ; i++) {
399  const Int_t off_i = i*n;
400  if (TMath::Abs(pLU[off_i+i]) < fTol) {
401  Error("TransSolve(TVectorD &b)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
402  return kFALSE;
403  }
404  }
405 
406  // Forward Substitution
407  for (i = 0; i < n; i++) {
408  const Int_t off_i = i*n;
409  Double_t r = pb[i];
410  for (Int_t j = 0; j < i ; j++) {
411  const Int_t off_j = j*n;
412  r -= pLU[off_j+i]*pb[j];
413  }
414  pb[i] = r/pLU[off_i+i];
415  }
416 
417  // Transform b allowing for leading zeros
418  Int_t nonzero = -1;
419  for (i = n-1 ; i >= 0; i--) {
420  Double_t r = pb[i];
421  if (nonzero >= 0) {
422  for (Int_t j = i+1; j <= nonzero; j++) {
423  const Int_t off_j = j*n;
424  r -= pLU[off_j+i]*pb[j];
425  }
426  } else if (r != 0.0)
427  nonzero = i;
428  const Int_t iperm = fIndex[i];
429  pb[i] = pb[iperm];
430  pb[iperm] = r;
431  }
432 
433  return kTRUE;
434 }
435 
436 ////////////////////////////////////////////////////////////////////////////////
437 /// Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has *not*
438 /// been transformed. Solution returned in b.
439 
441 {
442  TMatrixDBase *b = const_cast<TMatrixDBase *>(cb.GetMatrix());
443  R__ASSERT(b->IsValid());
444  if (TestBit(kSingular)) {
445  Error("TransSolve()","Matrix is singular");
446  return kFALSE;
447  }
448  if ( !TestBit(kDecomposed) ) {
449  if (!Decompose()) {
450  Error("TransSolve()","Decomposition failed");
451  return kFALSE;
452  }
453  }
454 
455  if (fLU.GetNrows() != b->GetNrows() || fLU.GetRowLwb() != b->GetRowLwb()) {
456  Error("TransSolve(TMatrixDColumn &","vector and matrix incompatible");
457  return kFALSE;
458  }
459 
460  const Int_t n = fLU.GetNrows();
461  const Int_t lwb = fLU.GetRowLwb();
462 
463  const Double_t *pLU = fLU.GetMatrixArray();
464 
465  Int_t i;
466 
467  // Check for zero diagonals
468  for (i = 0; i < n ; i++) {
469  const Int_t off_i = i*n;
470  if (TMath::Abs(pLU[off_i+i]) < fTol) {
471  Error("TransSolve(TMatrixDColumn &cb)","LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],fTol);
472  return kFALSE;
473  }
474  }
475 
476  // Forward Substitution
477  for (i = 0; i < n; i++) {
478  const Int_t off_i = i*n;
479  Double_t r = cb(i+lwb);
480  for (Int_t j = 0; j < i ; j++) {
481  const Int_t off_j = j*n;
482  r -= pLU[off_j+i]*cb(j+lwb);
483  }
484  cb(i+lwb) = r/pLU[off_i+i];
485  }
486 
487  // Transform b allowing for leading zeros
488  Int_t nonzero = -1;
489  for (i = n-1 ; i >= 0; i--) {
490  Double_t r = cb(i+lwb);
491  if (nonzero >= 0) {
492  for (Int_t j = i+1; j <= nonzero; j++) {
493  const Int_t off_j = j*n;
494  r -= pLU[off_j+i]*cb(j+lwb);
495  }
496  } else if (r != 0.0)
497  nonzero = i;
498  const Int_t iperm = fIndex[i];
499  cb(i+lwb) = cb(iperm+lwb);
500  cb(iperm+lwb) = r;
501  }
502 
503  return kTRUE;
504 }
505 
506 ////////////////////////////////////////////////////////////////////////////////
507 /// Calculate determinant det = d1*TMath::Power(2.,d2)
508 
510 {
511  if ( !TestBit(kDetermined) ) {
512  if ( !TestBit(kDecomposed) )
513  Decompose();
514  TDecompBase::Det(d1,d2);
515  fDet1 *= fSign;
517  }
518  d1 = fDet1;
519  d2 = fDet2;
520 }
521 
522 ////////////////////////////////////////////////////////////////////////////////
523 /// For a matrix A(m,m), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
524 /// (m x m) Ainv is returned .
525 
527 {
528  if (inv.GetNrows() != GetNrows() || inv.GetNcols() != GetNcols() ||
529  inv.GetRowLwb() != GetRowLwb() || inv.GetColLwb() != GetColLwb()) {
530  Error("Invert(TMatrixD &","Input matrix has wrong shape");
531  return kFALSE;
532  }
533 
534  inv.UnitMatrix();
535  const Bool_t status = MultiSolve(inv);
536 
537  return status;
538 }
539 
540 ////////////////////////////////////////////////////////////////////////////////
541 /// For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
542 /// (n x m) Ainv is returned .
543 
545 {
546  const Int_t rowLwb = GetRowLwb();
547  const Int_t rowUpb = rowLwb+GetNrows()-1;
548 
549  TMatrixD inv(rowLwb,rowUpb,rowLwb,rowUpb);
550  inv.UnitMatrix();
551  status = MultiSolve(inv);
552 
553  return inv;
554 }
555 
556 ////////////////////////////////////////////////////////////////////////////////
557 ///Print internals of this object
558 
559 void TDecompLU::Print(Option_t *opt) const
560 {
561  TDecompBase::Print(opt);
562  printf("fImplicitPivot = %d\n",fImplicitPivot);
563  printf("fSign = %f\n",fSign);
564  printf("fIndex:\n");
565  for (Int_t i = 0; i < fNIndex; i++)
566  printf("[%d] = %d\n",i,fIndex[i]);
567  fLU.Print("fLU");
568 }
569 
570 ////////////////////////////////////////////////////////////////////////////////
571 ///assignement operator
572 
574 {
575  if (this != &source) {
576  TDecompBase::operator=(source);
577  fLU.ResizeTo(source.fLU);
578  fLU = source.fLU;
579  fSign = source.fSign;
581  if (fNIndex != source.fNIndex) {
582  if (fIndex)
583  delete [] fIndex;
584  fNIndex = source.fNIndex;
585  fIndex = new Int_t[fNIndex];
586  }
587  if (fIndex) memcpy(fIndex,source.fIndex,fNIndex*sizeof(Int_t));
588  }
589  return *this;
590 }
591 
592 ////////////////////////////////////////////////////////////////////////////////
593 /// Crout/Doolittle algorithm of LU decomposing a square matrix, with implicit partial
594 /// pivoting. The decomposition is stored in fLU: U is explicit in the upper triag
595 /// and L is in multiplier form in the subdiagionals .
596 /// Row permutations are mapped out in fIndex. fSign, used for calculating the
597 /// determinant, is +/- 1 for even/odd row permutations. .
598 
600  Double_t tol,Int_t &nrZeros)
601 {
602  const Int_t n = lu.GetNcols();
603  Double_t *pLU = lu.GetMatrixArray();
604 
605  Double_t work[kWorkMax];
606  Bool_t isAllocated = kFALSE;
607  Double_t *scale = work;
608  if (n > kWorkMax) {
609  isAllocated = kTRUE;
610  scale = new Double_t[n];
611  }
612 
613  sign = 1.0;
614  nrZeros = 0;
615  // Find implicit scaling factors for each row
616  for (Int_t i = 0; i < n ; i++) {
617  const Int_t off_i = i*n;
618  Double_t max = 0.0;
619  for (Int_t j = 0; j < n; j++) {
620  const Double_t tmp = TMath::Abs(pLU[off_i+j]);
621  if (tmp > max)
622  max = tmp;
623  }
624  scale[i] = (max == 0.0 ? 0.0 : 1.0/max);
625  }
626 
627  for (Int_t j = 0; j < n; j++) {
628  const Int_t off_j = j*n;
629  // Run down jth column from top to diag, to form the elements of U.
630  for (Int_t i = 0; i < j; i++) {
631  const Int_t off_i = i*n;
632  Double_t r = pLU[off_i+j];
633  for (Int_t k = 0; k < i; k++) {
634  const Int_t off_k = k*n;
635  r -= pLU[off_i+k]*pLU[off_k+j];
636  }
637  pLU[off_i+j] = r;
638  }
639 
640  // Run down jth subdiag to form the residuals after the elimination of
641  // the first j-1 subdiags. These residuals divided by the appropriate
642  // diagonal term will become the multipliers in the elimination of the jth.
643  // subdiag. Find fIndex of largest scaled term in imax.
644 
645  Double_t max = 0.0;
646  Int_t imax = 0;
647  for (Int_t i = j; i < n; i++) {
648  const Int_t off_i = i*n;
649  Double_t r = pLU[off_i+j];
650  for (Int_t k = 0; k < j; k++) {
651  const Int_t off_k = k*n;
652  r -= pLU[off_i+k]*pLU[off_k+j];
653  }
654  pLU[off_i+j] = r;
655  const Double_t tmp = scale[i]*TMath::Abs(r);
656  if (tmp >= max) {
657  max = tmp;
658  imax = i;
659  }
660  }
661 
662  // Permute current row with imax
663  if (j != imax) {
664  const Int_t off_imax = imax*n;
665  for (Int_t k = 0; k < n; k++ ) {
666  const Double_t tmp = pLU[off_imax+k];
667  pLU[off_imax+k] = pLU[off_j+k];
668  pLU[off_j+k] = tmp;
669  }
670  sign = -sign;
671  scale[imax] = scale[j];
672  }
673  index[j] = imax;
674 
675  // If diag term is not zero divide subdiag to form multipliers.
676  if (pLU[off_j+j] != 0.0) {
677  if (TMath::Abs(pLU[off_j+j]) < tol)
678  nrZeros++;
679  if (j != n-1) {
680  const Double_t tmp = 1.0/pLU[off_j+j];
681  for (Int_t i = j+1; i < n; i++) {
682  const Int_t off_i = i*n;
683  pLU[off_i+j] *= tmp;
684  }
685  }
686  } else {
687  ::Error("TDecompLU::DecomposeLUCrout","matrix is singular");
688  if (isAllocated) delete [] scale;
689  return kFALSE;
690  }
691  }
692 
693  if (isAllocated)
694  delete [] scale;
695 
696  return kTRUE;
697 }
698 
699 ////////////////////////////////////////////////////////////////////////////////
700 /// LU decomposition using Gaussain Elimination with partial pivoting (See Golub &
701 /// Van Loan, Matrix Computations, Algorithm 3.4.1) of a square matrix .
702 /// The decomposition is stored in fLU: U is explicit in the upper triag and L is in
703 /// multiplier form in the subdiagionals . Row permutations are mapped out in fIndex.
704 /// fSign, used for calculating the determinant, is +/- 1 for even/odd row permutations.
705 /// Since this algorithm uses partial pivoting without scaling like in Crout/Doolitle.
706 /// it is somewhat faster but less precise .
707 
709  Double_t tol,Int_t &nrZeros)
710 {
711  const Int_t n = lu.GetNcols();
712  Double_t *pLU = lu.GetMatrixArray();
713 
714  sign = 1.0;
715  nrZeros = 0;
716 
717  index[n-1] = n-1;
718  for (Int_t j = 0; j < n-1; j++) {
719  const Int_t off_j = j*n;
720 
721  // Find maximum in the j-th column
722 
723  Double_t max = TMath::Abs(pLU[off_j+j]);
724  Int_t i_pivot = j;
725 
726  for (Int_t i = j+1; i < n; i++) {
727  const Int_t off_i = i*n;
728  const Double_t mLUij = TMath::Abs(pLU[off_i+j]);
729 
730  if (mLUij > max) {
731  max = mLUij;
732  i_pivot = i;
733  }
734  }
735 
736  if (i_pivot != j) {
737  const Int_t off_ipov = i_pivot*n;
738  for (Int_t k = 0; k < n; k++ ) {
739  const Double_t tmp = pLU[off_ipov+k];
740  pLU[off_ipov+k] = pLU[off_j+k];
741  pLU[off_j+k] = tmp;
742  }
743  sign = -sign;
744  }
745  index[j] = i_pivot;
746 
747  const Double_t mLUjj = pLU[off_j+j];
748 
749  if (mLUjj != 0.0) {
750  if (TMath::Abs(mLUjj) < tol)
751  nrZeros++;
752  for (Int_t i = j+1; i < n; i++) {
753  const Int_t off_i = i*n;
754  const Double_t mLUij = pLU[off_i+j]/mLUjj;
755  pLU[off_i+j] = mLUij;
756 
757  for (Int_t k = j+1; k < n; k++) {
758  const Double_t mLUik = pLU[off_i+k];
759  const Double_t mLUjk = pLU[off_j+k];
760  pLU[off_i+k] = mLUik-mLUij*mLUjk;
761  }
762  }
763  } else {
764  ::Error("TDecompLU::DecomposeLUGauss","matrix is singular");
765  return kFALSE;
766  }
767  }
768 
769  return kTRUE;
770 }
771 
772 ////////////////////////////////////////////////////////////////////////////////
773 /// Calculate matrix inversion through in place forward/backward substitution
774 
776 {
777  if (det)
778  *det = 0.0;
779 
780  if (lu.GetNrows() != lu.GetNcols() || lu.GetRowLwb() != lu.GetColLwb()) {
781  ::Error("TDecompLU::InvertLU","matrix should be square");
782  return kFALSE;
783  }
784 
785  const Int_t n = lu.GetNcols();
786  Double_t *pLU = lu.GetMatrixArray();
787 
788  Int_t worki[kWorkMax];
789  Bool_t isAllocatedI = kFALSE;
790  Int_t *index = worki;
791  if (n > kWorkMax) {
792  isAllocatedI = kTRUE;
793  index = new Int_t[n];
794  }
795 
796  Double_t sign = 1.0;
797  Int_t nrZeros = 0;
798  if (!DecomposeLUCrout(lu,index,sign,tol,nrZeros) || nrZeros > 0) {
799  if (isAllocatedI)
800  delete [] index;
801  ::Error("TDecompLU::InvertLU","matrix is singular, %d diag elements < tolerance of %.4e",nrZeros,tol);
802  return kFALSE;
803  }
804 
805  if (det) {
806  Double_t d1;
807  Double_t d2;
808  const TVectorD diagv = TMatrixDDiag_const(lu);
809  DiagProd(diagv,tol,d1,d2);
810  d1 *= sign;
811  *det = d1*TMath::Power(2.0,d2);
812  }
813 
814  // Form inv(U).
815 
816  Int_t j;
817 
818  for (j = 0; j < n; j++) {
819  const Int_t off_j = j*n;
820 
821  pLU[off_j+j] = 1./pLU[off_j+j];
822  const Double_t mLU_jj = -pLU[off_j+j];
823 
824 // Compute elements 0:j-1 of j-th column.
825 
826  Double_t *pX = pLU+j;
827  Int_t k;
828  for (k = 0; k <= j-1; k++) {
829  const Int_t off_k = k*n;
830  if ( pX[off_k] != 0.0 ) {
831  const Double_t tmp = pX[off_k];
832  for (Int_t i = 0; i <= k-1; i++) {
833  const Int_t off_i = i*n;
834  pX[off_i] += tmp*pLU[off_i+k];
835  }
836  pX[off_k] *= pLU[off_k+k];
837  }
838  }
839  for (k = 0; k <= j-1; k++) {
840  const Int_t off_k = k*n;
841  pX[off_k] *= mLU_jj;
842  }
843  }
844 
845  // Solve the equation inv(A)*L = inv(U) for inv(A).
846 
847  Double_t workd[kWorkMax];
848  Bool_t isAllocatedD = kFALSE;
849  Double_t *pWorkd = workd;
850  if (n > kWorkMax) {
851  isAllocatedD = kTRUE;
852  pWorkd = new Double_t[n];
853  }
854 
855  for (j = n-1; j >= 0; j--) {
856 
857  // Copy current column j of L to WORK and replace with zeros.
858  for (Int_t i = j+1; i < n; i++) {
859  const Int_t off_i = i*n;
860  pWorkd[i] = pLU[off_i+j];
861  pLU[off_i+j] = 0.0;
862  }
863 
864  // Compute current column of inv(A).
865 
866  if (j < n-1) {
867  const Double_t *mp = pLU+j+1; // Matrix row ptr
868  Double_t *tp = pLU+j; // Target vector ptr
869 
870  for (Int_t irow = 0; irow < n; irow++) {
871  Double_t sum = 0.;
872  const Double_t *sp = pWorkd+j+1; // Source vector ptr
873  for (Int_t icol = 0; icol < n-1-j ; icol++)
874  sum += *mp++ * *sp++;
875  *tp = -sum + *tp;
876  mp += j+1;
877  tp += n;
878  }
879  }
880  }
881 
882  if (isAllocatedD)
883  delete [] pWorkd;
884 
885  // Apply column interchanges.
886  for (j = n-1; j >= 0; j--) {
887  const Int_t jperm = index[j];
888  if (jperm != j) {
889  for (Int_t i = 0; i < n; i++) {
890  const Int_t off_i = i*n;
891  const Double_t tmp = pLU[off_i+jperm];
892  pLU[off_i+jperm] = pLU[off_i+j];
893  pLU[off_i+j] = tmp;
894  }
895  }
896  }
897 
898  if (isAllocatedI)
899  delete [] index;
900 
901  return kTRUE;
902 }
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:132
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual Int_t GetNrows() const
Definition: TDecompLU.h:52
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
Element * GetPtr() const
Double_t fSign
Definition: TDecompLU.h:33
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
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 fImplicitPivot
Definition: TDecompLU.h:29
virtual Bool_t Decompose()
Matrix A is decomposed in components U and L so that P * A = U * L If the decomposition succeeds...
Definition: TDecompLU.cxx:126
Int_t GetNrows() const
Definition: TVectorT.h:81
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
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=0)
Calculate matrix inversion through in place forward/backward substitution.
Definition: TDecompLU.cxx:775
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
int Int_t
Definition: RtypesCore.h:41
Int_t * fIndex
Definition: TDecompLU.h:32
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
Definition: TDecompLU.cxx:509
const Bool_t kFALSE
Definition: Rtypes.h:92
static Bool_t DecomposeLUCrout(TMatrixD &lu, Int_t *index, Double_t &sign, Double_t tol, Int_t &nrZeros)
Crout/Doolittle algorithm of LU decomposing a square matrix, with implicit partial pivoting...
Definition: TDecompLU.cxx:599
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:1201
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:732
ClassImp(TDecompLU) TDecompLU
Default constructor.
Definition: TDecompLU.cxx:15
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
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
Int_t GetColLwb() const
Definition: TDecompBase.h:77
Int_t fColLwb
Definition: TDecompBase.h:44
TDecompLU & operator=(const TDecompLU &source)
assignement operator
Definition: TDecompLU.cxx:573
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
Element * GetMatrixArray()
Definition: TVectorT.h:84
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:24
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Int_t fNIndex
Definition: TDecompLU.h:31
const double tol
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual Int_t GetNcols() const
Definition: TDecompLU.h:53
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:173
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:211
void ResetStatus()
Definition: TDecompBase.h:46
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetColLwb() const
Definition: TMatrixTBase.h:135
Double_t fDet2
Definition: TDecompBase.h:41
Double_t fDet1
Definition: TDecompBase.h:40
Int_t GetRowLwb() const
Definition: TDecompBase.h:76
TMatrixD fLU
Definition: TDecompLU.h:34
Int_t GetLwb() const
Definition: TVectorT.h:79
const TMatrixD GetMatrix()
Reconstruct the original matrix using the decomposition parts.
Definition: TDecompLU.cxx:151
double Double_t
Definition: RtypesCore.h:55
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
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)
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
Definition: TDecompLU.cxx:202
TMatrixD Invert()
Definition: TDecompLU.h:71
virtual Bool_t Solve(TVectorD &b)
Solve Ax=b assuming the LU form of A is stored in fLU, but assume b has not been transformed.
Definition: TDecompLU.cxx:233
Int_t GetInc() const
Int_t sign(Double_t x)
Definition: CsgOps.cxx:89
const Bool_t kTRUE
Definition: Rtypes.h:91
const Int_t n
Definition: legend1.C:16
virtual Bool_t TransSolve(TVectorD &b)
Solve A^T x=b assuming the LU form of A^T is stored in fLU, but assume b has not been transformed...
Definition: TDecompLU.cxx:371
void Print(Option_t *opt="") const
Print internals of this object.
Definition: TDecompLU.cxx:559
static Bool_t DecomposeLUGauss(TMatrixD &lu, Int_t *index, Double_t &sign, Double_t tol, Int_t &nrZeros)
LU decomposition using Gaussain Elimination with partial pivoting (See Golub & Van Loan...
Definition: TDecompLU.cxx:708