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


/*********************************************************************************/
/* 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 <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <float.h>
//
#include "TClass.h"
#include "TMath.h"
#include "AliSymBDMatrix.h"
//

ClassImp(AliSymBDMatrix)


//___________________________________________________________
AliSymBDMatrix::AliSymBDMatrix() 
: fElems(0)
{
  // def. c-tor
  fSymmetric = kTRUE;
}

//___________________________________________________________
AliSymBDMatrix::AliSymBDMatrix(Int_t size, Int_t w)
  : AliMatrixSq(),fElems(0)
{
  // c-tor for given size
  //
  fNcols = size;    // number of rows
  if (w<0) w = 0;
  if (w>=size) w = size-1;
  fNrows = w;
  fRowLwb = w+1;
  fSymmetric = kTRUE;
  //
  // total number of stored elements
  fNelems = size*(w+1) - w*(w+1)/2;
  //
  fElems = new Double_t[fNcols*fRowLwb];
  memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
  //
}

//___________________________________________________________
AliSymBDMatrix::AliSymBDMatrix(const AliSymBDMatrix &src) 
  : AliMatrixSq(src),fElems(0)
{
  // copy c-tor
  if (src.GetSize()<1) return;
  fNcols = src.GetSize();
  fNrows = src.GetBandHWidth();
  fRowLwb = fNrows+1;
  fNelems = src.GetNElemsStored();
  fElems = new Double_t[fNcols*fRowLwb];
  memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
  //
}

//___________________________________________________________
AliSymBDMatrix::~AliSymBDMatrix() 
{
  // d-tor
  Clear();
}

//___________________________________________________________
AliSymBDMatrix&  AliSymBDMatrix::operator=(const AliSymBDMatrix& src)
{
  // assignment
  //
  if (this != &src) {
    TObject::operator=(src);
    if (fNcols!=src.fNcols) {
      // recreate the matrix
      if (fElems) delete[] fElems;
      fNcols = src.GetSize();
      fNrows = src.GetBandHWidth();
      fNelems = src.GetNElemsStored();
      fRowLwb = src.fRowLwb;
      fElems = new Double_t[fNcols*fRowLwb];
    }
    memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));     
    fSymmetric = kTRUE;
  }
  //
  return *this;
  //
}

//___________________________________________________________
void AliSymBDMatrix::Clear(Option_t*)
{
  // clear dynamic part
  if (fElems) {delete[] fElems; fElems = 0;}
  //  
  fNelems = fNcols = fNrows = fRowLwb = 0;
  //
}

//___________________________________________________________
Float_t AliSymBDMatrix::GetDensity() const
{
  // get fraction of non-zero elements
  if (!fNelems) return 0;
  Int_t nel = 0;
  for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
  return nel/fNelems;
}

//___________________________________________________________
void AliSymBDMatrix::Print(Option_t* option) const
{
  // print data
  printf("Symmetric Band-Diagonal Matrix : Size = %d, half bandwidth = %d\n",
	 GetSize(),GetBandHWidth());
  TString opt = option; opt.ToLower();
  if (opt.IsNull()) return;
  opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|";
  for (Int_t i=0;i<GetSize();i++) {
    printf(opt,i);
    for (Int_t j=TMath::Max(0,i-GetBandHWidth());j<=i;j++) printf("%+.3e|",GetEl(i,j));
    printf("\n");
  }
}

//___________________________________________________________
void AliSymBDMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const
{
  // fill vecOut by matrix*vecIn
  // vector should be of the same size as the matrix
  if (IsDecomposed()) {
    for (int i=0;i<GetSize();i++) {
      double sm = 0;
      int jmax = TMath::Min(GetSize(),i+fRowLwb);
      for (int j=i+1;j<jmax;j++) sm += vecIn[j]*Query(j,i);
      vecOut[i] = QueryDiag(i)*(vecIn[i]+sm);
    }
    for (int i=GetSize();i--;) {
      double sm = 0;
      int jmin = TMath::Max(0,i - GetBandHWidth());
      int jmax = i-1;
      for (int j=jmin;j<jmax;j++) sm += vecOut[j]*Query(i,j);
      vecOut[i] += sm;
    }
  }
  else { // not decomposed
    for (int i=GetSize();i--;) {
      vecOut[i] = 0.0;
      int jmin = TMath::Max(0,i - GetBandHWidth());
      int jmax = TMath::Min(GetSize(),i + fRowLwb);
      for (int j=jmin;j<jmax;j++) vecOut[i] += vecIn[j]*Query(i,j);
    }
  }
  //
}

//___________________________________________________________
void AliSymBDMatrix::Reset()
{
  // set all elems to 0
  if (fElems) memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
  SetDecomposed(kFALSE);
  //
}

//___________________________________________________________
void AliSymBDMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
{
  // add list of elements to row r
  for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
}

//___________________________________________________________
void AliSymBDMatrix::DecomposeLDLT()
{
  // decomposition to L Diag L^T
  if (IsDecomposed()) return;
  //
  Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
  //
  Double_t dtmp,gamma=0.0,xi=0.0;
  int iDiag;
  //
  // find max diag and number of non-0 diag.elements
  for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
    if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
    if (gamma < dtmp) gamma = dtmp; 
  }
  //
  // find max. off-diag element
  for (int ir=1;ir<iDiag;ir++) {
    for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
      if (ic<0) continue;
      dtmp = TMath::Abs(Query(ir,ic));
      if (xi<dtmp) xi = dtmp;
    }
  }
  double delta = eps*TMath::Max(1.0,xi+gamma);
  //
  double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
  double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn))); 
  //
  for (int kr=1;kr<GetSize();kr++) {
    int colKmin = TMath::Max(0,kr - GetBandHWidth());
    double theta = 0.0;
    //
    for (int jr=colKmin;jr<=kr;jr++) {
      int colJmin = TMath::Max(0,jr - GetBandHWidth());
      //
      dtmp = 0.0;
      for (int i=TMath::Max(colKmin,colJmin);i<jr;i++) 
	dtmp += Query(kr,i)*QueryDiag(i)*Query(jr,i);
      dtmp = (*this)(kr,jr) -= dtmp;
      //
      theta = TMath::Max(theta, TMath::Abs(dtmp));
      //
      if (jr!=kr) {
	if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
	else                        (*this)(kr,jr) = 0.0;
      }
      else if (kr<iDiag) {
	dtmp = theta/beta;
	dtmp *= dtmp;
	dtmp = TMath::Max(dtmp, delta);
	(*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp); 
      }
    } // jr
  } // kr
  //
  for (int i=0;i<GetSize();i++) {
    dtmp = QueryDiag(i);
    if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
  }
  //
  SetDecomposed();
}

//___________________________________________________________
void AliSymBDMatrix::Solve(Double_t *rhs)
{
  // solve matrix equation
  //
  if (!IsDecomposed()) DecomposeLDLT();
  //
  for (int kr=0;kr<GetSize();kr++) 
    for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
  //
  for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
  //
  for (int kr=GetSize();kr--;)
    for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
  //
}

//___________________________________________________________
void AliSymBDMatrix::Solve(const Double_t *rhs,Double_t *sol)
{
  // solve matrix equation
  memcpy(sol,rhs,GetSize()*sizeof(Double_t));
  Solve(sol);
  //
}
Commit	Line	Data
7c3070ec	1
	2	/*********************************************************************************/
	3	/* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
	4	/* Only lower triangle is stored in the "profile" format */
	5	/* */
	6	/* */
	7	/* Author: ruben.shahoyan@cern.ch */
	8	/* */
	9	/*********************************************************************************/
	10	#include <stdlib.h>
	11	#include <stdio.h>
	12	#include <iostream>
	13	#include <float.h>
	14	//
	15	#include "TClass.h"
	16	#include "TMath.h"
	17	#include "AliSymBDMatrix.h"
	18	//
	19
7c3070ec	20	ClassImp(AliSymBDMatrix)
	21
	22
	23	//___________________________________________________________
	24	AliSymBDMatrix::AliSymBDMatrix()
	25	: fElems(0)
	26	{
	27	// def. c-tor
	28	fSymmetric = kTRUE;
	29	}
	30
	31	//___________________________________________________________
	32	AliSymBDMatrix::AliSymBDMatrix(Int_t size, Int_t w)
	33	: AliMatrixSq(),fElems(0)
	34	{
	35	// c-tor for given size
	36	//
	37	fNcols = size; // number of rows
	38	if (w<0) w = 0;
	39	if (w>=size) w = size-1;
	40	fNrows = w;
	41	fRowLwb = w+1;
	42	fSymmetric = kTRUE;
	43	//
	44	// total number of stored elements
	45	fNelems = size(w+1) - w(w+1)/2;
	46	//
	47	fElems = new Double_t[fNcols*fRowLwb];
	48	memset(fElems,0,fNcolsfRowLwbsizeof(Double_t));
	49	//
	50	}
	51
	52	//___________________________________________________________
	53	AliSymBDMatrix::AliSymBDMatrix(const AliSymBDMatrix &src)
	54	: AliMatrixSq(src),fElems(0)
	55	{
	56	// copy c-tor
	57	if (src.GetSize()<1) return;
	58	fNcols = src.GetSize();
	59	fNrows = src.GetBandHWidth();
	60	fRowLwb = fNrows+1;
	61	fNelems = src.GetNElemsStored();
	62	fElems = new Double_t[fNcols*fRowLwb];
	63	memcpy(fElems,src.fElems,fNcolsfRowLwbsizeof(Double_t));
	64	//
	65	}
	66
	67	//___________________________________________________________
	68	AliSymBDMatrix::~AliSymBDMatrix()
	69	{
	70	// d-tor
	71	Clear();
	72	}
	73
	74	//___________________________________________________________
	75	AliSymBDMatrix& AliSymBDMatrix::operator=(const AliSymBDMatrix& src)
	76	{
	77	// assignment
	78	//
	79	if (this != &src) {
	80	TObject::operator=(src);
	81	if (fNcols!=src.fNcols) {
	82	// recreate the matrix
	83	if (fElems) delete[] fElems;
84	fNcols = src.GetSize();
85	fNrows = src.GetBandHWidth();
86	fNelems = src.GetNElemsStored();
87	fRowLwb = src.fRowLwb;
88	fElems = new Double_t[fNcols*fRowLwb];
89	}
90	memcpy(fElems,src.fElems,fNcolsfRowLwbsizeof(Double_t));
91	fSymmetric = kTRUE;
92	}
93	//
94	return *this;
95	//
96	}
97
98	//___________________________________________________________
99	void AliSymBDMatrix::Clear(Option_t*)
100	{
101	// clear dynamic part
102	if (fElems) {delete[] fElems; fElems = 0;}
103	//
104	fNelems = fNcols = fNrows = fRowLwb = 0;
105	//
106	}
107
108	//___________________________________________________________
109	Float_t AliSymBDMatrix::GetDensity() const
110	{
111	// get fraction of non-zero elements
112	if (!fNelems) return 0;
113	Int_t nel = 0;
7c185459	114	for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
7c3070ec	115	return nel/fNelems;
	116	}
	117
	118	//___________________________________________________________
	119	void AliSymBDMatrix::Print(Option_t* option) const
	120	{
	121	// print data
	122	printf("Symmetric Band-Diagonal Matrix : Size = %d, half bandwidth = %d\n",
	123	GetSize(),GetBandHWidth());
	124	TString opt = option; opt.ToLower();
	125	if (opt.IsNull()) return;
	126	opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d\|";
	127	for (Int_t i=0;i<GetSize();i++) {
	128	printf(opt,i);
	129	for (Int_t j=TMath::Max(0,i-GetBandHWidth());j<=i;j++) printf("%+.3e\|",GetEl(i,j));
	130	printf("\n");
	131	}
	132	}
	133
	134	//___________________________________________________________
551c9e69	135	void AliSymBDMatrix::MultiplyByVec(const Double_t vecIn,Double_t vecOut) const
7c3070ec	136	{
	137	// fill vecOut by matrix*vecIn
	138	// vector should be of the same size as the matrix
	139	if (IsDecomposed()) {
	140	for (int i=0;i<GetSize();i++) {
	141	double sm = 0;
	142	int jmax = TMath::Min(GetSize(),i+fRowLwb);
	143	for (int j=i+1;j<jmax;j++) sm += vecIn[j]*Query(j,i);
	144	vecOut[i] = QueryDiag(i)*(vecIn[i]+sm);
	145	}
	146	for (int i=GetSize();i--;) {
	147	double sm = 0;
	148	int jmin = TMath::Max(0,i - GetBandHWidth());
	149	int jmax = i-1;
	150	for (int j=jmin;j<jmax;j++) sm += vecOut[j]*Query(i,j);
	151	vecOut[i] += sm;
	152	}
	153	}
	154	else { // not decomposed
	155	for (int i=GetSize();i--;) {
	156	vecOut[i] = 0.0;
	157	int jmin = TMath::Max(0,i - GetBandHWidth());
	158	int jmax = TMath::Min(GetSize(),i + fRowLwb);
	159	for (int j=jmin;j<jmax;j++) vecOut[i] += vecIn[j]*Query(i,j);
	160	}
	161	}
	162	//
	163	}
	164
	165	//___________________________________________________________
	166	void AliSymBDMatrix::Reset()
	167	{
	168	// set all elems to 0
	169	if (fElems) memset(fElems,0,fNcolsfRowLwbsizeof(Double_t));
c130d912	170	SetDecomposed(kFALSE);
7c3070ec	171	//
	172	}
	173
	174	//___________________________________________________________
	175	void AliSymBDMatrix::AddToRow(Int_t r, Double_t valc,Int_t indc,Int_t n)
	176	{
	177	// add list of elements to row r
	178	for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
	179	}
	180
	181	//___________________________________________________________
	182	void AliSymBDMatrix::DecomposeLDLT()
	183	{
	184	// decomposition to L Diag L^T
	185	if (IsDecomposed()) return;
	186	//
	187	Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
	188	//
	189	Double_t dtmp,gamma=0.0,xi=0.0;
	190	int iDiag;
	191	//
	192	// find max diag and number of non-0 diag.elements
	193	for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
	194	if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
	195	if (gamma < dtmp) gamma = dtmp;
	196	}
	197	//
	198	// find max. off-diag element
	199	for (int ir=1;ir<iDiag;ir++) {
	200	for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
	201	if (ic<0) continue;
	202	dtmp = TMath::Abs(Query(ir,ic));
	203	if (xi<dtmp) xi = dtmp;
	204	}
	205	}
	206	double delta = eps*TMath::Max(1.0,xi+gamma);
	207	//
	208	double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
	209	double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn)));
	210	//
	211	for (int kr=1;kr<GetSize();kr++) {
	212	int colKmin = TMath::Max(0,kr - GetBandHWidth());
	213	double theta = 0.0;
	214	//
	215	for (int jr=colKmin;jr<=kr;jr++) {
	216	int colJmin = TMath::Max(0,jr - GetBandHWidth());
	217	//
	218	dtmp = 0.0;
	219	for (int i=TMath::Max(colKmin,colJmin);i<jr;i++)
	220	dtmp += Query(kr,i)QueryDiag(i)Query(jr,i);
	221	dtmp = (*this)(kr,jr) -= dtmp;
	222	//
	223	theta = TMath::Max(theta, TMath::Abs(dtmp));
	224	//
	225	if (jr!=kr) {
7c185459	226	if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
7c185459	227	else (*this)(kr,jr) = 0.0;
7c3070ec	228	}
	229	else if (kr<iDiag) {
	230	dtmp = theta/beta;
	231	dtmp *= dtmp;
	232	dtmp = TMath::Max(dtmp, delta);
	233	(*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp);
	234	}
	235	} // jr
	236	} // kr
	237	//
	238	for (int i=0;i<GetSize();i++) {
	239	dtmp = QueryDiag(i);
7c185459	240	if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
7c3070ec	241	}
	242	//
	243	SetDecomposed();
	244	}
	245
	246	//___________________________________________________________
	247	void AliSymBDMatrix::Solve(Double_t *rhs)
	248	{
	249	// solve matrix equation
	250	//
	251	if (!IsDecomposed()) DecomposeLDLT();
	252	//
	253	for (int kr=0;kr<GetSize();kr++)
	254	for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
	255	//
	256	for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
	257	//
	258	for (int kr=GetSize();kr--;)
	259	for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
	260	//
	261	}
	262
	263	//___________________________________________________________
	264	void AliSymBDMatrix::Solve(const Double_t rhs,Double_t sol)
	265	{
	266	// solve matrix equation
	267	memcpy(sol,rhs,GetSize()*sizeof(Double_t));
	268	Solve(sol);
	269	//
	270	}