Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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;
28 const Double_t kInitThresholdPivoting = 1.0e-8;
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 
36 const Double_t kThresholdPivotingMax = 1.0e-2;
37 
38 // the factor in the range (1,inf) by which kThresholdPivoting is
39 // increased when it is found to be inadequate.
40 
41 const Double_t kThresholdPivotingFactor = 10.0;
42 
43 class TDecompSparse : public TDecompBase
44 {
45 protected :
46 
47  Int_t fVerbose;
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
58  TArrayI fIw;
59  TArrayI fIw1;
60  TArrayI fIw2;
61  Int_t fNsteps;
62  Int_t fMaxfrt;
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
70  Int_t fNrows;
71  Int_t fNnonZeros;
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.
74  TArrayI fRowFact;
75  TArrayI fColFact;
76 
77  static Int_t NonZerosUpperTriang(const TMatrixDSparse &a);
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);
146  TDecompSparse(const TMatrixDSparse &a,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