7 #include "AliSymMatrix.h"
12 ClassImp(AliSymMatrix)
15 AliSymMatrix* AliSymMatrix::fgBuffer = 0;
16 Int_t AliSymMatrix::fgCopyCnt = 0;
17 //___________________________________________________________
18 AliSymMatrix::AliSymMatrix()
19 : fElems(0),fElemsAdd(0)
25 //___________________________________________________________
26 AliSymMatrix::AliSymMatrix(Int_t size)
27 : AliMatrixSq(),fElems(0),fElemsAdd(0)
31 fNrowIndex = fNcols = size;
32 fElems = new Double_t[fNcols*(fNcols+1)/2];
39 //___________________________________________________________
40 AliSymMatrix::AliSymMatrix(const AliSymMatrix &src)
41 : AliMatrixSq(src),fElems(0),fElemsAdd(0)
43 fNrowIndex = fNcols = src.GetSize();
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));
66 //___________________________________________________________
67 AliSymMatrix::~AliSymMatrix()
70 if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;}
73 //___________________________________________________________
74 AliSymMatrix& AliSymMatrix::operator=(const AliSymMatrix& 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];
85 fNrowIndex = fNcols = src.GetSize();
87 fElems = new Double_t[fNcols*(fNcols+1)/2];
88 int nmainel = src.fNcols*(src.fNcols+1);
89 memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
90 if (src.fNrows) { // transfer extra rows to main matrix
91 Double_t *pnt = fElems + nmainel*sizeof(Double_t);
92 int ncl = src.fNcols + 1;
93 for (int ir=0;ir<src.fNrows;ir++) {
95 memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
96 pnt += ncl*sizeof(Double_t);
102 memcpy(fElems,src.fElems,fNcols*(fNcols+1)/2*sizeof(Double_t));
103 int ncl = fNcols + 1;
104 for (int ir=0;ir<fNrows;ir++) { // dynamic rows
106 memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
114 //___________________________________________________________
115 void AliSymMatrix::Clear(Option_t*)
117 if (fElems) {delete[] fElems; fElems = 0;}
120 for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
124 fNrowIndex = fNcols = fNrows = 0;
128 //___________________________________________________________
129 Float_t AliSymMatrix::GetDensity() const
131 // get fraction of non-zero elements
133 for (int i=GetSize();i--;) for (int j=i+1;j--;) if (GetEl(i,j)!=0) nel++;
134 return 2.*nel/( (GetSize()+1)*GetSize() );
137 //___________________________________________________________
138 void AliSymMatrix::Print(Option_t* option) const
140 printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
141 TString opt = option; opt.ToLower();
142 if (opt.IsNull()) return;
143 opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|";
144 for (Int_t i=0;i<fNrowIndex;i++) {
146 for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j));
151 //___________________________________________________________
152 void AliSymMatrix::MultiplyByVec(Double_t *vecIn,Double_t *vecOut) const
154 // fill vecOut by matrix*vecIn
155 // vector should be of the same size as the matrix
156 for (int i=fNrowIndex;i--;) {
158 for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
163 //___________________________________________________________
164 AliSymMatrix* AliSymMatrix::DecomposeChol()
166 // Return a matrix with Choleski decomposition
167 /*Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
168 consturcts Cholesky decomposition of SYMMETRIC and
169 POSITIVELY-DEFINED matrix a (a=L*Lt)
170 Only upper triangle of the matrix has to be filled.
171 In opposite to function from the book, the matrix is modified:
172 lower triangle and diagonal are refilled. */
174 if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
177 fgBuffer = new AliSymMatrix(*this);
180 printf("Failed to allocate memory for Choleski decompostions\n");
184 else (*fgBuffer) = *this;
186 AliSymMatrix& mchol = *fgBuffer;
188 for (int i=0;i<fNrowIndex;i++) {
189 for (int j=i;j<fNrowIndex;j++) {
190 double sum=(*this)(i,j);
191 for (int k=i-1;k>=0;k--) sum -= mchol(i,k)*mchol(j,k);
193 if (sum <= 0.0) { // not positive-definite
194 printf("The matrix is not positive definite [%e]\n"
195 "Choleski decomposition is not possible\n",sum);
198 mchol(i,i) = TMath::Sqrt(sum);
200 } else mchol(j,i) = sum/mchol(i,i);
206 //___________________________________________________________
207 Bool_t AliSymMatrix::InvertChol()
209 // Invert matrix using Choleski decomposition
211 AliSymMatrix* mchol = DecomposeChol();
213 printf("Failed to invert the matrix\n");
221 //___________________________________________________________
222 void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
224 // Invert matrix using Choleski decomposition, provided the Cholseki's L matrix
227 AliSymMatrix& mchol = *pmchol;
229 // Invert decomposed triangular L matrix (Lower triangle is filled)
230 for (i=0;i<fNrowIndex;i++) {
231 mchol(i,i) = 1.0/mchol(i,i);
232 for (j=i+1;j<fNrowIndex;j++) {
234 for (k=i;k<j;k++) sum -= mchol(j,k)*mchol(k,i);
235 mchol(j,i) = sum/mchol(j,j);
239 // take product of the inverted Choleski L matrix with its transposed
240 for (int i=fNrowIndex;i--;) {
241 for (int j=i+1;j--;) {
243 for (int k=i;k<fNrowIndex;k++) sum += mchol(i,k)*mchol(j,k);
251 //___________________________________________________________
252 Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert)
254 /* Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com
255 Solves the set of n linear equations A x = b,
256 where a is a positive-definite symmetric matrix.
257 a[1..n][1..n] is the output of the routine CholDecomposw.
258 Only the lower triangle of a is accessed. b[1..n] is input as the
259 right-hand side vector. The solution vector is returned in b[1..n].*/
264 AliSymMatrix *pmchol = DecomposeChol();
266 printf("SolveChol failed\n");
269 AliSymMatrix& mchol = *pmchol;
271 for (i=0;i<fNrowIndex;i++) {
272 for (sum=b[i],k=i-1;k>=0;k--) sum -= mchol(i,k)*b[k];
275 for (i=fNrowIndex-1;i>=0;i--) {
276 for (sum=b[i],k=i+1;k<fNrowIndex;k++) sum -= mchol(k,i)*b[k];
280 if (invert) InvertChol(pmchol);
285 //___________________________________________________________
286 Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert)
288 return SolveChol((Double_t*)b.GetMatrixArray(),invert);
292 //___________________________________________________________
293 Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert)
295 memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
296 return SolveChol(bsol,invert);
299 //___________________________________________________________
300 Bool_t AliSymMatrix::SolveChol(TVectorD &brhs, TVectorD &bsol,Bool_t invert)
303 return SolveChol(bsol,invert);
306 //___________________________________________________________
307 void AliSymMatrix::AddRows(int nrows)
310 Double_t **pnew = new Double_t*[nrows+fNrows];
311 for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows
312 for (int ir=0;ir<nrows;ir++) {
313 int ncl = GetSize()+1;
314 pnew[fNrows] = new Double_t[ncl];
315 memset(pnew[fNrows],0,ncl*sizeof(Double_t));
324 //___________________________________________________________
325 void AliSymMatrix::Reset()
327 // if additional rows exist, regularize it
330 for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i];
331 delete[] fElemsAdd; fElemsAdd = 0;
333 fElems = new Double_t[fNcols*(fNcols+1)/2];
336 if (fElems) memset(fElems,0,fNcols*(fNcols+1)/2*sizeof(Double_t));
340 //___________________________________________________________
341 int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
343 // Solution a la MP1: gaussian eliminations
344 /// Obtain solution of a system of linear equations with symmetric matrix
345 /// and the inverse (using 'singular-value friendly' GAUSS pivot)
351 double eps = 0.00000000000001;
352 int nGlo = GetSize();
353 bool *bUnUsed = new bool[nGlo];
354 double *rowMax,*colMax=0;
355 rowMax = new double[nGlo];
358 colMax = new double[nGlo];
359 for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0;
360 for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) {
361 double vl = TMath::Abs(Querry(i,j));
363 if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i
364 if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j
366 if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j
367 if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i
370 for (Int_t i=nGlo; i--;) {
371 if (0.0 != rowMax[i]) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
372 if (0.0 != colMax[i]) colMax[i] = 1./colMax[i]; // Max elemt of column i
377 for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
379 if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
382 fgBuffer = new AliSymMatrix(*this);
385 printf("Failed to allocate memory for matrix inversion buffer\n");
389 else (*fgBuffer) = *this;
391 if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning)
392 for (int j=0;j<=i; j++) {
393 double vl = Querry(i,j);
394 if (vl!=0) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
396 for (int j=i+1;j<nGlo;j++) {
397 double vl = Querry(j,i);
398 if (vl!=0) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
402 for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QuerryDiag(j)); // save diagonal elem absolute values
404 for (Int_t i=0; i<nGlo; i++) {
408 for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element
410 if (bUnUsed[j] && (TMath::Abs(vl=QuerryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QuerryDiag(j)))) {
416 if (iPivot >= 0) { // pivot found
418 bUnUsed[iPivot] = false; // This value is used
420 DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse
422 for (Int_t j=0; j<nGlo; j++) {
423 for (Int_t jj=0; jj<nGlo; jj++) {
424 if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!)
425 double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j);
426 r -= vPivot* ( j>iPivot ? Querry(j,iPivot) : fgBuffer->Querry(iPivot,j) )
427 * ( iPivot>jj ? Querry(iPivot,jj) : fgBuffer->Querry(jj,iPivot));
432 for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements
433 (*this)(j,iPivot) *= vPivot;
434 (*fgBuffer)(iPivot,j) *= vPivot;
438 else { // No more pivot value (clear those elements)
439 for (Int_t j=0; j<nGlo; j++) {
442 for (Int_t k=0; k<nGlo; k++) {
444 if (j!=k) (*fgBuffer)(j,k) = 0;
448 break; // No more pivots anyway, stop here
452 for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) {
453 double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V
454 if (i>=j) (*this)(i,j) *= vl;
455 else (*fgBuffer)(j,i) *= vl;
458 for (Int_t j=0; j<nGlo; j++) {
460 for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements
462 if (j>=jj) vl = (*this)(j,jj) = -Querry(j,jj);
463 else vl = (*fgBuffer)(j,jj) = -fgBuffer->Querry(j,jj);
464 rowMax[j] += vl*vecB[jj];
468 for (Int_t j=0; j<nGlo; j++) {
469 vecB[j] = rowMax[j]; // The final result
474 if (stabilize) delete [] colMax;