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

#include <stdlib.h>
#include <stdio.h>
#include <iostream>
//
#include "TClass.h"
#include "TMath.h"
#include "AliSymMatrix.h"
//

using namespace std;

ClassImp(AliSymMatrix)


AliSymMatrix* AliSymMatrix::fgBuffer = 0; 
Int_t         AliSymMatrix::fgCopyCnt = 0; 
//___________________________________________________________
AliSymMatrix::AliSymMatrix() 
: fElems(0),fElemsAdd(0)
{
  fSymmetric = kTRUE;
  fgCopyCnt++;
}

//___________________________________________________________
AliSymMatrix::AliSymMatrix(Int_t size)
  : AliMatrixSq(),fElems(0),fElemsAdd(0)
{
  //
  fNrows = 0;
  fNrowIndex = fNcols = size;
  fElems     = new Double_t[fNcols*(fNcols+1)/2];
  fSymmetric = kTRUE;
  Reset();
  fgCopyCnt++;
  //
}

//___________________________________________________________
AliSymMatrix::AliSymMatrix(const AliSymMatrix &src) 
  : AliMatrixSq(src),fElems(0),fElemsAdd(0)
{
  fNrowIndex = fNcols = src.GetSize();
  fNrows = 0;
  if (fNcols) {
    int nmainel = fNcols*(fNcols+1)/2;
    fElems     = new Double_t[nmainel];
    nmainel = src.fNcols*(src.fNcols+1)/2;
    memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
    if (src.fNrows) { // transfer extra rows to main matrix
      Double_t *pnt = fElems + nmainel;
      int ncl = src.fNcols + 1;
      for (int ir=0;ir<src.fNrows;ir++) {
	memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
	pnt += ncl;
	ncl++; 
      }
    }
  }
  else fElems = 0;
  fElemsAdd = 0;
  fgCopyCnt++;
  //
}

//___________________________________________________________
AliSymMatrix::~AliSymMatrix() 
{
  Clear();
  if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;}
}

//___________________________________________________________
AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
{
  //
  if (this != &src) {
    TObject::operator=(src);
    if (fNcols!=src.fNcols && fNrows!=src.fNrows) {
      // recreate the matrix
      if (fElems) delete[] fElems;
      for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
      delete[] fElemsAdd;
      //
      fNrowIndex = src.GetSize(); 
      fNcols = src.GetSize();
      fNrows = 0;
      fElems     = new Double_t[fNcols*(fNcols+1)/2];
      int nmainel = src.fNcols*(src.fNcols+1);
      memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
      if (src.fNrows) { // transfer extra rows to main matrix
	Double_t *pnt = fElems + nmainel*sizeof(Double_t);
	int ncl = src.fNcols + 1;
	for (int ir=0;ir<src.fNrows;ir++) {
	  ncl += ir; 
	  memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
	  pnt += ncl*sizeof(Double_t);
	}
      }
      //
    }
    else {
      memcpy(fElems,src.fElems,fNcols*(fNcols+1)/2*sizeof(Double_t));
      int ncl = fNcols + 1;
      for (int ir=0;ir<fNrows;ir++) { // dynamic rows
	ncl += ir; 
	memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
      }
    }
  }
  //
  return *this;
}

//___________________________________________________________
void AliSymMatrix::Clear(Option_t*)
{
  if (fElems) {delete[] fElems; fElems = 0;}
  //  
  if (fElemsAdd) {
    for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
    delete[] fElemsAdd;
    fElemsAdd = 0;
  }
  fNrowIndex = 0;
  fNcols = 0;
  fNrows = 0;
  //
}

//___________________________________________________________
Float_t AliSymMatrix::GetDensity() const
{
  // get fraction of non-zero elements
  Int_t nel = 0;
  for (int i=GetSize();i--;) for (int j=i+1;j--;) if (GetEl(i,j)!=0) nel++;
  return 2.*nel/( (GetSize()+1)*GetSize() );
}

//___________________________________________________________
void AliSymMatrix::Print(Option_t* option) const
{
  printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
  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<fNrowIndex;i++) {
    printf(opt,i);
    for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j));
    printf("\n");
  }
}

//___________________________________________________________
void AliSymMatrix::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
{
  // fill vecOut by matrix*vecIn
  // vector should be of the same size as the matrix
  for (int i=fNrowIndex;i--;) {
    vecOut[i] = 0.0;
    for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
  }
  //
}

//___________________________________________________________
AliSymMatrix* AliSymMatrix::DecomposeChol() 
{
  // Return a matrix with Choleski decomposition
  // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
  // consturcts Cholesky decomposition of SYMMETRIC and
  // POSITIVELY-DEFINED matrix a (a=L*Lt)
  // Only upper triangle of the matrix has to be filled.
  // In opposite to function from the book, the matrix is modified:
  // lower triangle and diagonal are refilled.
  //
  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
    delete fgBuffer; 
    try {
      fgBuffer = new AliSymMatrix(*this);
    }
    catch(bad_alloc&) {
      printf("Failed to allocate memory for Choleski decompostions\n");
      return 0;
    }
  }
  else (*fgBuffer) = *this;
  //
  AliSymMatrix& mchol = *fgBuffer;
  //
  for (int i=0;i<fNrowIndex;i++) {
    Double_t *rowi = mchol.GetRow(i);
    for (int j=i;j<fNrowIndex;j++) {
      Double_t *rowj = mchol.GetRow(j);
      double sum = rowj[i];
      for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k];
      if (i == j) {
	if (sum <= 0.0) { // not positive-definite
	  printf("The matrix is not positive definite [%e]\n"
		 "Choleski decomposition is not possible\n",sum);
	  return 0;
	}
	rowi[i] = TMath::Sqrt(sum);
	//
      } else rowj[i] = sum/rowi[i];
    }
  }
  return fgBuffer;
}

//___________________________________________________________
Bool_t AliSymMatrix::InvertChol() 
{
  // Invert matrix using Choleski decomposition
  //
  AliSymMatrix* mchol = DecomposeChol();
  if (!mchol) {
    printf("Failed to invert the matrix\n");
    return kFALSE;
  }
  //
  InvertChol(mchol);
  return kTRUE;
  //
}

//___________________________________________________________
void AliSymMatrix::InvertChol(AliSymMatrix* pmchol) 
{
  // Invert matrix using Choleski decomposition, provided the Cholseki's L matrix
  //
  Double_t sum;
  AliSymMatrix& mchol = *pmchol;
  //
  // Invert decomposed triangular L matrix (Lower triangle is filled)
  for (int i=0;i<fNrowIndex;i++) { 
    mchol(i,i) =  1.0/mchol(i,i);
    for (int j=i+1;j<fNrowIndex;j++) { 
      Double_t *rowj = mchol.GetRow(j);
      sum = 0.0; 
      for (int k=i;k<j;k++) if (rowj[k]) { 
	double &mki = mchol(k,i); if (mki) sum -= rowj[k]*mki;
      }
      rowj[i] = sum/rowj[j];
    } 
  }
  //
  // take product of the inverted Choleski L matrix with its transposed
  for (int i=fNrowIndex;i--;) {
    for (int j=i+1;j--;) {
      sum = 0;
      for (int k=i;k<fNrowIndex;k++) {
	double &mik = mchol(i,k); 
	if (mik) {
	  double &mjk = mchol(j,k);
	  if (mjk) sum += mik*mjk;
	}
      }
      (*this)(j,i) = sum;
    }
  }
  //
}


//___________________________________________________________
Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert) 
{
  // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
  // Solves the set of n linear equations A x = b, 
  // where a is a positive-definite symmetric matrix. 
  // a[1..n][1..n] is the output of the routine CholDecomposw. 
  // Only the lower triangle of a is accessed. b[1..n] is input as the 
  // right-hand side vector. The solution vector is returned in b[1..n].
  //
  Int_t i,k;
  Double_t sum;
  //
  AliSymMatrix *pmchol = DecomposeChol();
  if (!pmchol) {
    printf("SolveChol failed\n");
    //    Print("l");
    return kFALSE;
  }
  AliSymMatrix& mchol = *pmchol;
  //
  for (i=0;i<fNrowIndex;i++) {
    Double_t *rowi = mchol.GetRow(i);
    for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k];
    b[i]=sum/rowi[i];
  }
  //
  for (i=fNrowIndex-1;i>=0;i--) {
    for (sum=b[i],k=i+1;k<fNrowIndex;k++) if (b[k]) {
      double &mki=mchol(k,i); if (mki) sum -= mki*b[k];
    }
    b[i]=sum/mchol(i,i);
  }
  //
  if (invert) InvertChol(pmchol);
  return kTRUE;
  //
}

//___________________________________________________________
Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert) 
{
  return SolveChol((Double_t*)b.GetMatrixArray(),invert);
}


//___________________________________________________________
Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert) 
{
  memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
  return SolveChol(bsol,invert);  
}

//___________________________________________________________
Bool_t AliSymMatrix::SolveChol(TVectorD &brhs, TVectorD &bsol,Bool_t invert) 
{
  bsol = brhs;
  return SolveChol(bsol,invert);
}

//___________________________________________________________
void AliSymMatrix::AddRows(int nrows)
{
  if (nrows<1) return;
  Double_t **pnew = new Double_t*[nrows+fNrows];
  for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows
  for (int ir=0;ir<nrows;ir++) {
    int ncl = GetSize()+1;
    pnew[fNrows] = new Double_t[ncl];
    memset(pnew[fNrows],0,ncl*sizeof(Double_t));
    fNrows++;
    fNrowIndex++;
  }
  delete[] fElemsAdd;
  fElemsAdd = pnew;
  //
}

//___________________________________________________________
void AliSymMatrix::Reset()
{
  // if additional rows exist, regularize it
  if (fElemsAdd) {
    delete[] fElems;
    for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
    delete[] fElemsAdd; fElemsAdd = 0;
    fNcols = fNrowIndex;
    fElems = new Double_t[fNcols*(fNcols+1)/2];
    fNrows = 0;
  }
  if (fElems) memset(fElems,0,fNcols*(fNcols+1)/2*sizeof(Double_t));
  //
}

//___________________________________________________________
/*
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];
  //   }
  //   return;

  double *row;
  if (r>=fNrowIndex) {
    AddRows(r-fNrowIndex+1); 
    row = &((fElemsAdd[r-fNcols])[0]);
  }
  else row = &fElems[GetIndex(r,0)];
  //
  int nadd = 0;
  for (int i=n;i--;) {
    if (indc[i]>r) continue;
    row[indc[i]] += valc[i];
    nadd++;
  }
  if (nadd == n) return;
  //
  // add to col>row
  for (int i=n;i--;) {
    if (indc[i]>r) (*this)(indc[i],r) += valc[i];
  }
  //
}
*/

//___________________________________________________________
Double_t* AliSymMatrix::GetRow(Int_t r)
{
  if (r>=fNrowIndex) {
    int nn = fNrowIndex;
    AddRows(r-fNrowIndex+1); 
    printf("create %d of %d\n",r, nn);
    return &((fElemsAdd[r-fNcols])[0]);
  }
  else return &fElems[GetIndex(r,0)];
}


//___________________________________________________________
int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
{
  //   Solution a la MP1: gaussian eliminations
  ///  Obtain solution of a system of linear equations with symmetric matrix 
  ///  and the inverse (using 'singular-value friendly' GAUSS pivot)
  //

  Int_t nRank = 0;
  int iPivot;
  double vPivot = 0.;
  double eps = 0.00000000000001;
  int nGlo = GetSize();
  bool   *bUnUsed = new bool[nGlo];
  double *rowMax,*colMax=0;
  rowMax  = new double[nGlo];
  //
  if (stabilize) {
    colMax   = new double[nGlo];
    for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0;
    for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { 
	double vl = TMath::Abs(Query(i,j));
	if (vl==0) continue;
	if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i
	if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j
	if (i==j) continue;
	if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j
	if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i
      }
    //
    for (Int_t i=nGlo; i--;) {
      if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
      if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
    }
    //
  }
  //
  for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
  //  
  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
    delete fgBuffer; 
    try {
      fgBuffer = new AliSymMatrix(*this);
    }
    catch(bad_alloc&) {
      printf("Failed to allocate memory for matrix inversion buffer\n");
      return 0;
    }
  }
  else (*fgBuffer) = *this;
  //
  if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) 
      for (int j=0;j<=i; j++) {
	double vl = Query(i,j);
	if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
      }
      for (int j=i+1;j<nGlo;j++) {
	double vl = Query(j,i);
	if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
      }
    }
  //
  for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values 
  //
  for (Int_t i=0; i<nGlo; i++) {
    vPivot = 0.0;
    iPivot = -1;
    //
    for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element       
      double vl;
      if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) {    
	vPivot = vl;
	iPivot = j;
      }
    }
    //
    if (iPivot >= 0) {   // pivot found          
      nRank++;
      bUnUsed[iPivot] = false; // This value is used
      vPivot = 1.0/vPivot;
      DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse
      //
      for (Int_t j=0; j<nGlo; j++) {      
	for (Int_t jj=0; jj<nGlo; jj++) {  
	  if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)	  
	    double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j);
	    r -= vPivot* ( j>iPivot  ? Query(j,iPivot)  : fgBuffer->Query(iPivot,j) )
	      *          ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot));
	  }
	}
      }
      //
      for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements 
	  (*this)(j,iPivot)     *= vPivot;
	  (*fgBuffer)(iPivot,j) *= vPivot;
	}
      //
    }
    else {  // No more pivot value (clear those elements)
      for (Int_t j=0; j<nGlo; j++) {
	if (bUnUsed[j]) {
	  vecB[j] = 0.0;
	  for (Int_t k=0; k<nGlo; k++) {
	    (*this)(j,k) = 0.;
	    if (j!=k) (*fgBuffer)(j,k) = 0;
	  }
	}
      }
      break;  // No more pivots anyway, stop here
    }
  }
  //
  for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) {
      double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V
      if (i>=j) (*this)(i,j) *= vl;
      else      (*fgBuffer)(j,i) *= vl;
    }
  //
  for (Int_t j=0; j<nGlo; j++) {
    rowMax[j] = 0.0;
    for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
      double vl;
      if (j>=jj) vl = (*this)(j,jj)     = -Query(j,jj);
      else       vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj);
      rowMax[j] += vl*vecB[jj];
    }		
  }
  
  for (Int_t j=0; j<nGlo; j++) {
    vecB[j] = rowMax[j]; // The final result
  }
  //
  delete [] bUnUsed;
  delete [] rowMax;
  if (stabilize) delete [] colMax;

  return nRank;
}
Commit	Line	Data
8a9ab0eb	1	#include <stdlib.h>
	2	#include <stdio.h>
	3	#include <iostream>
	4	//
	5	#include "TClass.h"
	6	#include "TMath.h"
	7	#include "AliSymMatrix.h"
	8	//
	9
	10	using namespace std;
	11
	12	ClassImp(AliSymMatrix)
	13
	14
	15	AliSymMatrix* AliSymMatrix::fgBuffer = 0;
	16	Int_t AliSymMatrix::fgCopyCnt = 0;
	17	//___________________________________________________________
	18	AliSymMatrix::AliSymMatrix()
	19	: fElems(0),fElemsAdd(0)
	20	{
	21	fSymmetric = kTRUE;
	22	fgCopyCnt++;
	23	}
	24
	25	//___________________________________________________________
	26	AliSymMatrix::AliSymMatrix(Int_t size)
	27	: AliMatrixSq(),fElems(0),fElemsAdd(0)
	28	{
	29	//
	30	fNrows = 0;
	31	fNrowIndex = fNcols = size;
	32	fElems = new Double_t[fNcols*(fNcols+1)/2];
	33	fSymmetric = kTRUE;
	34	Reset();
	35	fgCopyCnt++;
	36	//
	37	}
	38
	39	//___________________________________________________________
	40	AliSymMatrix::AliSymMatrix(const AliSymMatrix &src)
	41	: AliMatrixSq(src),fElems(0),fElemsAdd(0)
	42	{
	43	fNrowIndex = fNcols = src.GetSize();
	44	fNrows = 0;
	45	if (fNcols) {
	46	int nmainel = fNcols*(fNcols+1)/2;
	47	fElems = new Double_t[nmainel];
	48	nmainel = src.fNcols*(src.fNcols+1)/2;
	49	memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
	50	if (src.fNrows) { // transfer extra rows to main matrix
	51	Double_t *pnt = fElems + nmainel;
	52	int ncl = src.fNcols + 1;
	53	for (int ir=0;ir<src.fNrows;ir++) {
	54	memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
	55	pnt += ncl;
	56	ncl++;
	57	}
	58	}
	59	}
	60	else fElems = 0;
	61	fElemsAdd = 0;
	62	fgCopyCnt++;
	63	//
	64	}
65
66	//___________________________________________________________
67	AliSymMatrix::~AliSymMatrix()
68	{
69	Clear();
70	if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;}
71	}
72
73	//___________________________________________________________
74	AliSymMatrix& AliSymMatrix::operator=(const AliSymMatrix& src)
75	{
76	//
77	if (this != &src) {
78	TObject::operator=(src);
79	if (fNcols!=src.fNcols && fNrows!=src.fNrows) {
80	// recreate the matrix
81	if (fElems) delete[] fElems;
82	for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
83	delete[] fElemsAdd;
84	//
de34b538	85	fNrowIndex = src.GetSize();
de34b538	86	fNcols = src.GetSize();
8a9ab0eb	87	fNrows = 0;
	88	fElems = new Double_t[fNcols*(fNcols+1)/2];
	89	int nmainel = src.fNcols*(src.fNcols+1);
	90	memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
	91	if (src.fNrows) { // transfer extra rows to main matrix
	92	Double_t pnt = fElems + nmainelsizeof(Double_t);
	93	int ncl = src.fNcols + 1;
	94	for (int ir=0;ir<src.fNrows;ir++) {
	95	ncl += ir;
	96	memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
	97	pnt += ncl*sizeof(Double_t);
	98	}
	99	}
	100	//
	101	}
	102	else {
	103	memcpy(fElems,src.fElems,fNcols(fNcols+1)/2sizeof(Double_t));
	104	int ncl = fNcols + 1;
	105	for (int ir=0;ir<fNrows;ir++) { // dynamic rows
	106	ncl += ir;
	107	memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
	108	}
	109	}
	110	}
	111	//
	112	return *this;
	113	}
	114
	115	//___________________________________________________________
	116	void AliSymMatrix::Clear(Option_t*)
	117	{
	118	if (fElems) {delete[] fElems; fElems = 0;}
	119	//
	120	if (fElemsAdd) {
	121	for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
	122	delete[] fElemsAdd;
	123	fElemsAdd = 0;
	124	}
de34b538	125	fNrowIndex = 0;
	126	fNcols = 0;
	127	fNrows = 0;
8a9ab0eb	128	//
	129	}
	130
	131	//___________________________________________________________
	132	Float_t AliSymMatrix::GetDensity() const
	133	{
	134	// get fraction of non-zero elements
	135	Int_t nel = 0;
	136	for (int i=GetSize();i--;) for (int j=i+1;j--;) if (GetEl(i,j)!=0) nel++;
	137	return 2.nel/( (GetSize()+1)GetSize() );
	138	}
	139
	140	//___________________________________________________________
	141	void AliSymMatrix::Print(Option_t* option) const
	142	{
	143	printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
	144	TString opt = option; opt.ToLower();
	145	if (opt.IsNull()) return;
	146	opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d\|";
	147	for (Int_t i=0;i<fNrowIndex;i++) {
	148	printf(opt,i);
	149	for (Int_t j=0;j<=i;j++) printf("%+.3e\|",GetEl(i,j));
	150	printf("\n");
	151	}
	152	}
	153
	154	//___________________________________________________________
	155	void AliSymMatrix::MultiplyByVec(Double_t vecIn,Double_t vecOut) const
	156	{
	157	// fill vecOut by matrix*vecIn
	158	// vector should be of the same size as the matrix
	159	for (int i=fNrowIndex;i--;) {
	160	vecOut[i] = 0.0;
	161	for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
	162	}
	163	//
	164	}
	165
	166	//___________________________________________________________
	167	AliSymMatrix* AliSymMatrix::DecomposeChol()
	168	{
	169	// Return a matrix with Choleski decomposition
de34b538	170	// Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
	171	// consturcts Cholesky decomposition of SYMMETRIC and
	172	// POSITIVELY-DEFINED matrix a (a=L*Lt)
	173	// Only upper triangle of the matrix has to be filled.
	174	// In opposite to function from the book, the matrix is modified:
	175	// lower triangle and diagonal are refilled.
8a9ab0eb	176	//
	177	if (!fgBuffer \|\| fgBuffer->GetSize()!=GetSize()) {
	178	delete fgBuffer;
	179	try {
	180	fgBuffer = new AliSymMatrix(*this);
	181	}
	182	catch(bad_alloc&) {
	183	printf("Failed to allocate memory for Choleski decompostions\n");
	184	return 0;
	185	}
	186	}
	187	else (fgBuffer) = this;
	188	//
	189	AliSymMatrix& mchol = *fgBuffer;
	190	//
	191	for (int i=0;i<fNrowIndex;i++) {
de34b538	192	Double_t *rowi = mchol.GetRow(i);
8a9ab0eb	193	for (int j=i;j<fNrowIndex;j++) {
de34b538	194	Double_t *rowj = mchol.GetRow(j);
	195	double sum = rowj[i];
	196	for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k];
8a9ab0eb	197	if (i == j) {
	198	if (sum <= 0.0) { // not positive-definite
	199	printf("The matrix is not positive definite [%e]\n"
	200	"Choleski decomposition is not possible\n",sum);
	201	return 0;
	202	}
de34b538	203	rowi[i] = TMath::Sqrt(sum);
8a9ab0eb	204	//
de34b538	205	} else rowj[i] = sum/rowi[i];
8a9ab0eb	206	}
	207	}
	208	return fgBuffer;
	209	}
	210
	211	//___________________________________________________________
	212	Bool_t AliSymMatrix::InvertChol()
	213	{
	214	// Invert matrix using Choleski decomposition
	215	//
	216	AliSymMatrix* mchol = DecomposeChol();
	217	if (!mchol) {
	218	printf("Failed to invert the matrix\n");
	219	return kFALSE;
	220	}
	221	//
	222	InvertChol(mchol);
	223	return kTRUE;
	224	//
	225	}
de34b538	226
8a9ab0eb	227	//___________________________________________________________
	228	void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
	229	{
	230	// Invert matrix using Choleski decomposition, provided the Cholseki's L matrix
5d88242b	231	//
8a9ab0eb	232	Double_t sum;
	233	AliSymMatrix& mchol = *pmchol;
	234	//
	235	// Invert decomposed triangular L matrix (Lower triangle is filled)
5d88242b	236	for (int i=0;i<fNrowIndex;i++) {
8a9ab0eb	237	mchol(i,i) = 1.0/mchol(i,i);
5d88242b	238	for (int j=i+1;j<fNrowIndex;j++) {
de34b538	239	Double_t *rowj = mchol.GetRow(j);
8a9ab0eb	240	sum = 0.0;
de34b538	241	for (int k=i;k<j;k++) if (rowj[k]) {
	242	double &mki = mchol(k,i); if (mki) sum -= rowj[k]*mki;
	243	}
	244	rowj[i] = sum/rowj[j];
8a9ab0eb	245	}
	246	}
	247	//
	248	// take product of the inverted Choleski L matrix with its transposed
	249	for (int i=fNrowIndex;i--;) {
	250	for (int j=i+1;j--;) {
	251	sum = 0;
de34b538	252	for (int k=i;k<fNrowIndex;k++) {
	253	double &mik = mchol(i,k);
	254	if (mik) {
	255	double &mjk = mchol(j,k);
	256	if (mjk) sum += mik*mjk;
	257	}
	258	}
8a9ab0eb	259	(*this)(j,i) = sum;
	260	}
	261	}
	262	//
	263	}
	264
	265
	266	//___________________________________________________________
	267	Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert)
	268	{
de34b538	269	// Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
	270	// Solves the set of n linear equations A x = b,
	271	// where a is a positive-definite symmetric matrix.
	272	// a[1..n][1..n] is the output of the routine CholDecomposw.
	273	// Only the lower triangle of a is accessed. b[1..n] is input as the
	274	// right-hand side vector. The solution vector is returned in b[1..n].
8a9ab0eb	275	//
	276	Int_t i,k;
	277	Double_t sum;
	278	//
	279	AliSymMatrix *pmchol = DecomposeChol();
	280	if (!pmchol) {
	281	printf("SolveChol failed\n");
5d88242b	282	// Print("l");
8a9ab0eb	283	return kFALSE;
	284	}
	285	AliSymMatrix& mchol = *pmchol;
	286	//
	287	for (i=0;i<fNrowIndex;i++) {
de34b538	288	Double_t *rowi = mchol.GetRow(i);
	289	for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k];
	290	b[i]=sum/rowi[i];
8a9ab0eb	291	}
de34b538	292	//
8a9ab0eb	293	for (i=fNrowIndex-1;i>=0;i--) {
de34b538	294	for (sum=b[i],k=i+1;k<fNrowIndex;k++) if (b[k]) {
	295	double &mki=mchol(k,i); if (mki) sum -= mki*b[k];
	296	}
8a9ab0eb	297	b[i]=sum/mchol(i,i);
	298	}
	299	//
	300	if (invert) InvertChol(pmchol);
	301	return kTRUE;
	302	//
	303	}
	304
	305	//___________________________________________________________
	306	Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert)
	307	{
	308	return SolveChol((Double_t*)b.GetMatrixArray(),invert);
	309	}
	310
	311
	312	//___________________________________________________________
	313	Bool_t AliSymMatrix::SolveChol(Double_t brhs, Double_t bsol,Bool_t invert)
	314	{
	315	memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
	316	return SolveChol(bsol,invert);
	317	}
	318
	319	//___________________________________________________________
	320	Bool_t AliSymMatrix::SolveChol(TVectorD &brhs, TVectorD &bsol,Bool_t invert)
	321	{
	322	bsol = brhs;
	323	return SolveChol(bsol,invert);
	324	}
	325
	326	//___________________________________________________________
	327	void AliSymMatrix::AddRows(int nrows)
	328	{
	329	if (nrows<1) return;
	330	Double_t *pnew = new Double_t[nrows+fNrows];
	331	for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows
	332	for (int ir=0;ir<nrows;ir++) {
	333	int ncl = GetSize()+1;
	334	pnew[fNrows] = new Double_t[ncl];
	335	memset(pnew[fNrows],0,ncl*sizeof(Double_t));
	336	fNrows++;
	337	fNrowIndex++;
	338	}
	339	delete[] fElemsAdd;
	340	fElemsAdd = pnew;
	341	//
	342	}
	343
	344	//___________________________________________________________
	345	void AliSymMatrix::Reset()
	346	{
	347	// if additional rows exist, regularize it
	348	if (fElemsAdd) {
	349	delete[] fElems;
	350	for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
	351	delete[] fElemsAdd; fElemsAdd = 0;
	352	fNcols = fNrowIndex;
	353	fElems = new Double_t[fNcols*(fNcols+1)/2];
	354	fNrows = 0;
	355	}
	356	if (fElems) memset(fElems,0,fNcols(fNcols+1)/2sizeof(Double_t));
	357	//
	358	}
	359
de34b538	360	//___________________________________________________________
	361	/*
	362	void AliSymMatrix::AddToRow(Int_t r, Double_t valc,Int_t indc,Int_t n)
	363	{
	364	// for (int i=n;i--;) {
	365	// (*this)(indc[i],r) += valc[i];
	366	// }
	367	// return;
	368
	369	double *row;
	370	if (r>=fNrowIndex) {
	371	AddRows(r-fNrowIndex+1);
	372	row = &((fElemsAdd[r-fNcols])[0]);
	373	}
	374	else row = &fElems[GetIndex(r,0)];
	375	//
	376	int nadd = 0;
	377	for (int i=n;i--;) {
	378	if (indc[i]>r) continue;
	379	row[indc[i]] += valc[i];
	380	nadd++;
	381	}
	382	if (nadd == n) return;
	383	//
	384	// add to col>row
	385	for (int i=n;i--;) {
	386	if (indc[i]>r) (*this)(indc[i],r) += valc[i];
	387	}
	388	//
	389	}
	390	*/
	391
	392	//___________________________________________________________
	393	Double_t* AliSymMatrix::GetRow(Int_t r)
	394	{
	395	if (r>=fNrowIndex) {
	396	int nn = fNrowIndex;
	397	AddRows(r-fNrowIndex+1);
	398	printf("create %d of %d\n",r, nn);
	399	return &((fElemsAdd[r-fNcols])[0]);
	400	}
	401	else return &fElems[GetIndex(r,0)];
	402	}
	403
	404
8a9ab0eb	405	//___________________________________________________________
	406	int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
	407	{
	408	// Solution a la MP1: gaussian eliminations
	409	/// Obtain solution of a system of linear equations with symmetric matrix
	410	/// and the inverse (using 'singular-value friendly' GAUSS pivot)
	411	//
	412
	413	Int_t nRank = 0;
	414	int iPivot;
	415	double vPivot = 0.;
	416	double eps = 0.00000000000001;
	417	int nGlo = GetSize();
	418	bool *bUnUsed = new bool[nGlo];
	419	double rowMax,colMax=0;
	420	rowMax = new double[nGlo];
	421	//
	422	if (stabilize) {
	423	colMax = new double[nGlo];
	424	for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0;
	425	for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) {
de34b538	426	double vl = TMath::Abs(Query(i,j));
8a9ab0eb	427	if (vl==0) continue;
	428	if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i
	429	if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j
	430	if (i==j) continue;
	431	if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j
	432	if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i
	433	}
	434	//
	435	for (Int_t i=nGlo; i--;) {
	436	if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
	437	if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
	438	}
	439	//
	440	}
	441	//
	442	for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
	443	//
	444	if (!fgBuffer \|\| fgBuffer->GetSize()!=GetSize()) {
	445	delete fgBuffer;
	446	try {
	447	fgBuffer = new AliSymMatrix(*this);
	448	}
	449	catch(bad_alloc&) {
	450	printf("Failed to allocate memory for matrix inversion buffer\n");
	451	return 0;
	452	}
	453	}
	454	else (fgBuffer) = this;
	455	//
	456	if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning)
	457	for (int j=0;j<=i; j++) {
de34b538	458	double vl = Query(i,j);
8a9ab0eb	459	if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])vlTMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
	460	}
	461	for (int j=i+1;j<nGlo;j++) {
de34b538	462	double vl = Query(j,i);
8a9ab0eb	463	if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])vlTMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
	464	}
	465	}
	466	//
de34b538	467	for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values
8a9ab0eb	468	//
	469	for (Int_t i=0; i<nGlo; i++) {
	470	vPivot = 0.0;
	471	iPivot = -1;
	472	//
	473	for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
	474	double vl;
de34b538	475	if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) {
8a9ab0eb	476	vPivot = vl;
	477	iPivot = j;
	478	}
	479	}
	480	//
	481	if (iPivot >= 0) { // pivot found
	482	nRank++;
	483	bUnUsed[iPivot] = false; // This value is used
	484	vPivot = 1.0/vPivot;
	485	DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse
	486	//
	487	for (Int_t j=0; j<nGlo; j++) {
	488	for (Int_t jj=0; jj<nGlo; jj++) {
	489	if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
	490	double &r = j>=jj ? (this)(j,jj) : (fgBuffer)(jj,j);
de34b538	491	r -= vPivot* ( j>iPivot ? Query(j,iPivot) : fgBuffer->Query(iPivot,j) )
de34b538	492	* ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot));
8a9ab0eb	493	}
	494	}
	495	}
	496	//
	497	for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements
	498	(this)(j,iPivot) = vPivot;
	499	(fgBuffer)(iPivot,j) = vPivot;
	500	}
	501	//
	502	}
	503	else { // No more pivot value (clear those elements)
	504	for (Int_t j=0; j<nGlo; j++) {
	505	if (bUnUsed[j]) {
	506	vecB[j] = 0.0;
	507	for (Int_t k=0; k<nGlo; k++) {
	508	(*this)(j,k) = 0.;
	509	if (j!=k) (*fgBuffer)(j,k) = 0;
	510	}
	511	}
	512	}
	513	break; // No more pivots anyway, stop here
	514	}
	515	}
	516	//
	517	for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) {
	518	double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V
	519	if (i>=j) (this)(i,j) = vl;
	520	else (fgBuffer)(j,i) = vl;
	521	}
	522	//
	523	for (Int_t j=0; j<nGlo; j++) {
	524	rowMax[j] = 0.0;
	525	for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
	526	double vl;
de34b538	527	if (j>=jj) vl = (*this)(j,jj) = -Query(j,jj);
de34b538	528	else vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj);
8a9ab0eb	529	rowMax[j] += vl*vecB[jj];
	530	}
	531	}
	532
	533	for (Int_t j=0; j<nGlo; j++) {
	534	vecB[j] = rowMax[j]; // The final result
	535	}
	536	//
	537	delete [] bUnUsed;
	538	delete [] rowMax;
	539	if (stabilize) delete [] colMax;
	540
	541	return nRank;
	542	}
5d88242b	543
5d88242b	544