library: libMatrix
#include "TDecompBase.h"

TDecompBase


class description - source file - inheritance tree (.pdf)

class TDecompBase : public TObject

Inheritance Chart:
TObject
<-
TDecompBase
<-
TDecompBK
TDecompChol
TDecompLU
TDecompQRH
TDecompSparse
TDecompSVD
 
    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)

Data Members


    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

Class Description

                                                                       
 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.