[u/mrichter/AliRoot.git] / STEER / AliSymMatrix.h

#ifndef ALISYMMATRIX_H
#define ALISYMMATRIX_H
/**********************************************************************************************/
/* Fast symmetric matrix with dynamically expandable size.                                    */
/* Only part can be used for matrix operations. It is defined as:                             */ 
/* fNCols: rows built by constructor (GetSizeBooked)                                          */ 
/* fNRows: number of rows added dynamically (automatically added on assignment to row)        */ 
/*         GetNRowAdded                                                                       */ 
/* fNRowIndex: total size (fNCols+fNRows), GetSize                                            */ 
/* fRowLwb   : actual size to used for given operation, by default = total size, GetSizeUsed  */ 
/*                                                                                            */ 
/* Author: ruben.shahoyan@cern.ch                                                             */
/*                                                                                            */ 
/**********************************************************************************************/

#include <TVectorD.h>
#include "AliMatrixSq.h"


class AliSymMatrix : public AliMatrixSq {
  //
 public:
  AliSymMatrix();
  AliSymMatrix(Int_t size);
  AliSymMatrix(const AliSymMatrix &mat);
  virtual ~AliSymMatrix();
  //
  void          Clear(Option_t* option="");
  void          Reset();
  //
  Int_t         GetSize()                                        const {return fNrowIndex;}
  Int_t         GetSizeUsed()                                    const {return fRowLwb;}
  Int_t         GetSizeBooked()                                  const {return fNcols;}
  Int_t         GetSizeAdded()                                   const {return fNrows;}
  Float_t       GetDensity()                                     const;
  AliSymMatrix& operator=(const AliSymMatrix& src);
  AliSymMatrix& operator+=(const AliSymMatrix& src);
  Double_t      operator()(Int_t rown, Int_t coln)               const;
  Double_t&     operator()(Int_t rown, Int_t coln);
  //
  Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymMatrix*)this)(r,r);}
  Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
  //
  Double_t*     GetRow(Int_t r);
  //
  void          Print(const Option_t* option="")                 const;
  void          AddRows(int nrows=1);
  void          SetSizeUsed(Int_t sz)                                  {fRowLwb = sz;}
  //
  void          Scale(Double_t coeff);
  void          MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
  void          MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
  void          AddToRow(Int_t r, Double_t *valc, Int_t *indc,Int_t n);
  //
  // ---------------------------------- Dummy methods of MatrixBase
  virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
  virtual             Double_t   *GetMatrixArray  ()       {return (Double_t*)fElems;}
  virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
  virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
  virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
  virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
  virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
  virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
  virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
  virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  //
  // ----------------------------- Choleski methods ----------------------------------------
  AliSymMatrix*       DecomposeChol();                  //Obtain Cholesky decomposition L matrix 
  void                InvertChol(AliSymMatrix* mchol);  //Invert using provided Choleski decomposition
  Bool_t              InvertChol();                     //Invert 
  Bool_t              SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
  Bool_t              SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert=kFALSE);
  Bool_t              SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
  Bool_t              SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
  //
  int                 SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);

 protected:
  virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
  Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
  void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
  //
 protected:
  Double_t*  fElems;     //   Elements booked by constructor
  Double_t** fElemsAdd;  //   Elements (rows) added dynamicaly
  //
  static AliSymMatrix* fgBuffer;  // buffer for fast solution
  static Int_t         fgCopyCnt;  // matrix copy counter
  ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
};


//___________________________________________________________
inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
{
  // lower triangle is actually filled
  return ((row*(row+1))>>1) + col;
}

//___________________________________________________________
inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
{
  //
  if (row<col) Swap(row,col);
  if (row>=fNrowIndex) return 0;
  return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
}

//___________________________________________________________
inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
{
  if (row<col) Swap(row,col);
  if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
  return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
}

//___________________________________________________________
inline void AliSymMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
{
  MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
}

//___________________________________________________________
inline void AliSymMatrix::Scale(Double_t coeff)
{
  for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
}

//___________________________________________________________
inline void AliSymMatrix::AddToRow(Int_t r, Double_t *valc, Int_t *indc,Int_t n)
{
  for (int i=n;i--;) (*this)(indc[i],r) += valc[i];
}

#endif
Commit	Line	Data
8a9ab0eb	1	#ifndef ALISYMMATRIX_H
8a9ab0eb	2	#define ALISYMMATRIX_H
7c3070ec	3	/**********************************************************************************************/
	4	/* Fast symmetric matrix with dynamically expandable size. */
	5	/* Only part can be used for matrix operations. It is defined as: */
	6	/* fNCols: rows built by constructor (GetSizeBooked) */
	7	/* fNRows: number of rows added dynamically (automatically added on assignment to row) */
	8	/* GetNRowAdded */
	9	/* fNRowIndex: total size (fNCols+fNRows), GetSize */
	10	/* fRowLwb : actual size to used for given operation, by default = total size, GetSizeUsed */
	11	/* */
	12	/* Author: ruben.shahoyan@cern.ch */
	13	/* */
	14	/**********************************************************************************************/
8a9ab0eb	15
8a9ab0eb	16	#include <TVectorD.h>
	17	#include "AliMatrixSq.h"
	18
	19
	20
	21	class AliSymMatrix : public AliMatrixSq {
	22	//
	23	public:
	24	AliSymMatrix();
	25	AliSymMatrix(Int_t size);
	26	AliSymMatrix(const AliSymMatrix &mat);
	27	virtual ~AliSymMatrix();
	28	//
	29	void Clear(Option_t* option="");
	30	void Reset();
	31	//
	32	Int_t GetSize() const {return fNrowIndex;}
7c3070ec	33	Int_t GetSizeUsed() const {return fRowLwb;}
	34	Int_t GetSizeBooked() const {return fNcols;}
	35	Int_t GetSizeAdded() const {return fNrows;}
8a9ab0eb	36	Float_t GetDensity() const;
8a9ab0eb	37	AliSymMatrix& operator=(const AliSymMatrix& src);
7c3070ec	38	AliSymMatrix& operator+=(const AliSymMatrix& src);
8a9ab0eb	39	Double_t operator()(Int_t rown, Int_t coln) const;
	40	Double_t& operator()(Int_t rown, Int_t coln);
	41	//
	42	Double_t DiagElem(Int_t r) const {return ((const AliSymMatrix)this)(r,r);}
	43	Double_t& DiagElem(Int_t r) {return (*this)(r,r);}
	44	//
de34b538	45	Double_t* GetRow(Int_t r);
de34b538	46	//
68f76645	47	void Print(const Option_t* option="") const;
8a9ab0eb	48	void AddRows(int nrows=1);
7c3070ec	49	void SetSizeUsed(Int_t sz) {fRowLwb = sz;}
8a9ab0eb	50	//
8a9ab0eb	51	void Scale(Double_t coeff);
551c9e69	52	void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
551c9e69	53	void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
68f76645	54	void AddToRow(Int_t r, Double_t valc, Int_t indc,Int_t n);
8a9ab0eb	55	//
	56	// ---------------------------------- Dummy methods of MatrixBase
	57	virtual const Double_t *GetMatrixArray () const {return fElems;};
	58	virtual Double_t GetMatrixArray () {return (Double_t)fElems;}
	59	virtual const Int_t *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
	60	virtual Int_t *GetRowIndexArray() {Error("GetRowIndexArray","Dummy"); return 0;};
	61	virtual const Int_t *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
	62	virtual Int_t *GetColIndexArray() {Error("GetColIndexArray","Dummy"); return 0;};
	63	virtual TMatrixDBase &SetRowIndexArray(Int_t ) {Error("SetRowIndexArray","Dummy"); return this;}
	64	virtual TMatrixDBase &SetColIndexArray(Int_t ) {Error("SetColIndexArray","Dummy"); return this;}
	65	virtual TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t ) const {Error("GetSub","Dummy"); return ((TMatrixDBase*)this);}
	66	virtual TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
	67	virtual TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
	68	virtual TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
	69	//
	70	// ----------------------------- Choleski methods ----------------------------------------
	71	AliSymMatrix* DecomposeChol(); //Obtain Cholesky decomposition L matrix
	72	void InvertChol(AliSymMatrix* mchol); //Invert using provided Choleski decomposition
	73	Bool_t InvertChol(); //Invert
	74	Bool_t SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
	75	Bool_t SolveChol(Double_t brhs, Double_t bsol,Bool_t invert=kFALSE);
	76	Bool_t SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
	77	Bool_t SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
	78	//
	79	int SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);
	80
	81	protected:
	82	virtual Int_t GetIndex(Int_t row,Int_t col) const;
	83	Double_t GetEl(Int_t row, Int_t col) const {return operator()(row,col);}
	84	void SetEl(Int_t row, Int_t col,Double_t val) {operator()(row,col) = val;}
	85	//
	86	protected:
	87	Double_t* fElems; // Elements booked by constructor
	88	Double_t** fElemsAdd; // Elements (rows) added dynamicaly
	89	//
	90	static AliSymMatrix* fgBuffer; // buffer for fast solution
	91	static Int_t fgCopyCnt; // matrix copy counter
	92	ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
	93	};
	94
	95
	96	//___________________________________________________________
	97	inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
	98	{
	99	// lower triangle is actually filled
	100	return ((row*(row+1))>>1) + col;
	101	}
	102
	103	//___________________________________________________________
	104	inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
	105	{
	106	//
	107	if (row<col) Swap(row,col);
	108	if (row>=fNrowIndex) return 0;
	109	return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
	110	}
	111
	112	//___________________________________________________________
	113	inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
	114	{
	115	if (row<col) Swap(row,col);
	116	if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
	117	return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
	118	}
119
120	//___________________________________________________________
551c9e69	121	inline void AliSymMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
8a9ab0eb	122	{
	123	MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
	124	}
	125
	126	//___________________________________________________________
	127	inline void AliSymMatrix::Scale(Double_t coeff)
	128	{
	129	for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
	130	}
	131
de34b538	132	//___________________________________________________________
68f76645	133	inline void AliSymMatrix::AddToRow(Int_t r, Double_t valc, Int_t indc,Int_t n)
de34b538	134	{
	135	for (int i=n;i--;) (*this)(indc[i],r) += valc[i];
	136	}
8a9ab0eb	137
8a9ab0eb	138	#endif