include/root/TDecompBase.h

0001 // @(#)root/matrix:$Id$
0002 // Authors: Fons Rademakers, Eddy Offermann   Dec 2003
0003
0004 /*************************************************************************
0005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOT_TDecompBase
0013 #define ROOT_TDecompBase
0014
0015 ///////////////////////////////////////////////////////////////////////////
0016 //                                                                       //
0017 // Decomposition Base class                                              //
0018 //                                                                       //
0019 // This class forms the base for all the decompositions methods in the   //
0020 // linear algebra package .                                              //
0021 //                                                                       //
0022 ///////////////////////////////////////////////////////////////////////////
0023
0024 #include "Rtypes.h"
0025
0026 #include "TMatrixD.h"
0027 #include "TMatrixDUtils.h"
0028 #include "TObject.h"
0029 #include "TVectorD.h"
0030
0031 #include <limits>
0032
0033 class TDecompBase : public TObject
0034 {
0035 protected :
0036    Double_t fTol;       // sqrt(epsilon); epsilon is smallest number number so that  1+epsilon > 1
0037    Double_t fDet1;      // determinant mantissa
0038    Double_t fDet2;      // determinant exponent for powers of 2
0039    Double_t fCondition; // matrix condition number
0040    Int_t    fRowLwb;    // Row    lower bound of decomposed matrix
0041    Int_t    fColLwb;    // Column lower bound of decomposed matrix
0042
0043    void          ResetStatus() { for (Int_t i = 14; i < 22; i++) ResetBit(BIT(i)); }
0044    Int_t         Hager      (Double_t& est,Int_t iter=5);
0045    static  void  DiagProd   (const TVectorD &diag,Double_t tol,Double_t &d1,Double_t &d2);
0046
0047    virtual const TMatrixDBase &GetDecompMatrix() const = 0;
0048
0049    enum EMatrixDecompStat {
0050 // clang++ <v20 (-Wshadow) complains about shadowing Getline.h global enum EGetLineMode. Let's silence warning:
0051 #if defined(__clang__) && __clang_major__ < 20
0052 #pragma clang diagnostic push
0053 #pragma clang diagnostic ignored "-Wshadow"
0054 #endif
0055       kInit       = BIT(14),
0056 #if defined(__clang__) && __clang_major__ < 20
0057 #pragma clang diagnostic pop
0058 #endif
0059       kPatternSet = BIT(15),
0060       kValuesSet  = BIT(16),
0061       kMatrixSet  = BIT(17),
0062       kDecomposed = BIT(18),
0063       kDetermined = BIT(19),
0064       kCondition  = BIT(20),
0065       kSingular   = BIT(21)
0066    };
0067
0068    enum {kWorkMax = 100}; // size of work array's in several routines
0069
0070 public :
0071    TDecompBase();
0072    TDecompBase(const TDecompBase &another);
0073    ~TDecompBase() override {};
0074
0075    inline  Double_t GetTol       () const { return fTol; }
0076    inline  Double_t GetDet1      () const { return fDet1; }
0077    inline  Double_t GetDet2      () const { return fDet2; }
0078    inline  Double_t GetCondition () const { return fCondition; }
0079    virtual Int_t    GetNrows     () const = 0;
0080    virtual Int_t    GetNcols     () const = 0;
0081    Int_t            GetRowLwb    () const { return fRowLwb; }
0082    Int_t            GetColLwb    () const { return fColLwb; }
0083    inline Double_t  SetTol       (Double_t tol);
0084
0085    virtual Double_t Condition  ();
0086    virtual void     Det        (Double_t &d1,Double_t &d2);
0087    virtual Bool_t   Decompose  ()                             = 0;
0088    virtual Bool_t   Solve      (      TVectorD &b)            = 0;
0089    virtual TVectorD Solve      (const TVectorD& b,Bool_t &ok) = 0;
0090    virtual Bool_t   Solve      (      TMatrixDColumn& b)      = 0;
0091    virtual Bool_t   TransSolve (      TVectorD &b)            = 0;
0092    virtual TVectorD TransSolve (const TVectorD &b,Bool_t &ok) = 0;
0093    virtual Bool_t   TransSolve (      TMatrixDColumn& b)      = 0;
0094
0095    virtual Bool_t   MultiSolve (TMatrixD &B);
0096
0097    void Print(Option_t *opt="") const override;
0098
0099    TDecompBase &operator= (const TDecompBase &source);
0100
0101    ClassDefOverride(TDecompBase,2) // Matrix Decomposition Base
0102 };
0103
0104 Double_t TDecompBase::SetTol(Double_t newTol)
0105 {
0106    const Double_t oldTol = fTol;
0107    if (newTol >= 0.0)
0108       fTol = newTol;
0109    return oldTol;
0110 }
0111
0112 Bool_t DefHouseHolder  (const TVectorD &vc,Int_t     lp,Int_t     l,Double_t &up,Double_t &b,Double_t tol=0.0);
0113 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TMatrixDRow &cr);
0114 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TMatrixDColumn &cc);
0115 void   ApplyHouseHolder(const TVectorD &vc,Double_t  up,Double_t  b,Int_t     lp,Int_t     l,TVectorD &cv);
0116 void   DefGivens       (      Double_t  v1,Double_t  v2,Double_t &c,Double_t &s);
0117 void   DefAplGivens    (      Double_t &v1,Double_t &v2,Double_t &c,Double_t &s);
0118 void   ApplyGivens     (      Double_t &z1,Double_t &z2,Double_t  c,Double_t  s);
0119
0120 #endif