ROOT_528-00b_version: include/TDecompBase.h Source File

00001 // @(#)root/matrix:$Id: TDecompBase.h 20882 2007-11-19 11:31:26Z rdm $
00002 // Authors: Fons Rademakers, Eddy Offermann   Dec 2003
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #ifndef ROOT_TDecompBase
00013 #define ROOT_TDecompBase
00014 
00015 ///////////////////////////////////////////////////////////////////////////
00016 //                                                                       //
00017 // Decomposition Base class                                              //
00018 //                                                                       //
00019 // This class forms the base for all the decompositions methods in the   //
00020 // linear algebra package .                                              //
00021 //                                                                       //
00022 ///////////////////////////////////////////////////////////////////////////
00023 
00024 #include <limits>
00025 
00026 #ifndef ROOT_TMatrixD
00027 #include "TMatrixD.h"
00028 #endif
00029 #ifndef ROOT_TMatrixDUtils
00030 #include "TMatrixDUtils.h"
00031 #endif
00032 #ifndef ROOT_TVectorD
00033 #include "TVectorD.h"
00034 #endif
00035 
00036 class TDecompBase : public TObject
00037 {
00038 protected :
00039    Double_t fTol;       // sqrt(epsilon); epsilon is smallest number number so that  1+epsilon > 1
00040    Double_t fDet1;      // determinant mantissa
00041    Double_t fDet2;      // determinant exponent for powers of 2
00042    Double_t fCondition; // matrix condition number
00043    Int_t    fRowLwb;    // Row    lower bound of decomposed matrix
00044    Int_t    fColLwb;    // Column lower bound of decomposed matrix
00045 
00046    void          ResetStatus() { for (Int_t i = 14; i < 22; i++) ResetBit(BIT(i)); }
00047    Int_t         Hager      (Double_t& est,Int_t iter=5);
00048    static  void  DiagProd   (const TVectorD &diag,Double_t tol,Double_t &d1,Double_t &d2);
00049 
00050    virtual const TMatrixDBase &GetDecompMatrix() const = 0;
00051 
00052    enum EMatrixDecompStat {
00053       kInit       = BIT(14),
00054       kPatternSet = BIT(15),
00055       kValuesSet  = BIT(16),
00056       kMatrixSet  = BIT(17),
00057       kDecomposed = BIT(18),
00058       kDetermined = BIT(19),
00059       kCondition  = BIT(20),
00060       kSingular   = BIT(21)
00061    };
00062 
00063    enum {kWorkMax = 100}; // size of work array's in several routines
00064 
00065 public :
00066    TDecompBase();
00067    TDecompBase(const TDecompBase &another);
00068    virtual ~TDecompBase() {};
00069 
00070    inline  Double_t GetTol       () const { return fTol; }
00071    inline  Double_t GetDet1      () const { return fDet1; }
00072    inline  Double_t GetDet2      () const { return fDet2; }
00073    inline  Double_t GetCondition () const { return fCondition; }
00074    virtual Int_t    GetNrows     () const = 0;
00075    virtual Int_t    GetNcols     () const = 0;
00076    Int_t            GetRowLwb    () const { return fRowLwb; }
00077    Int_t            GetColLwb    () const { return fColLwb; }
00078    inline Double_t  SetTol       (Double_t tol);
00079 
00080    virtual Double_t Condition  ();
00081    virtual void     Det        (Double_t &d1,Double_t &d2);
00082    virtual Bool_t   Decompose  ()                             = 0;
00083    virtual Bool_t   Solve      (      TVectorD &b)            = 0;
00084    virtual TVectorD Solve      (const TVectorD& b,Bool_t &ok) = 0;
00085    virtual Bool_t   Solve      (      TMatrixDColumn& b)      = 0;
00086    virtual Bool_t   TransSolve (      TVectorD &b)            = 0;
00087    virtual TVectorD TransSolve (const TVectorD &b,Bool_t &ok) = 0;
00088    virtual Bool_t   TransSolve (      TMatrixDColumn& b)      = 0;
00089 
00090    virtual Bool_t   MultiSolve (TMatrixD &B);
00091 
00092    void Print(Option_t *opt="") const;
00093 
00094    TDecompBase &operator= (const TDecompBase &source);
00095 
00096    ClassDef(TDecompBase,2) // Matrix Decomposition Base
00097 };
00098 
00099 Double_t TDecompBase::SetTol(Double_t newTol) 
00100 {
00101    const Double_t oldTol = fTol; 
00102    if (newTol >= 0.0) 
00103       fTol = newTol; 
00104    return oldTol; 
00105 }
00106 
00107 Bool_t DefHouseHolder  (const TVectorD &vc,Int_t     lp,Int_t     l,Double_t &up,Double_t &b,Double_t tol=0.0);
00108 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TMatrixDRow &cr);
00109 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TMatrixDColumn &cc);
00110 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TVectorD &cv);
00111 void   DefGivens       (      Double_t  v1,Double_t  v2,Double_t &c,Double_t &s);
00112 void   DefAplGivens    (      Double_t &v1,Double_t &v2,Double_t &c,Double_t &s);
00113 void   ApplyGivens     (      Double_t &z1,Double_t &z2,Double_t  c,Double_t  s);
00114 
00115 #endif