ROOT   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/*************************************************************************
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/** \class TDecompLU
18 \ingroup Matrix
19
20 LU Decomposition class
21
22 Decompose a general n x n matrix A into P A = L U
23
24 where P is a permutation matrix, L is unit lower triangular and U
25 is upper triangular.
26 L is stored in the strict lower triangular part of the matrix fLU.
27 The diagonal elements of L are unity and are not stored.
28 U is stored in the diagonal and upper triangular part of the matrix
29 fU.
30 P is stored in the index array fIndex : j = fIndex[i] indicates that
31 row j and row i should be swapped .
32
33 fSign gives the sign of the permutation, (-1)^n, where n is the
34 number of interchanges in the permutation.
35
36 fLU has the same indexing range as matrix A .
37
38 The decomposition fails if a diagonal element of abs(fLU) is == 0,
39 The matrix fUL is made invalid .
40*/
41
42////////////////////////////////////////////////////////////////////////////////
43/// Default constructor
44
46{
47 fSign = 0.0;
48 fNIndex = 0;
49 fIndex = 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));
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;
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
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
560{
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///assignment operator
572
574{
575 if (this != &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 Gaussian 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}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
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.
void Print(Option_t *opt="") const
Print class members.
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
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
LU Decomposition class.
Definition: TDecompLU.h:24
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
Definition: TDecompLU.cxx:509
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
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
Definition: TDecompLU.cxx:202
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,...
Definition: TDecompLU.cxx:708
virtual Int_t GetNcols() const
Definition: TDecompLU.h:51
TDecompLU()
Default constructor.
Definition: TDecompLU.cxx:45
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
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
virtual Int_t GetNrows() const
Definition: TDecompLU.h:50
Double_t fSign
Definition: TDecompLU.h:31
const TMatrixD GetMatrix()
Reconstruct the original matrix using the decomposition parts.
Definition: TDecompLU.cxx:151
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.
Definition: TDecompLU.cxx:599
TMatrixD Invert()
Definition: TDecompLU.h:69
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
TMatrixD fLU
Definition: TDecompLU.h:32
void Print(Option_t *opt="") const
Print internals of this object.
Definition: TDecompLU.cxx:559
Int_t * fIndex
Definition: TDecompLU.h:30
Int_t fImplicitPivot
Definition: TDecompLU.h:27
TDecompLU & operator=(const TDecompLU &source)
assignment operator
Definition: TDecompLU.cxx:573
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:122
Int_t GetColLwb() const
Definition: TMatrixTBase.h:125
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
const TMatrixTBase< Element > * GetMatrix() const
Int_t GetInc() const
Element * GetPtr() const
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:1213
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
const Int_t n
Definition: legend1.C:16
static constexpr double mL
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
auto * a
Definition: textangle.C:12
static long int sum(long int i)
Definition: Factory.cxx:2275