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