STEER/STEER/AliSymMatrix.cxx

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