Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/root/TDecompSparse.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // @(#)root/matrix:$Id$
0002 // Authors: Fons Rademakers, Eddy Offermann   Apr 2004
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_TDecompSparse
0013 #define ROOT_TDecompSparse
0014 
0015 ///////////////////////////////////////////////////////////////////////////
0016 //                                                                       //
0017 // Sparse Decomposition class                                            //
0018 //                                                                       //
0019 ///////////////////////////////////////////////////////////////////////////
0020 
0021 #include "TDecompBase.h"
0022 #include "TMatrixDSparse.h"
0023 #include "TArrayD.h"
0024 #include "TArrayI.h"
0025 
0026 // globals
0027 const Double_t kInitTreatAsZero         = 1.0e-12;
0028 const Double_t kInitThresholdPivoting   = 1.0e-8;
0029 const Double_t kInitPrecision           = 1.0e-7;
0030 
0031 // the Threshold Pivoting parameter may need to be increased during
0032 // the algorithm if poor precision is obtained from the linear
0033 // solves.  kThresholdPivoting indicates the largest value we are
0034 // willing to tolerate.
0035 
0036 const Double_t kThresholdPivotingMax    = 1.0e-2;
0037 
0038 // the factor in the range (1,inf) by which kThresholdPivoting is
0039 // increased when it is found to be inadequate.
0040 
0041 const Double_t kThresholdPivotingFactor = 10.0;
0042 
0043 class TDecompSparse : public TDecompBase
0044 {
0045 protected :
0046 
0047    Int_t     fVerbose;
0048 
0049    Int_t     fIcntl[31]; // integer control numbers
0050    Double_t  fCntl[6];   // float control numbers
0051    Int_t     fInfo[21];  // array used for communication between programs
0052 
0053    Double_t  fPrecision; // precision we demand from the linear system solver. If it isn't
0054                          // attained on the first solve, we use iterative refinement and
0055                          // possibly refactorization with a higher value of kThresholdPivoting.
0056 
0057    TArrayI   fIkeep;     // pivot sequence and temporary storage information
0058    TArrayI   fIw;
0059    TArrayI   fIw1;
0060    TArrayI   fIw2;
0061    Int_t     fNsteps;
0062    Int_t     fMaxfrt;
0063    TArrayD   fW;         // temporary storage for the factorization
0064 
0065    Double_t  fIPessimism; // amounts by which to increase allocated factorization space when
0066    Double_t  fRPessimism; // inadequate space is detected. fIPessimism is for array "fIw",
0067                           // fRPessimism is for the array "fact".
0068 
0069    TMatrixDSparse fA; // original matrix; needed for the iterative solving procedure
0070    Int_t          fNrows;
0071    Int_t          fNnonZeros;
0072    TArrayD        fFact; // size of fFact array; may be increased during the numerical factorization
0073                          // if the estimate obtained during the symbolic factorization proves to be inadequate.
0074    TArrayI        fRowFact;
0075    TArrayI        fColFact;
0076 
0077    static Int_t NonZerosUpperTriang(const TMatrixDSparse &a);
0078    static void  CopyUpperTriang    (const TMatrixDSparse &a,Double_t *b);
0079 
0080           void  InitParam();
0081    static void  InitPivot (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,
0082                            TArrayI &Aiw,TArrayI &Aikeep,TArrayI &Aiw1,Int_t &nsteps,
0083                            const Int_t iflag,Int_t *icntl,Double_t *cntl,Int_t *info,Double_t &ops);
0084    static void   Factor   (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,TArrayD &Aa,
0085                            TArrayI &Aiw,TArrayI &Aikeep,const Int_t nsteps,Int_t &maxfrt,
0086                            TArrayI &Aiw1,Int_t *icntl,Double_t *cntl,Int_t *info);
0087    static void   Solve    (const Int_t n,TArrayD &Aa,TArrayI &Aiw,TArrayD &Aw,const Int_t maxfrt,
0088                            TVectorD &b,TArrayI &Aiw1,const Int_t nsteps,Int_t *icntl,Int_t *info);
0089 
0090    static void   InitPivot_sub1 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *iw,Int_t *ipe,
0091                                  Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
0092    static void   InitPivot_sub2 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *nv,
0093                                  Int_t *nxt,Int_t *lst,Int_t *ipd,Int_t *flag,const Int_t iovflo,Int_t &ncmpa,
0094                                  const Double_t fratio);
0095    static void   InitPivot_sub2a(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t &ncmpa);
0096    static void   InitPivot_sub3 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *iw,
0097                                  Int_t *ipe,Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
0098    static void   InitPivot_sub4 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *ips,
0099                                  Int_t *ipv,Int_t *nv,Int_t *flag,Int_t &ncmpa);
0100    static void   InitPivot_sub5 (const Int_t n,Int_t *ipe,Int_t *nv,Int_t *ips,Int_t *ne,Int_t *na,Int_t *nd,
0101                                  Int_t &nsteps,const Int_t nemin);
0102    static void   InitPivot_sub6 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *na,
0103                                  Int_t *ne,Int_t *nd,const Int_t nsteps,Int_t *lstki,Int_t *lstkr,Int_t *iw,
0104                                  Int_t *info,Double_t &ops);
0105 
0106    static void   Factor_sub1    (const Int_t n,const Int_t nz,Int_t &nz1,Double_t *a,const Int_t la,
0107                                  Int_t *irn,Int_t *icn,Int_t *iw,const Int_t liw,Int_t *perm,Int_t *iw2,
0108                                  Int_t *icntl,Int_t *info);
0109    static void   Factor_sub2    (const Int_t n,const Int_t nz,Double_t *a,const Int_t la,Int_t *iw,
0110                                  const Int_t liw,Int_t *perm,Int_t *nstk,const Int_t nsteps,Int_t &maxfrt,
0111                                  Int_t *nelim,Int_t *iw2,Int_t *icntl,Double_t *cntl,Int_t *info);
0112    static void   Factor_sub3    (Double_t *a,Int_t *iw,Int_t &j1,Int_t &j2,const Int_t itop,const Int_t ireal,
0113                                  Int_t &ncmpbr,Int_t &ncmpbi);
0114 
0115    static void   Solve_sub1     (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
0116                                  const Int_t nblk,Int_t &latop,Int_t *icntl);
0117    static void   Solve_sub2     (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
0118                                  const Int_t nblk,const Int_t latop,Int_t *icntl);
0119    static Int_t  IDiag          (Int_t ix,Int_t iy) { return ((iy-1)*(2*ix-iy+2))/2; }
0120 
0121    inline Int_t IError          () { return fInfo[2]; }
0122    inline Int_t MinRealWorkspace() { return fInfo[5]; }
0123    inline Int_t MinIntWorkspace () { return fInfo[6]; }
0124    inline Int_t ErrorFlag       () { return fInfo[1]; }
0125 
0126    // Takes values in the range [0,1]. Larger values enforce greater stability in
0127    // the factorization as they insist on larger pivots. Smaller values preserve
0128    // sparsity at the cost of using smaller pivots.
0129 
0130    inline Double_t GetThresholdPivoting() { return fCntl[1]; }
0131    inline Double_t GetTreatAsZero      () { return fCntl[3]; }
0132 
0133    // The factorization will not accept a pivot whose absolute value is less than fCntl[3] as
0134    // a 1x1 pivot or as the off-diagonal in a 2x2 pivot.
0135 
0136    inline void     SetThresholdPivoting(Double_t piv) { fCntl[1] = piv; }
0137    inline void     SetTreatAsZero      (Double_t tol) { fCntl[3] = tol; }
0138 
0139    const TMatrixDBase &GetDecompMatrix() const override { MayNotUse("GetDecompMatrix()"); return fA; }
0140 
0141 public :
0142 
0143    TDecompSparse();
0144    TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose);
0145    TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose);
0146    TDecompSparse(const TMatrixDSparse &a,Int_t verbose);
0147    TDecompSparse(const TDecompSparse &another);
0148    ~TDecompSparse() override {}
0149 
0150    inline  void     SetVerbose (Int_t v) { fVerbose = (v) ? 1 : 0;
0151                                             if (fVerbose) { fIcntl[1] = fIcntl[2] = 1; fIcntl[3] = 2; }
0152                                             else          { fIcntl[1] = fIcntl[2] = fIcntl[3] = 0; }
0153                                          }
0154    Int_t    GetNrows   () const override { return fA.GetNrows(); }
0155    Int_t    GetNcols   () const override { return fA.GetNcols(); }
0156 
0157    virtual void     SetMatrix  (const TMatrixDSparse &a);
0158 
0159    Bool_t   Decompose  () override;
0160    Bool_t   Solve      (      TVectorD &b) override;
0161    TVectorD Solve      (const TVectorD& b,Bool_t &ok) override { TVectorD x = b; ok = Solve(x); return x; }
0162    Bool_t   Solve      (      TMatrixDColumn & /*b*/) override
0163                                { MayNotUse("Solve(TMatrixDColumn &)"); return kFALSE; }
0164    Bool_t   TransSolve (      TVectorD &b) override            { return Solve(b); }
0165    TVectorD TransSolve (const TVectorD& b,Bool_t &ok) override { TVectorD x = b; ok = Solve(x); return x; }
0166    Bool_t   TransSolve (      TMatrixDColumn & /*b*/) override
0167                                { MayNotUse("TransSolve(TMatrixDColumn &)"); return kFALSE; }
0168 
0169    void     Det        (Double_t &/*d1*/,Double_t &/*d2*/) override
0170                                 { MayNotUse("Det(Double_t&,Double_t&)"); }
0171 
0172    void Print(Option_t *opt ="") const override; // *MENU*
0173 
0174    TDecompSparse &operator= (const TDecompSparse &source);
0175 
0176    ClassDefOverride(TDecompSparse,1) // Matrix Decompositition LU
0177 };
0178 
0179 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!