Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
27const Double_t kInitTreatAsZero = 1.0e-12;
29const 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
44{
45protected :
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
141public :
142
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; }
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
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:101
bool Bool_t
Definition RtypesCore.h:63
double Double_t
Definition RtypesCore.h:59
const char Option_t
Definition RtypesCore.h:66
#define ClassDef(name, id)
Definition Rtypes.h:325
const Double_t kThresholdPivotingFactor
const Double_t kThresholdPivotingMax
const Double_t kInitThresholdPivoting
const Double_t kInitTreatAsZero
const Double_t kInitPrecision
int * iq
Array of doubles (64 bits per element).
Definition TArrayD.h:27
Array of integers (32 bits per element).
Definition TArrayI.h:27
Decomposition Base class.
Definition TDecompBase.h:34
Sparse Symmetric Decomposition class.
virtual ~TDecompSparse()
Double_t fRPessimism
virtual Bool_t Decompose()
Decomposition engine .
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.
TDecompSparse()
Default constructor.
Double_t GetTreatAsZero()
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.
static Int_t IDiag(Int_t ix, Int_t iy)
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_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.
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a)
Static function, returning the number of non-zero entries in the upper triangular matrix .
virtual void Det(Double_t &, Double_t &)
Matrix determinant det = d1*TMath::Power(2.,d2)
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 Int_t GetNrows() const
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.
static void CopyUpperTriang(const TMatrixDSparse &a, Double_t *b)
Static function, copying the non-zero entries in the upper triangle to array b .
Int_t MinRealWorkspace()
Double_t GetThresholdPivoting()
void SetTreatAsZero(Double_t tol)
void Print(Option_t *opt="") const
Print class members.
TDecompSparse & operator=(const TDecompSparse &source)
Assignment operator.
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.
Int_t MinIntWorkspace()
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.
virtual void SetMatrix(const TMatrixDSparse &a)
Set matrix to be decomposed .
Double_t fPrecision
void SetThresholdPivoting(Double_t piv)
Int_t fIcntl[31]
virtual const TMatrixDBase & GetDecompMatrix() const
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.
virtual Bool_t Solve(TMatrixDColumn &)
virtual Int_t GetNcols() const
void SetVerbose(Int_t v)
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.
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.
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_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.
Double_t fIPessimism
virtual Bool_t TransSolve(TVectorD &b)
void InitParam()
initializing control parameters
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.
TMatrixDSparse fA
virtual TVectorD TransSolve(const TVectorD &b, Bool_t &ok)
Int_t fInfo[21]
virtual Bool_t TransSolve(TMatrixDColumn &)
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.
Double_t fCntl[6]
Int_t GetNrows() const
Int_t GetNcols() const
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:1016
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16