// @(#)root/matrix:$Name: $:$Id: TDecompSparse.h,v 1.4 2004/10/16 18:09:16 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann Apr 2004
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
#ifndef ROOT_TDecompSparse
#define ROOT_TDecompSparse
///////////////////////////////////////////////////////////////////////////
// //
// Sparse Decomposition class //
// //
///////////////////////////////////////////////////////////////////////////
#ifndef ROOT_TDecompBase
#include "TDecompBase.h"
#endif
#ifndef ROOT_TMatrixDSparse
#include "TMatrixDSparse.h"
#endif
#ifndef ROOT_TArrayD
#include "TArrayD.h"
#endif
#ifndef ROOT_TArrayI
#include "TArrayI.h"
#endif
// globals
const Double_t kInitTreatAsZero = 1.0e-12;
const Double_t kInitThresholdPivoting = 1.0e-8;
const Double_t kInitPrecision = 1.0e-7;
// the Threshold Pivoting parameter may need to be increased during
// the algorithm if poor precision is obtained from the linear
// solves. kThresholdPivoting indicates the largest value we are
// willing to tolerate.
const Double_t kThresholdPivotingMax = 1.0e-2;
// the factor in the range (1,inf) by which kThresholdPivoting is
// increased when it is found to be inadequate.
const Double_t kThresholdPivotingFactor = 10.0;
class TDecompSparse : public TDecompBase
{
protected :
Int_t fVerbose;
Int_t fIcntl[31]; // integer control numbers
Double_t fCntl[6]; // float control numbers
Int_t fInfo[21]; // array used for communication between programs
Double_t fPrecision; // precision we demand from the linear system solver. If it isn't
// attained on the first solve, we use iterative refinement and
// possibly refactorization with a higher value of kThresholdPivoting.
TArrayI fIkeep; // pivot sequence and temporary storage information
TArrayI fIw;
TArrayI fIw1;
TArrayI fIw2;
Int_t fNsteps;
Int_t fMaxfrt;
TArrayD fW; // temporary storage for the factorization
Double_t fIPessimism; // amounts by which to increase allocated factorization space when
Double_t fRPessimism; // inadequate space is detected. fIPessimism is for array "fIw",
// fRPessimism is for the array "fact".
TMatrixDSparse fA; // original matrix; needed for the iterative solving procedure
Int_t fNrows;
Int_t fNnonZeros;
TArrayD fFact; // size of fFact array; may be increased during the numerical factorization
// if the estimate obtained during the symbolic factorization proves to be inadequate.
TArrayI fRowFact;
TArrayI fColFact;
static Int_t NonZerosUpperTriang(const TMatrixDSparse &a);
static void CopyUpperTriang (const TMatrixDSparse &a,Double_t *b);
void InitParam();
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);
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);
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);
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);
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);
static void InitPivot_sub2a(const Int_t n,Int_t *ipe,Int_t *iw,const Int_t lw,Int_t &iwfr,Int_t &ncmpa);
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);
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);
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);
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);
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);
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);
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);
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);
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);
static Int_t IDiag (Int_t ix,Int_t iy) { return ((iy-1)*(2*ix-iy+2))/2; }
inline Int_t IError () { return fInfo[2]; }
inline Int_t MinRealWorkspace() { return fInfo[5]; }
inline Int_t MinIntWorkspace () { return fInfo[6]; }
inline Int_t ErrorFlag () { return fInfo[1]; }
// Takes values in the range [0,1]. Larger values enforce greater stability in
// the factorization as they insist on larger pivots. Smaller values preserve
// sparsity at the cost of using smaller pivots.
inline Double_t GetThresholdPivoting() { return fCntl[1]; }
inline Double_t GetTreatAsZero () { return fCntl[3]; }
// The factorization will not accept a pivot whose absolute value is less than fCntl[3] as
// a 1x1 pivot or as the off-diagonal in a 2x2 pivot.
inline void SetThresholdPivoting(Double_t piv) { fCntl[1] = piv; }
inline void SetTreatAsZero (Double_t tol) { fCntl[3] = tol; }
virtual const TMatrixDBase &GetDecompMatrix() const { MayNotUse("GetDecompMatrix()"); return fA; }
public :
TDecompSparse();
TDecompSparse(Int_t nRows,Int_t nr_nonZeros,Int_t verbose);
TDecompSparse(Int_t row_lwb,Int_t row_upb,Int_t nr_nonZeros,Int_t verbose);
TDecompSparse(const TMatrixDSparse &a,Int_t verbose);
TDecompSparse(const TDecompSparse &another);
virtual ~TDecompSparse() {}
inline void SetVerbose (Int_t v) { fVerbose = (v) ? 1 : 0;
if (fVerbose) { fIcntl[1] = fIcntl[2] = 1; fIcntl[3] = 2; }
else { fIcntl[1] = fIcntl[2] = fIcntl[3] = 0; }
}
virtual Int_t GetNrows () const { return fA.GetNrows(); }
virtual Int_t GetNcols () const { return fA.GetNcols(); }
virtual void SetMatrix (const TMatrixDSparse &a);
virtual Bool_t Decompose ();
virtual Bool_t Solve ( TVectorD &b);
virtual TVectorD Solve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
virtual Bool_t Solve ( TMatrixDColumn & /*b*/)
{ MayNotUse("Solve(TMatrixDColumn &)"); return kFALSE; }
virtual Bool_t TransSolve ( TVectorD &b) { return Solve(b); }
virtual TVectorD TransSolve (const TVectorD& b,Bool_t &ok) { TVectorD x = b; ok = Solve(x); return x; }
virtual Bool_t TransSolve ( TMatrixDColumn & /*b*/)
{ MayNotUse("TransSolve(TMatrixDColumn &)"); return kFALSE; }
virtual void Det (Double_t &/*d1*/,Double_t &/*d2*/)
{ MayNotUse("Det(Double_t&,Double_t&)"); }
void Print(Option_t *opt ="") const; // *MENU*
TDecompSparse &operator= (const TDecompSparse &source);
ClassDef(TDecompSparse,1) // Matrix Decompositition LU
};
#endif
syntax highlighted by Code2HTML, v. 0.9.1