library: libMatrix
#include "TDecompBase.h" |
TDecompBase
class description - source file - inheritance tree (.pdf)
This is an abstract class, constructors will not be documented.
Look at the header to check for available constructors.
protected:
static void DiagProd(const TVectorD& diag, Double_t tol, Double_t& d1, Double_t& d2)
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)
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 TClass* IsA() const
virtual Bool_t MultiSolve(TMatrixD& 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.
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)
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 Double_t max_abs = Abs(a).Max();
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
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.16 2004/10/16 18:09:16 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.