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