Logo ROOT   6.21/01
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 #include "TDecompBase.h"
22 #include "TMatrixDSparse.h"
23 #include "TArrayD.h"
24 #include "TArrayI.h"
25 
26 // globals
27 const Double_t kInitTreatAsZero = 1.0e-12;
29 const Double_t kInitPrecision = 1.0e-7;
30 
31 // the Threshold Pivoting parameter may need to be increased during
32 // the algorithm if poor precision is obtained from the linear
33 // solves. kThresholdPivoting indicates the largest value we are
34 // willing to tolerate.
35 
37 
38 // the factor in the range (1,inf) by which kThresholdPivoting is
39 // increased when it is found to be inadequate.
40 
42 
43 class TDecompSparse : public TDecompBase
44 {
45 protected :
46 
48 
49  Int_t fIcntl[31]; // integer control numbers
50  Double_t fCntl[6]; // float control numbers
51  Int_t fInfo[21]; // array used for communication between programs
52 
53  Double_t fPrecision; // precision we demand from the linear system solver. If it isn't
54  // attained on the first solve, we use iterative refinement and
55  // possibly refactorization with a higher value of kThresholdPivoting.
56 
57  TArrayI fIkeep; // pivot sequence and temporary storage information
63  TArrayD fW; // temporary storage for the factorization
64 
65  Double_t fIPessimism; // amounts by which to increase allocated factorization space when
66  Double_t fRPessimism; // inadequate space is detected. fIPessimism is for array "fIw",
67  // fRPessimism is for the array "fact".
68 
69  TMatrixDSparse fA; // original matrix; needed for the iterative solving procedure
72  TArrayD fFact; // size of fFact array; may be increased during the numerical factorization
73  // if the estimate obtained during the symbolic factorization proves to be inadequate.
76 
78  static void CopyUpperTriang (const TMatrixDSparse &a,Double_t *b);
79 
80  void InitParam();
81  static void InitPivot (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,
82  TArrayI &Aiw,TArrayI &Aikeep,TArrayI &Aiw1,Int_t &nsteps,
83  const Int_t iflag,Int_t *icntl,Double_t *cntl,Int_t *info,Double_t &ops);
84  static void Factor (const Int_t n,const Int_t nz,TArrayI &Airn,TArrayI &Aicn,TArrayD &Aa,
85  TArrayI &Aiw,TArrayI &Aikeep,const Int_t nsteps,Int_t &maxfrt,
86  TArrayI &Aiw1,Int_t *icntl,Double_t *cntl,Int_t *info);
87  static void Solve (const Int_t n,TArrayD &Aa,TArrayI &Aiw,TArrayD &Aw,const Int_t maxfrt,
88  TVectorD &b,TArrayI &Aiw1,const Int_t nsteps,Int_t *icntl,Int_t *info);
89 
90  static void InitPivot_sub1 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *iw,Int_t *ipe,
91  Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
92  static void InitPivot_sub2 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *nv,
93  Int_t *nxt,Int_t *lst,Int_t *ipd,Int_t *flag,const Int_t iovflo,Int_t &ncmpa,
94  const Double_t fratio);
95  static void InitPivot_sub2a(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t &ncmpa);
96  static void InitPivot_sub3 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *iw,
97  Int_t *ipe,Int_t *iq,Int_t *flag,Int_t &iwfr,Int_t *icntl,Int_t *info);
98  static void InitPivot_sub4 (const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t *ips,
99  Int_t *ipv,Int_t *nv,Int_t *flag,Int_t &ncmpa);
100  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,
101  Int_t &nsteps,const Int_t nemin);
102  static void InitPivot_sub6 (const Int_t n,const Int_t nz,Int_t *irn,Int_t *icn,Int_t *perm,Int_t *na,
103  Int_t *ne,Int_t *nd,const Int_t nsteps,Int_t *lstki,Int_t *lstkr,Int_t *iw,
104  Int_t *info,Double_t &ops);
105 
106  static void Factor_sub1 (const Int_t n,const Int_t nz,Int_t &nz1,Double_t *a,const Int_t la,
107  Int_t *irn,Int_t *icn,Int_t *iw,const Int_t liw,Int_t *perm,Int_t *iw2,
108  Int_t *icntl,Int_t *info);
109  static void Factor_sub2 (const Int_t n,const Int_t nz,Double_t *a,const Int_t la,Int_t *iw,
110  const Int_t liw,Int_t *perm,Int_t *nstk,const Int_t nsteps,Int_t &maxfrt,
111  Int_t *nelim,Int_t *iw2,Int_t *icntl,Double_t *cntl,Int_t *info);
112  static void Factor_sub3 (Double_t *a,Int_t *iw,Int_t &j1,Int_t &j2,const Int_t itop,const Int_t ireal,
113  Int_t &ncmpbr,Int_t &ncmpbi);
114 
115  static void Solve_sub1 (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
116  const Int_t nblk,Int_t &latop,Int_t *icntl);
117  static void Solve_sub2 (const Int_t n,Double_t *a,Int_t *iw,Double_t *w,Double_t *rhs,Int_t *iw2,
118  const Int_t nblk,const Int_t latop,Int_t *icntl);
119  static Int_t IDiag (Int_t ix,Int_t iy) { return ((iy-1)*(2*ix-iy+2))/2; }
120 
121  inline Int_t IError () { return fInfo[2]; }
122  inline Int_t MinRealWorkspace() { return fInfo[5]; }
123  inline Int_t MinIntWorkspace () { return fInfo[6]; }
124  inline Int_t ErrorFlag () { return fInfo[1]; }
125 
126  // Takes values in the range [0,1]. Larger values enforce greater stability in
127  // the factorization as they insist on larger pivots. Smaller values preserve
128  // sparsity at the cost of using smaller pivots.
129 
130  inline Double_t GetThresholdPivoting() { return fCntl[1]; }
131  inline Double_t GetTreatAsZero () { return fCntl[3]; }
132 
133  // The factorization will not accept a pivot whose absolute value is less than fCntl[3] as
134  // a 1x1 pivot or as the off-diagonal in a 2x2 pivot.
135 
136  inline void SetThresholdPivoting(Double_t piv) { fCntl[1] = piv; }
137  inline void SetTreatAsZero (Double_t tol) { fCntl[3] = tol; }
138 
139  virtual const TMatrixDBase &GetDecompMatrix() const { MayNotUse("GetDecompMatrix()"); return fA; }
140 
141 public :
142 
143  TDecompSparse();
144  TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose);
145  TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose);
147  TDecompSparse(const TDecompSparse &another);
148  virtual ~TDecompSparse() {}
149 
150  inline void SetVerbose (Int_t v) { fVerbose = (v) ? 1 : 0;
151  if (fVerbose) { fIcntl[1] = fIcntl[2] = 1; fIcntl[3] = 2; }
152  else { fIcntl[1] = fIcntl[2] = fIcntl[3] = 0; }
153  }
154  virtual Int_t GetNrows () const { return fA.GetNrows(); }
155  virtual Int_t GetNcols () const { return fA.GetNcols(); }
156 
157  virtual void SetMatrix (const TMatrixDSparse &a);
158 
159  virtual Bool_t Decompose ();
160  virtual Bool_t Solve ( TVectorD &b);
161  virtual TVectorD Solve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
162  virtual Bool_t Solve ( TMatrixDColumn & /*b*/)
163  { MayNotUse("Solve(TMatrixDColumn &)"); return kFALSE; }
164  virtual Bool_t TransSolve ( TVectorD &b) { return Solve(b); }
165  virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
166  virtual Bool_t TransSolve ( TMatrixDColumn & /*b*/)
167  { MayNotUse("TransSolve(TMatrixDColumn &)"); return kFALSE; }
168 
169  virtual void Det (Double_t &/*d1*/,Double_t &/*d2*/)
170  { MayNotUse("Det(Double_t&,Double_t&)"); }
171 
172  void Print(Option_t *opt ="") const; // *MENU*
173 
174  TDecompSparse &operator= (const TDecompSparse &source);
175 
176  ClassDef(TDecompSparse,1) // Matrix Decompositition LU
177 };
178 
179 #endif
void SetTreatAsZero(Double_t tol)
virtual Int_t GetNcols() const
Double_t fCntl[6]
Definition: TDecompSparse.h:50
TMatrixDSparse fA
Definition: TDecompSparse.h:69
TArrayI fColFact
Definition: TDecompSparse.h:75
Double_t fRPessimism
Definition: TDecompSparse.h:66
virtual Int_t GetNrows() 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:65
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:41
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
TArrayI fIkeep
Definition: TDecompSparse.h:57
virtual Bool_t TransSolve(TMatrixDColumn &)
virtual const TMatrixDBase & GetDecompMatrix() const
void Print(Option_t *opt="") const
Print class members.
TDecompSparse()
Default constructor.
Decomposition Base class.
Definition: TDecompBase.h:33
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
Array of integers (32 bits per element).
Definition: TArrayI.h:27
virtual Bool_t Decompose()
Decomposition engine .
Int_t MinIntWorkspace()
Double_t x[n]
Definition: legend1.C:17
#define ClassDef(name, id)
Definition: Rtypes.h:326
static Int_t IDiag(Int_t ix, Int_t iy)
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 decomposition step.
void InitParam()
initializing control parameters
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:933
TArrayI fRowFact
Definition: TDecompSparse.h:74
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix ...
Int_t fInfo[21]
Definition: TDecompSparse.h:51
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.
auto * a
Definition: textangle.C:12
Double_t GetThresholdPivoting()
Sparse Symmetric Decomposition class.
Definition: TDecompSparse.h:43
void SetThresholdPivoting(Double_t piv)
Linear Algebra Package.
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:49
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t fPrecision
Definition: TDecompSparse.h:53
virtual TVectorD TransSolve(const TVectorD &b, Bool_t &ok)
const Double_t kInitPrecision
Definition: TDecompSparse.h:29
double Double_t
Definition: RtypesCore.h:55
virtual ~TDecompSparse()
const Double_t kThresholdPivotingMax
Definition: TDecompSparse.h:36
Int_t MinRealWorkspace()
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
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:27
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()
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
const Double_t kInitThresholdPivoting
Definition: TDecompSparse.h:28
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.