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

#ifndef ALISYMBDMATRIX_H
#define ALISYMBDMATRIX_H

/*********************************************************************************/
/* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal)       */
/* Only lower triangle is stored in the "profile" format                         */ 
/*                                                                               */ 
/*                                                                               */
/* Author: ruben.shahoyan@cern.ch                                                */
/*                                                                               */ 
/*********************************************************************************/

//#include <string.h>
#include <TObject.h>
#include <TVectorD.h>
#include "AliMatrixSq.h"


class AliSymBDMatrix : public AliMatrixSq {
  //
 public:
  enum {kDecomposedBit = 0x1};
  //
  AliSymBDMatrix();
  AliSymBDMatrix(Int_t size, Int_t w = 0);
  AliSymBDMatrix(const AliSymBDMatrix &mat);
  virtual ~AliSymBDMatrix();
  //
  Int_t         GetBandHWidth()                                  const {return fNrows;}
  Int_t         GetNElemsStored()                                const {return fNelems;}
  void          Clear(Option_t* option="");
  void          Reset();
  //
  Float_t       GetDensity()                                     const;
  AliSymBDMatrix& operator=(const AliSymBDMatrix& src);
  Double_t      operator()(Int_t rown, Int_t coln)               const;
  Double_t&     operator()(Int_t rown, Int_t coln);
  Double_t      operator()(Int_t rown)                           const;
  Double_t&     operator()(Int_t rown);
  //
  Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymBDMatrix*)this)(r,r);}
  Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
  void          DecomposeLDLT();
  void          Solve(Double_t *rhs);
  void          Solve(const Double_t *rhs,Double_t *sol);
  void          Solve(TVectorD &rhs)                                   {Solve(rhs.GetMatrixArray());}
  void          Solve(const TVectorD &rhs,TVectorD &sol)               {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());}
  //
  void          Print(Option_t* option="")                       const;
  void          SetDecomposed(Bool_t v=kTRUE)                          {SetBit(kDecomposedBit,v);}
  Bool_t        IsDecomposed()                                   const {return TestBit(kDecomposedBit);}
  //
  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);
  //
  // protected:
  virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
  virtual             Int_t GetIndex(Int_t diagID)             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
  //
  ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
};


//___________________________________________________________
inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
{
  // lower triangle band is actually filled
  if (row<col) Swap(row,col);
  col -= row;
  if (col < -GetBandHWidth()) return -1;
  return GetIndex(row) + col;
}

//___________________________________________________________
inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
{
  // Get index of the diagonal element on row diagID
  return (diagID+1)*fRowLwb-1;
}

//___________________________________________________________
inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
{
  // query element
  int idx = GetIndex(row,col);
  return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
}

//___________________________________________________________
inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
{
  // get element for assingment; assignment outside of the stored range has no effect
  int idx = GetIndex(row,col);  
  if (idx>=0) return fElems[idx];
  fTol = 0; 
  return fTol;
}

//___________________________________________________________
inline Double_t AliSymBDMatrix::operator()(Int_t row) const
{
  // query diagonal 
  return (const Double_t&) fElems[GetIndex(row)];
}

//___________________________________________________________
inline Double_t& AliSymBDMatrix::operator()(Int_t row)
{
  // get diagonal for assingment; assignment outside of the stored range has no effect
  return fElems[GetIndex(row)];
}

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

#endif
Commit	Line	Data
7c3070ec	1	#ifndef ALISYMBDMATRIX_H
	2	#define ALISYMBDMATRIX_H
	3
	4	/*********************************************************************************/
	5	/* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
	6	/* Only lower triangle is stored in the "profile" format */
	7	/* */
	8	/* */
	9	/* Author: ruben.shahoyan@cern.ch */
	10	/* */
	11	/*********************************************************************************/
	12
	13	//#include <string.h>
	14	#include <TObject.h>
	15	#include <TVectorD.h>
	16	#include "AliMatrixSq.h"
	17
	18
	19	class AliSymBDMatrix : public AliMatrixSq {
	20	//
	21	public:
	22	enum {kDecomposedBit = 0x1};
	23	//
	24	AliSymBDMatrix();
	25	AliSymBDMatrix(Int_t size, Int_t w = 0);
	26	AliSymBDMatrix(const AliSymBDMatrix &mat);
	27	virtual ~AliSymBDMatrix();
	28	//
	29	Int_t GetBandHWidth() const {return fNrows;}
	30	Int_t GetNElemsStored() const {return fNelems;}
	31	void Clear(Option_t* option="");
	32	void Reset();
	33	//
	34	Float_t GetDensity() const;
	35	AliSymBDMatrix& operator=(const AliSymBDMatrix& src);
	36	Double_t operator()(Int_t rown, Int_t coln) const;
	37	Double_t& operator()(Int_t rown, Int_t coln);
	38	Double_t operator()(Int_t rown) const;
	39	Double_t& operator()(Int_t rown);
	40	//
	41	Double_t DiagElem(Int_t r) const {return ((const AliSymBDMatrix)this)(r,r);}
	42	Double_t& DiagElem(Int_t r) {return (*this)(r,r);}
	43	void DecomposeLDLT();
	44	void Solve(Double_t *rhs);
	45	void Solve(const Double_t rhs,Double_t sol);
	46	void Solve(TVectorD &rhs) {Solve(rhs.GetMatrixArray());}
	47	void Solve(const TVectorD &rhs,TVectorD &sol) {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());}
	48	//
	49	void Print(Option_t* option="") const;
c130d912	50	void SetDecomposed(Bool_t v=kTRUE) {SetBit(kDecomposedBit,v);}
7c3070ec	51	Bool_t IsDecomposed() const {return TestBit(kDecomposedBit);}
7c3070ec	52	//
551c9e69	53	void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
551c9e69	54	void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
7c3070ec	55	void AddToRow(Int_t r, Double_t valc,Int_t indc,Int_t n);
	56	//
	57	// protected:
	58	virtual Int_t GetIndex(Int_t row,Int_t col) const;
	59	virtual Int_t GetIndex(Int_t diagID) const;
	60	Double_t GetEl(Int_t row, Int_t col) const {return operator()(row,col);}
	61	void SetEl(Int_t row, Int_t col,Double_t val) {operator()(row,col) = val;}
	62	//
	63	protected:
	64	Double_t* fElems; // Elements booked by constructor
	65	//
	66	ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
	67	};
	68
	69
	70	//___________________________________________________________
	71	inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
	72	{
	73	// lower triangle band is actually filled
	74	if (row<col) Swap(row,col);
	75	col -= row;
	76	if (col < -GetBandHWidth()) return -1;
	77	return GetIndex(row) + col;
	78	}
	79
	80	//___________________________________________________________
	81	inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
	82	{
	83	// Get index of the diagonal element on row diagID
	84	return (diagID+1)*fRowLwb-1;
	85	}
	86
	87	//___________________________________________________________
	88	inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
	89	{
	90	// query element
	91	int idx = GetIndex(row,col);
	92	return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
	93	}
	94
	95	//___________________________________________________________
	96	inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
	97	{
	98	// get element for assingment; assignment outside of the stored range has no effect
	99	int idx = GetIndex(row,col);
	100	if (idx>=0) return fElems[idx];
	101	fTol = 0;
	102	return fTol;
	103	}
	104
	105	//___________________________________________________________
	106	inline Double_t AliSymBDMatrix::operator()(Int_t row) const
	107	{
	108	// query diagonal
	109	return (const Double_t&) fElems[GetIndex(row)];
	110	}
	111
	112	//___________________________________________________________
	113	inline Double_t& AliSymBDMatrix::operator()(Int_t row)
	114	{
	115	// get diagonal for assingment; assignment outside of the stored range has no effect
	116	return fElems[GetIndex(row)];
	117	}
	118
119	//___________________________________________________________
551c9e69	120	inline void AliSymBDMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
7c3070ec	121	{
	122	MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
	123	}
	124
	125	#endif