71 const Int_t nrows = row_upb-row_lwb+1;
79 fLU.
ResizeTo(row_lwb,row_lwb+nrows-1,row_lwb,row_lwb+nrows-1);
90 Error(
"TDecompLU(const TMatrixD &",
"matrix should be square");
131 Error(
"Decompose()",
"Matrix has not been set");
154 Error(
"GetMatrix()",
"Matrix is singular");
159 Error(
"GetMatrix()",
"Decomposition failed");
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.;
183 for (
Int_t i = n-1; i >= 0; i--) {
188 for (
Int_t k = 0; k <
n; k++) {
190 pA[off_j+k] = pA[off_i+k];
208 Error(
"TDecompLU(const TMatrixD &",
"matrix should be square");
237 Error(
"Solve()",
"Matrix is singular");
242 Error(
"Solve()",
"Decomposition failed");
248 Error(
"Solve(TVectorD &",
"vector and matrix incompatible");
260 for (i = 0; i <
n ; i++) {
263 Error(
"Solve(TVectorD &b)",
"LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],
fTol);
270 for (i = 0; i <
n; i++) {
276 for (
Int_t j = nonzero; j < i; j++)
277 r -= pLU[off_i+j]*pb[j];
284 for (i = n-1; i >= 0; 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];
304 Error(
"Solve()",
"Matrix is singular");
309 Error(
"Solve()",
"Decomposition failed");
315 Error(
"Solve(TMatrixDColumn &",
"vector and matrix incompatible");
325 for (i = 0; i <
n ; i++) {
328 Error(
"Solve(TMatrixDColumn &cb)",
"LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],
fTol);
338 for (i = 0; i <
n; i++) {
340 const Int_t off_i2 = i*inc;
342 const Int_t off_iperm = iperm*inc;
344 pcb[off_iperm] = pcb[off_i2];
346 for (
Int_t j = nonzero; j <= i-1; j++)
347 r -= pLU[off_i+j]*pcb[j*inc];
354 for (i = n-1; i >= 0; i--) {
356 const Int_t off_i2 = i*inc;
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];
375 Error(
"TransSolve()",
"Matrix is singular");
380 Error(
"TransSolve()",
"Decomposition failed");
386 Error(
"TransSolve(TVectorD &",
"vector and matrix incompatible");
398 for (i = 0; i <
n ; i++) {
401 Error(
"TransSolve(TVectorD &b)",
"LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],
fTol);
407 for (i = 0; i <
n; i++) {
410 for (
Int_t j = 0; j < i ; j++) {
412 r -= pLU[off_j+i]*pb[j];
414 pb[i] = r/pLU[off_i+i];
419 for (i = n-1 ; i >= 0; i--) {
422 for (
Int_t j = i+1; j <= nonzero; j++) {
424 r -= pLU[off_j+i]*pb[j];
445 Error(
"TransSolve()",
"Matrix is singular");
450 Error(
"TransSolve()",
"Decomposition failed");
456 Error(
"TransSolve(TMatrixDColumn &",
"vector and matrix incompatible");
468 for (i = 0; i <
n ; i++) {
471 Error(
"TransSolve(TMatrixDColumn &cb)",
"LU[%d,%d]=%.4e < %.4e",i,i,pLU[off_i+i],
fTol);
477 for (i = 0; i <
n; i++) {
480 for (
Int_t j = 0; j < i ; j++) {
482 r -= pLU[off_j+i]*cb(j+lwb);
484 cb(i+lwb) = r/pLU[off_i+i];
489 for (i = n-1 ; i >= 0; i--) {
492 for (
Int_t j = i+1; j <= nonzero; j++) {
494 r -= pLU[off_j+i]*cb(j+lwb);
499 cb(i+lwb) = cb(iperm+lwb);
530 Error(
"Invert(TMatrixD &",
"Input matrix has wrong shape");
563 printf(
"fSign = %f\n",
fSign);
566 printf(
"[%d] = %d\n",i,
fIndex[i]);
575 if (
this != &source) {
616 for (
Int_t i = 0; i <
n ; i++) {
619 for (
Int_t j = 0; j <
n; j++) {
624 scale[i] = (max == 0.0 ? 0.0 : 1.0/max);
627 for (
Int_t j = 0; j <
n; j++) {
630 for (
Int_t i = 0; i < j; i++) {
633 for (
Int_t k = 0; k < i; k++) {
635 r -= pLU[off_i+k]*pLU[off_k+j];
647 for (
Int_t i = j; i <
n; i++) {
650 for (
Int_t k = 0; k < j; k++) {
652 r -= pLU[off_i+k]*pLU[off_k+j];
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];
671 scale[imax] = scale[j];
676 if (pLU[off_j+j] != 0.0) {
680 const Double_t tmp = 1.0/pLU[off_j+j];
681 for (
Int_t i = j+1; i <
n; i++) {
687 ::Error(
"TDecompLU::DecomposeLUCrout",
"matrix is singular");
688 if (isAllocated)
delete [] scale;
718 for (
Int_t j = 0; j < n-1; j++) {
726 for (
Int_t i = j+1; i <
n; i++) {
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];
747 const Double_t mLUjj = pLU[off_j+j];
752 for (
Int_t i = j+1; i <
n; i++) {
754 const Double_t mLUij = pLU[off_i+j]/mLUjj;
755 pLU[off_i+j] = mLUij;
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;
764 ::Error(
"TDecompLU::DecomposeLUGauss",
"matrix is singular");
781 ::Error(
"TDecompLU::InvertLU",
"matrix should be square");
790 Int_t *index = worki;
792 isAllocatedI =
kTRUE;
801 ::Error(
"TDecompLU::InvertLU",
"matrix is singular, %d diag elements < tolerance of %.4e",nrZeros,tol);
818 for (j = 0; j <
n; j++) {
821 pLU[off_j+j] = 1./pLU[off_j+j];
822 const Double_t mLU_jj = -pLU[off_j+j];
828 for (k = 0; k <= j-1; k++) {
830 if ( pX[off_k] != 0.0 ) {
832 for (
Int_t i = 0; i <= k-1; i++) {
834 pX[off_i] += tmp*pLU[off_i+k];
836 pX[off_k] *= pLU[off_k+k];
839 for (k = 0; k <= j-1; k++) {
851 isAllocatedD =
kTRUE;
855 for (j = n-1; j >= 0; j--) {
858 for (
Int_t i = j+1; i <
n; i++) {
860 pWorkd[i] = pLU[off_i+j];
870 for (
Int_t irow = 0; irow <
n; irow++) {
873 for (
Int_t icol = 0; icol < n-1-j ; icol++)
874 sum += *mp++ * *sp++;
886 for (j = n-1; j >= 0; j--) {
887 const Int_t jperm = index[j];
889 for (
Int_t i = 0; i <
n; i++) {
891 const Double_t tmp = pLU[off_i+jperm];
892 pLU[off_i+jperm] = pLU[off_i+j];
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
static long int sum(long int i)
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
void Print(Option_t *opt="") const
Print internals of this object.
virtual Bool_t Decompose()
Matrix A is decomposed in components U and L so that P * A = U * L If the decomposition succeeds...
TDecompLU()
Default constructor.
virtual const Element * GetMatrixArray() const
Decomposition Base class.
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
static Bool_t InvertLU(TMatrixD &a, Double_t tol, Double_t *det=0)
Calculate matrix inversion through in place forward/backward substitution.
virtual void Det(Double_t &d1, Double_t &d2)
Calculate determinant det = d1*TMath::Power(2.,d2)
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...
static constexpr double mL
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...
virtual Int_t GetNcols() const
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
TDecompLU & operator=(const TDecompLU &source)
assignment operator
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
Element * GetMatrixArray()
TMatrixT< Double_t > TMatrixD
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
void Print(Option_t *opt="") const
Print class members.
const TMatrixD GetMatrix()
Reconstruct the original matrix using the decomposition parts.
virtual Int_t GetNrows() const
virtual void Det(Double_t &d1, Double_t &d2)
Matrix determinant det = d1*TMath::Power(2.,d2)
void Print(Option_t *name="") const
Print the matrix as a table of elements.
const TMatrixTBase< Element > * GetMatrix() const
virtual void SetMatrix(const TMatrixD &a)
Set matrix to be decomposed.
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.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
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...
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...