ROOT  6.06/09
Reference Guide
TDecompSparse.h
Go to the documentation of this file.
1 // @(#)root/matrix:$Id$
2 // Authors: Fons Rademakers, Eddy Offermann Apr 2004
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #ifndef ROOT_TDecompSparse
13 #define ROOT_TDecompSparse
14 
15 ///////////////////////////////////////////////////////////////////////////
16 // //
17 // Sparse Decomposition class //
18 // //
19 ///////////////////////////////////////////////////////////////////////////
20 
21 #ifndef ROOT_TDecompBase
22 #include "TDecompBase.h"
23 #endif
24 #ifndef ROOT_TMatrixDSparse
25 #include "TMatrixDSparse.h"
26 #endif
27 #ifndef ROOT_TArrayD
28 #include "TArrayD.h"
29 #endif
30 #ifndef ROOT_TArrayI
31 #include "TArrayI.h"
32 #endif
33 
34 // globals
35 const Double_t kInitTreatAsZero = 1.0e-12;
37 const Double_t kInitPrecision = 1.0e-7;
38 
39 // the Threshold Pivoting parameter may need to be increased during
40 // the algorithm if poor precision is obtained from the linear
41 // solves. kThresholdPivoting indicates the largest value we are
42 // willing to tolerate.
43 
45 
46 // the factor in the range (1,inf) by which kThresholdPivoting is
47 // increased when it is found to be inadequate.
48 
50 
51 class TDecompSparse : public TDecompBase
52 {
53 protected :
54 
56 
57  Int_t fIcntl[31]; // integer control numbers
58  Double_t fCntl[6]; // float control numbers
59  Int_t fInfo[21]; // array used for communication between programs
60 
61  Double_t fPrecision; // precision we demand from the linear system solver. If it isn't
62  // attained on the first solve, we use iterative refinement and
63  // possibly refactorization with a higher value of kThresholdPivoting.
64 
65  TArrayI fIkeep; // pivot sequence and temporary storage information
71  TArrayD fW; // temporary storage for the factorization
72 
73  Double_t fIPessimism; // amounts by which to increase allocated factorization space when
74  Double_t fRPessimism; // inadequate space is detected. fIPessimism is for array "fIw",
75  // fRPessimism is for the array "fact".
76 
77  TMatrixDSparse fA; // original matrix; needed for the iterative solving procedure
80  TArrayD fFact; // size of fFact array; may be increased during the numerical factorization
81  // if the estimate obtained during the symbolic factorization proves to be inadequate.
84 
86  static void CopyUpperTriang (const TMatrixDSparse &a,Double_t *b);
87 
88  void InitParam();
89  static void InitPivot (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,
90  TArrayI &Aiw,TArrayI &Aikeep,TArrayI &Aiw1,Int_t &nsteps,
91  const Int_t iflag,Int_t *icntl,Double_t *cntl,Int_t *info,Double_t &ops);
92  static void Factor (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,TArrayD &Aa,
93  TArrayI &Aiw,TArrayI &Aikeep,const Int_t nsteps,Int_t &maxfrt,
94  TArrayI &Aiw1,Int_t *icntl,Double_t *cntl,Int_t *info);
95  static void Solve (const Int_t n,TArrayD &Aa,TArrayI &Aiw,TArrayD &Aw,const Int_t maxfrt,
96  TVectorD &b,TArrayI &Aiw1,const Int_t nsteps,Int_t *icntl,Int_t *info);
97 
98  static void InitPivot_sub1 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *iw,Int_t *ipe,
99  Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
100  static void InitPivot_sub2 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *nv,
101  Int_t *nxt,Int_t *lst,Int_t *ipd,Int_t *flag,const Int_t iovflo,Int_t &ncmpa,
102  const Double_t fratio);
103  static void InitPivot_sub2a(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t &ncmpa);
104  static void InitPivot_sub3 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *iw,
105  Int_t *ipe,Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
106  static void InitPivot_sub4 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *ips,
107  Int_t *ipv,Int_t *nv,Int_t *flag,Int_t &ncmpa);
108  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,
109  Int_t &nsteps,const Int_t nemin);
110  static void InitPivot_sub6 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *na,
111  Int_t *ne,Int_t *nd,const Int_t nsteps,Int_t *lstki,Int_t *lstkr,Int_t *iw,
112  Int_t *info,Double_t &ops);
113 
114  static void Factor_sub1 (const Int_t n,const Int_t nz,Int_t &nz1,Double_t *a,const Int_t la,
115  Int_t *irn,Int_t *icn,Int_t *iw,const Int_t liw,Int_t *perm,Int_t *iw2,
116  Int_t *icntl,Int_t *info);
117  static void Factor_sub2 (const Int_t n,const Int_t nz,Double_t *a,const Int_t la,Int_t *iw,
118  const Int_t liw,Int_t *perm,Int_t *nstk,const Int_t nsteps,Int_t &maxfrt,
119  Int_t *nelim,Int_t *iw2,Int_t *icntl,Double_t *cntl,Int_t *info);
120  static void Factor_sub3 (Double_t *a,Int_t *iw,Int_t &j1,Int_t &j2,const Int_t itop,const Int_t ireal,
121  Int_t &ncmpbr,Int_t &ncmpbi);
122 
123  static void Solve_sub1 (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
124  const Int_t nblk,Int_t &latop,Int_t *icntl);
125  static void Solve_sub2 (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
126  const Int_t nblk,const Int_t latop,Int_t *icntl);
127  static Int_t IDiag (Int_t ix,Int_t iy) { return ((iy-1)*(2*ix-iy+2))/2; }
128 
129  inline Int_t IError () { return fInfo[2]; }
130  inline Int_t MinRealWorkspace() { return fInfo[5]; }
131  inline Int_t MinIntWorkspace () { return fInfo[6]; }
132  inline Int_t ErrorFlag () { return fInfo[1]; }
133 
134  // Takes values in the range [0,1]. Larger values enforce greater stability in
135  // the factorization as they insist on larger pivots. Smaller values preserve
136  // sparsity at the cost of using smaller pivots.
137 
138  inline Double_t GetThresholdPivoting() { return fCntl[1]; }
139  inline Double_t GetTreatAsZero () { return fCntl[3]; }
140 
141  // The factorization will not accept a pivot whose absolute value is less than fCntl[3] as
142  // a 1x1 pivot or as the off-diagonal in a 2x2 pivot.
143 
144  inline void SetThresholdPivoting(Double_t piv) { fCntl[1] = piv; }
145  inline void SetTreatAsZero (Double_t tol) { fCntl[3] = tol; }
146 
147  virtual const TMatrixDBase &GetDecompMatrix() const { MayNotUse("GetDecompMatrix()"); return fA; }
148 
149 public :
150 
151  TDecompSparse();
152  TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose);
153  TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose);
155  TDecompSparse(const TDecompSparse &another);
156  virtual ~TDecompSparse() {}
157 
158  inline void SetVerbose (Int_t v) { fVerbose = (v) ? 1 : 0;
159  if (fVerbose) { fIcntl[1] = fIcntl[2] = 1; fIcntl[3] = 2; }
160  else { fIcntl[1] = fIcntl[2] = fIcntl[3] = 0; }
161  }
162  virtual Int_t GetNrows () const { return fA.GetNrows(); }
163  virtual Int_t GetNcols () const { return fA.GetNcols(); }
164 
165  virtual void SetMatrix (const TMatrixDSparse &a);
166 
167  virtual Bool_t Decompose ();
168  virtual Bool_t Solve ( TVectorD &b);
169  virtual TVectorD Solve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
170  virtual Bool_t Solve ( TMatrixDColumn & /*b*/)
171  { MayNotUse("Solve(TMatrixDColumn &)"); return kFALSE; }
172  virtual Bool_t TransSolve ( TVectorD &b) { return Solve(b); }
173  virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
174  virtual Bool_t TransSolve ( TMatrixDColumn & /*b*/)
175  { MayNotUse("TransSolve(TMatrixDColumn &)"); return kFALSE; }
176 
177  virtual void Det (Double_t &/*d1*/,Double_t &/*d2*/)
178  { MayNotUse("Det(Double_t&,Double_t&)"); }
179 
180  void Print(Option_t *opt ="") const; // *MENU*
181 
182  TDecompSparse &operator= (const TDecompSparse &source);
183 
184  ClassDef(TDecompSparse,1) // Matrix Decompositition LU
185 };
186 
187 #endif
void SetTreatAsZero(Double_t tol)
Double_t fCntl[6]
Definition: TDecompSparse.h:58
TMatrixDSparse fA
Definition: TDecompSparse.h:77
TArrayI fColFact
Definition: TDecompSparse.h:83
Double_t fRPessimism
Definition: TDecompSparse.h:74
virtual const TMatrixDBase & GetDecompMatrix() const
virtual Bool_t TransSolve(TVectorD &b)
virtual void SetMatrix(const TMatrixDSparse &a)
Set matrix to be decomposed .
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, Int_t &nsteps, const Int_t nemin)
Help routine for pivoting setup.
Double_t fIPessimism
Definition: TDecompSparse.h:73
const char Option_t
Definition: RtypesCore.h:62
static void CopyUpperTriang(const TMatrixDSparse &a, Double_t *b)
Static function, copying the non-zero entries in the upper triangle to array b .
static void InitPivot_sub1(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *iw, Int_t *ipe, Int_t *iq, Int_t *flag, Int_t &iwfr, Int_t *icntl, Int_t *info)
Help routine for pivoting setup.
const Double_t kThresholdPivotingFactor
Definition: TDecompSparse.h:49
TArrayI fIkeep
Definition: TDecompSparse.h:65
void MayNotUse(const char *method) const
Use this method to signal that a method (defined in a base class) may not be called in a derived clas...
Definition: TObject.cxx:971
virtual Bool_t TransSolve(TMatrixDColumn &)
static void Factor_sub3(Double_t *a, Int_t *iw, Int_t &j1, Int_t &j2, const Int_t itop, const Int_t ireal, Int_t &ncmpbr, Int_t &ncmpbi)
Help routine for factorization.
TDecompSparse & operator=(const TDecompSparse &source)
Assignment operator.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
Array of integers (32 bits per element).
Definition: TArrayI.h:29
virtual Bool_t Decompose()
Decomposition engine .
Int_t MinIntWorkspace()
Double_t x[n]
Definition: legend1.C:17
#define ClassDef(name, id)
Definition: Rtypes.h:254
static Int_t IDiag(Int_t ix, Int_t iy)
virtual Int_t GetNcols() const
void Print(Option_t *opt="") const
Print class members.
static void Solve(const Int_t n, TArrayD &Aa, TArrayI &Aiw, TArrayD &Aw, const Int_t maxfrt, TVectorD &b, TArrayI &Aiw1, const Int_t nsteps, Int_t *icntl, Int_t *info)
Main routine for solving Ax=b.
static void Solve_sub1(const Int_t n, Double_t *a, Int_t *iw, Double_t *w, Double_t *rhs, Int_t *iw2, const Int_t nblk, Int_t &latop, Int_t *icntl)
Help routine for solving.
virtual void Det(Double_t &, Double_t &)
Matrix determinant det = d1*TMath::Power(2.,d2)
static void Factor(const Int_t n, const Int_t nz, TArrayI &Airn, TArrayI &Aicn, TArrayD &Aa, TArrayI &Aiw, TArrayI &Aikeep, const Int_t nsteps, Int_t &maxfrt, TArrayI &Aiw1, Int_t *icntl, Double_t *cntl, Int_t *info)
Factorization routine, the workhorse for the decompostion step.
void InitParam()
initializing control parameters
TArrayI fRowFact
Definition: TDecompSparse.h:82
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix ...
const double tol
Int_t fInfo[21]
Definition: TDecompSparse.h:59
static void Solve_sub2(const Int_t n, Double_t *a, Int_t *iw, Double_t *w, Double_t *rhs, Int_t *iw2, const Int_t nblk, const Int_t latop, Int_t *icntl)
Help routine for solving.
SVector< double, 2 > v
Definition: Dict.h:5
bool verbose
Double_t GetThresholdPivoting()
void SetThresholdPivoting(Double_t piv)
static void Factor_sub2(const Int_t n, const Int_t nz, Double_t *a, const Int_t la, Int_t *iw, const Int_t liw, Int_t *perm, Int_t *nstk, const Int_t nsteps, Int_t &maxfrt, Int_t *nelim, Int_t *iw2, Int_t *icntl, Double_t *cntl, Int_t *info)
Help routine for factorization.
Int_t fIcntl[31]
Definition: TDecompSparse.h:57
Double_t fPrecision
Definition: TDecompSparse.h:61
virtual TVectorD TransSolve(const TVectorD &b, Bool_t &ok)
const Double_t kInitPrecision
Definition: TDecompSparse.h:37
virtual Int_t GetNrows() const
double Double_t
Definition: RtypesCore.h:55
virtual ~TDecompSparse()
const Double_t kThresholdPivotingMax
Definition: TDecompSparse.h:44
Int_t MinRealWorkspace()
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
Array of doubles (64 bits per element).
Definition: TArrayD.h:29
Double_t GetTreatAsZero()
virtual Bool_t Solve(TMatrixDColumn &)
static void InitPivot_sub2(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t *nv, Int_t *nxt, Int_t *lst, Int_t *ipd, Int_t *flag, const Int_t iovflo, Int_t &ncmpa, const Double_t fratio)
Help routine for pivoting setup.
static void InitPivot_sub4(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t *ips, Int_t *ipv, Int_t *nv, Int_t *flag, Int_t &ncmpa)
Help routine for pivoting setup.
const Double_t kInitTreatAsZero
Definition: TDecompSparse.h:35
static void InitPivot_sub2a(const Int_t n, Int_t *ipe, Int_t *iw, const Int_t lw, Int_t &iwfr, Int_t &ncmpa)
Help routine for pivoting setup.
Int_t ErrorFlag()
const Double_t kInitThresholdPivoting
Definition: TDecompSparse.h:36
int * iq
Definition: THbookFile.cxx:86
void SetVerbose(Int_t v)
virtual TVectorD Solve(const TVectorD &b, Bool_t &ok)
static void InitPivot(const Int_t n, const Int_t nz, TArrayI &Airn, TArrayI &Aicn, TArrayI &Aiw, TArrayI &Aikeep, TArrayI &Aiw1, Int_t &nsteps, const Int_t iflag, Int_t *icntl, Double_t *cntl, Int_t *info, Double_t &ops)
Setup Pivoting variables.
static void InitPivot_sub6(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *perm, Int_t *na, Int_t *ne, Int_t *nd, const Int_t nsteps, Int_t *lstki, Int_t *lstkr, Int_t *iw, Int_t *info, Double_t &ops)
Help routine for pivoting setup.
const Int_t n
Definition: legend1.C:16
static void Factor_sub1(const Int_t n, const Int_t nz, Int_t &nz1, Double_t *a, const Int_t la, Int_t *irn, Int_t *icn, Int_t *iw, const Int_t liw, Int_t *perm, Int_t *iw2, Int_t *icntl, Int_t *info)
Help routine for factorization.
static void InitPivot_sub3(const Int_t n, const Int_t nz, Int_t *irn, Int_t *icn, Int_t *perm, Int_t *iw, Int_t *ipe, Int_t *iq, Int_t *flag, Int_t &iwfr, Int_t *icntl, Int_t *info)
Help routine for pivoting setup.