89 if (
a.GetNrows() !=
a.GetNcols() ||
a.GetRowLwb() !=
a.GetColLwb()) {
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");
183 for (
Int_t i =
n-1; i >= 0; i--) {
188 for (
Int_t k = 0; k <
n; k++) {
207 if (
a.GetNrows() !=
a.GetNcols() ||
a.GetRowLwb() !=
a.GetColLwb()) {
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++) {
270 for (i = 0; i <
n; i++) {
284 for (i =
n-1; i >= 0; 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++) {
354 for (i =
n-1; i >= 0; 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++) {
407 for (i = 0; i <
n; i++) {
419 for (i =
n-1 ; i >= 0; i--) {
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++) {
489 for (i =
n-1 ; i >= 0; i--) {
530 Error(
"Invert(TMatrixD &",
"Input matrix has wrong shape");
616 for (
Int_t i = 0; i <
n ; i++) {
624 scale[i] = (max == 0.0 ? 0.0 : 1.0/max);
630 for (
Int_t i = 0; i <
j; i++) {
633 for (
Int_t k = 0; k < i; k++) {
647 for (
Int_t i =
j; i <
n; i++) {
650 for (
Int_t k = 0; k <
j; k++) {
665 for (
Int_t k = 0; k <
n; k++ ) {
681 for (
Int_t i =
j+1; i <
n; i++) {
687 ::Error(
"TDecompLU::DecomposeLUCrout",
"matrix is singular");
726 for (
Int_t i =
j+1; i <
n; i++) {
738 for (
Int_t k = 0; k <
n; k++ ) {
752 for (
Int_t i =
j+1; i <
n; i++) {
757 for (
Int_t k =
j+1; k <
n; k++) {
764 ::Error(
"TDecompLU::DecomposeLUGauss",
"matrix is singular");
780 if (
lu.GetNrows() !=
lu.GetNcols() ||
lu.GetRowLwb() !=
lu.GetColLwb()) {
781 ::Error(
"TDecompLU::InvertLU",
"matrix should be square");
801 ::Error(
"TDecompLU::InvertLU",
"matrix is singular, %d diag elements < tolerance of %.4e",
nrZeros,
tol);
818 for (
j = 0;
j <
n;
j++) {
828 for (k = 0; k <=
j-1; k++) {
832 for (
Int_t i = 0; i <= k-1; i++) {
839 for (k = 0; k <=
j-1; k++) {
855 for (
j =
n-1;
j >= 0;
j--) {
858 for (
Int_t i =
j+1; i <
n; i++) {
886 for (
j =
n-1;
j >= 0;
j--) {
889 for (
Int_t i = 0; i <
n; i++) {
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.
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
Decomposition Base class.
static void DiagProd(const TVectorD &diag, Double_t tol, Double_t &d1, Double_t &d2)
virtual Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B.
TDecompBase & operator=(const TDecompBase &source)
Assignment operator.
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)
Int_t GetNrows() const override
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
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.
const TMatrixD GetMatrix()
Reconstruct the original matrix using the decomposition parts.
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.
void Det(Double_t &d1, Double_t &d2) override
Calculate determinant det = d1*TMath::Power(2.,d2)
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.
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.
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
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
const TMatrixTBase< Element > * GetMatrix() const
const Element * GetMatrixArray() const override
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
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
static uint64_t sum(uint64_t i)