library: libMatrix #include "TDecompBase.h"

# TDecompBase

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

## class TDecompBase : public TObject

| Inheritance Chart:

`TObject`

<-

`TDecompBase`

<-

`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:

          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)

### 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.

     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. *