// @(#)root/matrix:$Id$
// 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