library: libMatrix
#include "TDecompBase.h" |
TDecompBase
class description - source file - inheritance tree (.ps)
This is an abstract class, constructors will not be documented.
Look at the header to check for available constructors.
protected:
virtual const TMatrixDBase& GetDecompMatrix() const
Int_t Hager(Double_t& est, Int_t iter = 5)
void ResetStatus()
public:
virtual ~TDecompBase()
static TClass* Class()
virtual Double_t Condition()
virtual Bool_t Decompose()
virtual void Det(Double_t& d1, Double_t& d2)
static void DiagProd(const TVectorD& diag, Double_t tol, Double_t& d1, Double_t& d2)
Int_t GetColLwb() const
Double_t GetCondition() const
Double_t GetDet1() const
Double_t GetDet2() const
virtual Int_t GetNcols() const
virtual Int_t GetNrows() const
Int_t GetRowLwb() const
Double_t GetTol() const
virtual void Invert(TMatrixD& inv)
virtual TMatrixD Invert()
virtual TClass* IsA() const
virtual Bool_t MultiSolve(TMatrixD& B)
virtual Bool_t MultiSolve(TMatrixDSym& B)
TDecompBase& operator=(const TDecompBase& source)
virtual void Print(Option_t* opt) const
Double_t SetTol(Double_t tol)
virtual void ShowMembers(TMemberInspector& insp, char* parent)
virtual Bool_t Solve(TVectorD& b)
virtual TVectorD Solve(const TVectorD& b, Bool_t& ok)
virtual Bool_t Solve(TMatrixDColumn& b)
virtual void Streamer(TBuffer& b)
void StreamerNVirtual(TBuffer& b)
virtual Bool_t TransSolve(TVectorD& b)
virtual TVectorD TransSolve(const TVectorD& b, Bool_t& ok)
virtual Bool_t TransSolve(TMatrixDColumn& b)
protected:
Double_t fTol sqrt(epsilon); epsilon is smallest number number so that 1+epsilon > 1
Double_t fDet1 determinant mantissa
Double_t fDet2 determinant exponent for powers of 2
Double_t fCondition matrix condition number
Int_t fRowLwb Row lower bound of decomposed matrix
Int_t fColLwb Column lower bound of decomposed matrix
public:
static const TDecompBase::EMatrixDecompStat kInit
static const TDecompBase::EMatrixDecompStat kPatternSet
static const TDecompBase::EMatrixDecompStat kValuesSet
static const TDecompBase::EMatrixDecompStat kMatrixSet
static const TDecompBase::EMatrixDecompStat kDecomposed
static const TDecompBase::EMatrixDecompStat kDetermined
static const TDecompBase::EMatrixDecompStat kCondition
static const TDecompBase::EMatrixDecompStat kSingular
static const enum TDecompBase:: kWorkMax
Decomposition Base class
This class forms the base for all the decompositions methods in the
linear algebra package .
It or its derived classes have installed the methods to solve
equations,invert matrices and calculate determinants while monitoring
the accuracy.
When the constructor is called with a "const" matrix, the original
matrix "survives" (of course) . Some classes (LU) have also a non-
"const" constructor available . Here the original matrix is adopted
by the decomposition class (See Use in TMatrixDBase) to store the
decomposed matrix, thereby avoiding new memory allocation.
The decomposition (which is called by the constructor) fails when the
matrix is singular or not positive-definite in case of Cholesky
This can be checked before applying the decomposition by checking the
matrix through GetDecompMatrix()
Each derived class has always the following methods available:
Condition() :
In an iterative scheme the condition number for matrix inversion is
calculated . This number is of interest for estimating the accuracy
of x in the equation Ax=b
For example:
A is a (10x10) Hilbert matrix which looks deceivingly innocent
and simple, A(i,j) = 1/(i+j+1)
b(i) = Sum_j A(i,j), so a sum of a row in A
the solution is x(i) = 1. i=0,.,9
However,
TMatrixD m....; TVectorD b.....
TDecompLU lu(m); lu.SetTol(1.0e-12); lu.Solve(b); b.Print()
gives,
{1.000,1.000,1.000,1.000,0.998,1.000,0.993,1.001,0.996,1.000}
Looking at the condition number, this is in line with expected the
accuracy . The condition number is 3.957e+12 . As a simple rule of
thumb, a condition number of 1.0e+n means that you lose up to n
digits of accuracy in a solution . Since doubles are stored with 15
digits, we can expect the accuracy to be as small as 3 digits .
Det(Double_t &d1,Double_t &d2)
The determinant is d1*TMath::Power(2.,d2)
Expressing the determinant this way makes under/over-flow very
unlikely .
Decompose()
Here the actually decomposition is performed . One can change the
matrix A after the decomposition constructor has been called
without effecting the decomposition result
Solve(TVectorD &b)
Solve A x = b . x is supplied through the argument and replaced with
the solution .
TransSolve(TVectorD &b)
Solve A^T x = b . x is supplied through the argument and replaced
with the solution .
MultiSolve(TMatrixD &B)
MultiSolve(TMatrixDSym &B)
Solve A X = B . where X and are now matrices . X is supplied through
the argument and replaced with the solution .
Invert(TMatrixD &inv)
This is of course just a call to MultiSolve with as input argument
the unit matrix . Note that for a matrix a(m,n) with m > n a
pseudo-inverse is calculated .
Tolerances and Scaling
----------------------
The tolerance parameter (which is a member of this base class) plays
a crucial role in all operations of the decomposition classes . It
gives the user a powerful tool to monitor and steer the operations
Its default value is sqrt(epsilon) where 1+epsilon = 1
If you do not want to be bothered by the following considerations,
like in most other linear algebra packages, just set the tolerance
with SetTol to an arbitrary small number .
The tolerance number is used by each decomposition method to decide
whether the matrix is near singular, except of course SVD which can
handle singular matrices .
For each decomposition this will be checked in a different way; in LU
the matrix is considered singular when, at some point in the
decomposition, a diagonal element < fTol . Therefore, we had to set in
the example above of the (10x10) Hilbert, which is near singular, the
tolerance on 10e-12 . (The fact that we have to set the tolerance <
sqrt(epsilon) is a clear indication that we are losing precision .)
If the matrix is flagged as being singular, operations with the
decomposition will fail and will return matrices/vectors that are
invalid .
The observant reader will notice that by scaling the complete matrix
by some small number the decomposition will detect a singular matrix .
In this case the user will have to reduce the tolerance number by this
factor . (For CPU time saving we decided not to make this an automatic
procedure) .
Code for this could look as follows:
const TMatrixD ab = Abs(a);
const Int_t imax = TMath::LocMax(ab.GetNoElements(),
ab.GetMatrixArray());
const Double_t max_abs = ab.GetMatrixArray()[imax];
const Double_t scale = TMath::Min(max_abs,1.);
a.SetTol(a.GetTol()*scale);
For usage examples see $ROOTSYS/test/stressLinear.cxx
Int_t Hager(Double_t &est,Int_t iter)
void DiagProd(const TVectorD &diag,Double_t tol,Double_t &d1,Double_t &d2)
Double_t Condition()
Bool_t MultiSolve(TMatrixD &B)
Solve set of equations with RHS in columns of B
Bool_t MultiSolve(TMatrixDSym &B)
Solve set of equations with RHS in columns of B
void Invert(TMatrixD &inv)
For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
The user should always supply a matrix of size (m x m) !
If m > n , only the (n x m) part of the returned (pseudo inverse) matrix
should be used .
TMatrixD Invert()
For a matrix A(m,n), its inverse A_inv is defined as A * A_inv = A_inv * A = unit
(n x m) Ainv is returned .
void Det(Double_t &d1,Double_t &d2)
void Print(Option_t * /*opt*/) const
Inline Functions
void ~TDecompBase()
void ResetStatus()
const TMatrixDBase& GetDecompMatrix() const
Double_t GetTol() const
Double_t GetDet1() const
Double_t GetDet2() const
Double_t GetCondition() const
Int_t GetNrows() const
Int_t GetNcols() const
Int_t GetRowLwb() const
Int_t GetColLwb() const
Double_t SetTol(Double_t tol)
Bool_t Decompose()
Bool_t Solve(TVectorD& b)
TVectorD Solve(const TVectorD& b, Bool_t& ok)
Bool_t Solve(TMatrixDColumn& b)
Bool_t TransSolve(TVectorD& b)
TVectorD TransSolve(const TVectorD& b, Bool_t& ok)
Bool_t TransSolve(TMatrixDColumn& b)
TDecompBase& operator=(const TDecompBase& source)
TClass* Class()
TClass* IsA() const
void ShowMembers(TMemberInspector& insp, char* parent)
void Streamer(TBuffer& b)
void StreamerNVirtual(TBuffer& b)
Last update: root/matrix:$Name: $:$Id: TDecompBase.cxx,v 1.14 2004/06/13 14:53:15 brun Exp $
Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.