Adding Sym.Band-Diagonal matrix class + new preconditioners + coding conventions
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 18 Oct 2009 18:05:28 +0000 (18:05 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 18 Oct 2009 18:05:28 +0000 (18:05 +0000)
STEER/AliMatrixSq.cxx
STEER/AliMatrixSq.h
STEER/AliMillePede2.cxx
STEER/AliMillePede2.h
STEER/AliMinResSolve.cxx
STEER/AliMinResSolve.h
STEER/AliSymBDMatrix.cxx [new file with mode: 0644]
STEER/AliSymBDMatrix.h [new file with mode: 0644]
STEER/AliSymMatrix.cxx
STEER/AliSymMatrix.h

index f5e7232..f00439a 100644 (file)
@@ -1,10 +1,15 @@
+/**********************************************************************************************/
+/* Abstract class for matrix used for millepede2 operation.                                   */
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
 #include <stdlib.h>
 #include <stdio.h>
 #include <iostream>
 //
 #include "TClass.h"
 #include "TMath.h"
-#include "TVectorD.h"
 #include "AliMatrixSq.h"
 //
 
index 0cc7701..1689439 100644 (file)
@@ -1,5 +1,10 @@
 #ifndef ALIMATRIXSQ_H
 #define ALIMATRIXSQ_H
+/**********************************************************************************************/
+/* Abstract class for matrix used for millepede2 operation.                                   */
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
 
 #include <TMatrixDBase.h>
 #include <TVectorD.h>
index 62a3f94..67f0a45 100644 (file)
@@ -12,6 +12,7 @@
 #include <TStopwatch.h>
 #include <TFile.h>
 #include <TMath.h>
+#include <TVectorD.h>
 #include "AliMatrixSq.h"
 #include "AliSymMatrix.h"
 #include "AliRectMatrix.h"
@@ -449,6 +450,8 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   }
   double vl;
   //
+  Int_t maxLocUsed = 0;
+  //
   for (int ip=nPoints;ip--;) {  // Transfer the measurement records to matrices
     double  resid  = fRecord->GetValue( refLoc[ip]-1 );
     double  weight = fRecord->GetValue( refGlo[ip]-1 );    
@@ -467,11 +470,14 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     // Symmetric matrix, don't bother j>i coeffs
     for (int i=nrefLoc[ip];i--;) {         // Fill local matrix and vector
       fVecBLoc[ indLoc[i] ] += weight*resid*derLoc[i];
+      if (indLoc[i]>maxLocUsed) maxLocUsed = indLoc[i];  
       for (int j=i+1;j--;) matCLoc(indLoc[i] ,indLoc[j]) += weight*derLoc[i]*derLoc[j];
-    }   
+    }
     //
   } // end of the transfer of the measurement record to matrices
   //
+  matCLoc.SetSizeUsed(++maxLocUsed);   // data with B=0 may use less than declared nLocals 
+  //
   // first try to solve by faster Cholesky decomposition, then by Gaussian elimination
   if (!matCLoc.SolveChol(fVecBLoc,kTRUE)) {
     AliInfo("Failed to solve locals by Cholesky, traying Gaussian Elimination");
@@ -484,7 +490,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
   //
   // If requested, store the track params and errors
   //  printf("locfit: "); for (int i=0;i<fNLocPar;i++) printf("%+e |",fVecBLoc[i]); printf("\n");
-  if (localParams) for (int i=fNLocPar; i--;) {
+  if (localParams) for (int i=maxLocUsed; i--;) {
       localParams[2*i]   = fVecBLoc[i];
       localParams[2*i+1] = TMath::Sqrt(TMath::Abs(matCLoc.QueryDiag(i)));
     }
@@ -522,7 +528,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     nEq++;                        // number of equations                       
   } // end of Calculate residuals
   //
-  int nDoF = nEq-fNLocPar;
+  int nDoF = nEq-maxLocUsed;
   lChi2 = (nDoF>0) ? lChi2/nDoF : 0;  // Chi^2/dof  
   //
   if (fNStdDev != 0 && nDoF>0 && lChi2 > Chi2DoFLim(fNStdDev,nDoF)*fChi2CutFactor) { // check final chi2
@@ -581,7 +587,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
       int iCIDg = fGlo2CGlo[iIDg];  // compressed Index of index          
       if (iCIDg == -1) {
        Double_t *rowGL = matCGloLoc(nGloInFit);
-       for (int k=fNLocPar;k--;) rowGL[k] = 0.0;  // reset the row
+       for (int k=maxLocUsed;k--;) rowGL[k] = 0.0;  // reset the row
        iCIDg = fGlo2CGlo[iIDg] = nGloInFit;
        fCGlo2Glo[nGloInFit++] = iIDg;
       }
@@ -603,7 +609,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
     //
     vl = 0;
     Double_t *rowGLIDg =  matCGloLoc(iCIDg);
-    for (int kl=0;kl<fNLocPar;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
+    for (int kl=0;kl<maxLocUsed;kl++) if (rowGLIDg[kl]) vl += rowGLIDg[kl]*fVecBLoc[kl];
     if  (vl!=0) fVecBGlo[iIDg] -= fLocFitAdd ? vl : -vl;
     //
     int nfill = 0;
@@ -612,7 +618,7 @@ Int_t AliMillePede2::LocalFit(double *localParams)
       //
       vl = 0;
       Double_t *rowGLJDg =  matCGloLoc(jCIDg);
-      for (int kl=0;kl<fNLocPar;kl++) {
+      for (int kl=0;kl<maxLocUsed;kl++) {
        // diag terms
        if ( (vll=rowGLIDg[kl]*rowGLJDg[kl])!=0 ) vl += matCLoc.QueryDiag(kl)*vll;
        //
@@ -736,7 +742,7 @@ Int_t AliMillePede2::GlobalFitIteration()
   for (Long_t i=0;i<ndr;i++) {
     ReadRecordData(i);
     LocalFit();
-    if ( (i%int(0.1*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
+    if ( (i%int(0.2*ndr)) == 0) printf("%.1f%% of local fits done\n", double(100.*i)/ndr);
   }
   swt.Stop();
   printf("%ld local fits done: ", ndr);
@@ -792,6 +798,7 @@ Int_t AliMillePede2::GlobalFitIteration()
   // add large number to diagonal of fixed params  
   //
   for (int i=fNGloPar;i--;) { // // Reset row and column of fixed params and add 1/sig^2 to free ones
+    //    printf("#%3d : Nproc : %5d   grp: %d\n",i,fProcPnt[i],fParamGrID[i]);
     if (fProcPnt[i]<1) {
       fNGloFix++; 
       fVecBGlo[i] = 0.;
index 9d81293..1b107ae 100644 (file)
@@ -11,6 +11,7 @@
 /**********************************************************************************************/
 
 #include <TObject.h>
+#include <TString.h>
 #include <TTree.h>
 #include "AliMinResSolve.h"
 #include "AliMillePedeRecord.h"
index 54d47f0..740020e 100644 (file)
@@ -1,6 +1,19 @@
+/**********************************************************************************************/
+/* General class for solving large system of linear equations                                 */
+/* Includes MINRES, FGMRES methods as well as a few precondiotiong methods                    */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
 #include "AliMinResSolve.h"
 #include "AliLog.h"
+#include "AliMatrixSq.h"
+#include "AliMatrixSparse.h"
+#include "AliSymBDMatrix.h"
 #include <TStopwatch.h>
+#include <float.h>
+#include <TMath.h>
 
 using namespace std;
 
@@ -11,8 +24,10 @@ AliMinResSolve::AliMinResSolve() :
 fSize(0),fPrecon(0),fMatrix(0),fRHS(0),
   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
   fPvv(0),fPvz(0),fPhh(0),
-  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
-{}
+  fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+  // default constructor
+}
 
 //______________________________________________________
 AliMinResSolve::AliMinResSolve(const AliMinResSolve& src) : 
@@ -20,34 +35,42 @@ AliMinResSolve::AliMinResSolve(const AliMinResSolve& src) :
   fSize(src.fSize),fPrecon(src.fPrecon),fMatrix(src.fMatrix),fRHS(src.fRHS),
   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
   fPvv(0),fPvz(0),fPhh(0),
-  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
-{}
+  fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+  // copy constructor
+}
 
 //______________________________________________________
 AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) :
-  fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs->GetMatrixArray()),
+  fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs->GetMatrixArray()),
   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
   fPvv(0),fPvz(0),fPhh(0),
-  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
-{}
+  fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+  // copy accepting equation
+}
 
 //______________________________________________________
 AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) :
-  fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs),
+  fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs),
   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
   fPvv(0),fPvz(0),fPhh(0),
-  fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0)
-{}
+  fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
+{
+  // copy accepting equation
+}
 
 //______________________________________________________
 AliMinResSolve::~AliMinResSolve()
 {
+  // destructor
   ClearAux();
 }
 
 //______________________________________________________
 AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
 {
+  // assignment op.
   if (this != &src) {
     fSize    = src.fSize;
     fPrecon  = src.fPrecon;
@@ -60,27 +83,13 @@ AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
 //_______________________________________________________________
 Int_t AliMinResSolve::BuildPrecon(Int_t prec)
 {
-  const Double_t kTiny = 1E-12;
+  // preconditioner building
+  //  const Double_t kTiny = 1E-12;
   fPrecon = prec;
   //
-  if (fPrecon == kPreconDiag) return 0;     // nothing to do
-  //
-  else if (fPrecon == kPreconDILU) {
-    fPrecDiag = new Double_t[ fSize ];
-    fPrecAux  = new Double_t[ fSize ];
-    //
-    // copy inverse diagonal
-    for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QueryDiag(i);
-    //
-    for (int i=0;i<fSize;i++) {
-      fPrecDiag[i] = TMath::Abs(fPrecDiag[i])>kTiny  ? 1./fPrecDiag[i] : 1./TMath::Sign(kTiny,fPrecDiag[i]);
-      for (int j=i+1;j<fSize;j++) {
-       double vl = fMatrix->Query(j,i);
-       if (vl!=0) fPrecDiag[j] -= fPrecDiag[i]*vl*vl; 
-      }
-    }
-    return 0;
-  } // precon DILU
+  if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
+    return BuildPreconBD(fPrecon - kPreconBD);     // with halfbandwidth + diagonal = fPrecon
+  }
   //
   if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
     if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(fPrecon-kPreconILU0);
@@ -90,10 +99,10 @@ Int_t AliMinResSolve::BuildPrecon(Int_t prec)
   return -1;
 }
 
-
 //________________________________ FGMRES METHODS ________________________________
 Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
 {
+  // solve by fgmres
   return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov);
 }
 
@@ -156,7 +165,7 @@ Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double
     for (l=fSize;l--;) beta += fPvv[0][l]*fPvv[0][l];
     beta = TMath::Sqrt(beta);
     //
-    if (beta == 0.0) break;                                      // success? 
+    if (beta < epsmac) break;                                      // success? 
     t = 1.0 / beta;
     //--------------------   normalize:  fPvv[0] = fPvv[0] / beta 
     for (l=fSize;l--;) fPvv[0][l] *= t;
@@ -190,7 +199,7 @@ Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double
       // -------------------- h_{j+1,j} = ||w||_{2}
       for (t=0,l=fSize;l--;) t += fPvv[i1][l]*fPvv[i1][l]; t = TMath::Sqrt(t); 
       fPhh[i][i1] = t;
-      if (t != 0.0) for (t=1./t, l=0; l<fSize; l++) fPvv[i1][l] *= t;  //  v_{j+1} = w / h_{j+1,j}
+      if (t 0) for (t=1./t, l=0; l<fSize; l++) fPvv[i1][l] *= t;  //  v_{j+1} = w / h_{j+1,j}
       //
       // done with modified gram schimdt and arnoldi step
       // now  update factorization of fPhh
@@ -207,7 +216,7 @@ Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double
       //
       // if gamma is zero then any small value will do...
       // will affect only residual estimate
-      if (gam == 0.0) gam = epsmac;
+      if (gam < epsmac) gam = epsmac;
       //  get  next plane rotation
       fPVecR1[i] = fPhh[i][i]/gam;
       fPVecR2[i]  = fPhh[i][i1]/gam;
@@ -255,6 +264,7 @@ Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double
 //________________________________ MINRES METHODS ________________________________
 Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
 {
+  // solve by minres
   return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
 }
 
@@ -297,8 +307,8 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
   // ------------------------ initialization  ---------------------->>>>
   memset(VecSol,0,fSize*sizeof(double));
   int status=0, itn=0;
-  double Anorm = 0;
-  double Acond = 0;
+  double normA = 0;
+  double condA = 0;
   double ynorm = 0;
   double rnorm = 0;
   double gam,gmax=1,gmin=1,gbar,oldeps,epsa,epsx,epsr,diag, delta,phi,denom,z;
@@ -323,7 +333,7 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
     return kFALSE;
   }
   //
-  if (beta1 == 0) {
+  if (beta1 < eps) {
     AliInfo(Form("RHS is zero or is the nullspace of the Preconditioner: Solution is {0}"));
     ClearAux();
     return kTRUE;
@@ -446,11 +456,11 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
     rhs2   =       -  epsln * z;
     //
     //   Estimate various norms and test for convergence.
-    Anorm  = TMath::Sqrt( tnorm2 );
+    normA  = TMath::Sqrt( tnorm2 );
     ynorm  = TMath::Sqrt( ynorm2 );
-    epsa   = Anorm * eps;
-    epsx   = Anorm * ynorm * eps;
-    epsr   = Anorm * ynorm * rtol;
+    epsa   = normA * eps;
+    epsx   = normA * ynorm * eps;
+    epsr   = normA * ynorm * rtol;
     diag   = gbar;
     if (diag == 0) diag = epsa;
     //
@@ -463,14 +473,17 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
       where H is the tridiagonal matrix from Lanczos with one
       extra row, beta(k+1) e_k^T.
     */
-    Acond  = gmax / gmin;
+    condA  = gmax / gmin;
     //
     // See if any of the stopping criteria are satisfied.
     // In rare cases, istop is already -1 from above (Abar = const*I).
     //
+    AliInfo(Form("#%5d |qnrm: %+.2e Anrm:%+.2e Cnd:%+.2e Rnrm:%+.2e Ynrm:%+.2e EpsR:%+.2e EpsX:%+.2e Beta1:%+.2e",
+                itn,qrnorm, normA,condA,rnorm,ynorm,epsr,epsx,beta1));
+
     if (status == 0) {
       if (itn    >= itnlim    ) {status = 5; AliInfo(Form("%d iterations limit exceeded",itnlim));}
-      if (Acond  >= 0.1/eps   ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",Acond,0.1/eps));}
+      if (condA  >= 0.1/eps   ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",condA,0.1/eps));}
       if (epsx   >= beta1     ) {status = 3; AliInfo(Form("Approximate convergence"));}
       if (qrnorm <= epsx      ) {status = 2; AliInfo(Form("Converged within machine precision"));}
       if (qrnorm <= epsr      ) {status = 1; AliInfo(Form("Converged"));}
@@ -488,7 +501,7 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
               "Condition :  %+e\n"
               "Res.Norm  :  %+e\n"
               "Sol.Norm  :  %+e",
-              timer.CpuTime(),status,itn,Anorm,Acond,rnorm,ynorm));
+              timer.CpuTime(),status,itn,normA,condA,rnorm,ynorm));
   //
   return status>=0 && status<=3;
   //
@@ -497,6 +510,7 @@ Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double
 //______________________________________________________________
 void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
 {
+  // apply precond.
   ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
 }
 
@@ -505,26 +519,9 @@ void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
 {
   // Application of the preconditioner matrix:
   // implicitly defines the matrix solving the M*VecOut = VecRHS 
-  const Double_t kTiny = 1E-12;
-  if (fPrecon==kPreconDiag) {
-    for (int i=fSize;i--;) {
-      double d = fMatrix->QueryDiag(i);
-      vecOut[i] =  vecRHS[i]/(TMath::Abs(d)>kTiny ? d : kTiny);
-    }
-    //    return;
-  }
-  //
-  else if (fPrecon==kPreconDILU) {
-    for (int i=0;i<fSize;i++) {
-      double el = 0;
-      for (int j=i;j--;) {double vl = fMatrix->Query(i,j);if (vl!=0) el += fPrecAux[j]*vl;}
-      fPrecAux[i] = fPrecDiag[i]*(vecRHS[i]-el); 
-    }
-    for (int i=fSize;i--;) {
-      double el = 0;
-      for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Query(i,j);if (vl!=0) el += vl*vecOut[j];}
-      vecOut[i] = fPrecAux[i] - fPrecDiag[i]*el;
-    }
+  //  const Double_t kTiny = 1E-12;
+  if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
+    fMatBD->Solve(vecRHS,vecOut);
     //    return;
   }
   //
@@ -542,7 +539,7 @@ void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
       AliVectorSparse &rowUi = *fMatU->GetRow(i);
       int n = rowUi.GetNElems();
       for(int j=0;j<n;j++ ) vecOut[i] -= vecOut[ rowUi.GetIndex(j) ] * rowUi.GetElem(j);
-      vecOut[i] *= fPrecDiag[i];
+      vecOut[i] *= fDiagLU[i];
     }
     //
   }
@@ -553,6 +550,7 @@ void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
 //___________________________________________________________
 Bool_t AliMinResSolve::InitAuxMinRes()
 {
+  // init auxiliary space for minres
   try {
     fPVecY   = new double[fSize];   
     fPVecR1  = new double[fSize];   
@@ -577,6 +575,7 @@ Bool_t AliMinResSolve::InitAuxMinRes()
 //___________________________________________________________
 Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
 {
+  // init auxiliary space for fgmres
   try {
     fPvv     = new double*[nkrylov+1];
     fPvz     = new double*[nkrylov];
@@ -604,6 +603,7 @@ Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
 //___________________________________________________________
 void AliMinResSolve::ClearAux()
 {
+  // clear aux. space
   if (fPVecY)      delete[] fPVecY;  fPVecY   = 0;
   if (fPVecR1)     delete[] fPVecR1; fPVecR1  = 0; 
   if (fPVecR2)     delete[] fPVecR2; fPVecR2  = 0; 
@@ -611,10 +611,41 @@ void AliMinResSolve::ClearAux()
   if (fPVecW)      delete[] fPVecW;  fPVecW   = 0;
   if (fPVecW1)     delete[] fPVecW1; fPVecW1  = 0;
   if (fPVecW2)     delete[] fPVecW2; fPVecW2  = 0;
-  if (fPrecDiag)   delete[] fPrecDiag; fPrecDiag = 0;
-  if (fPrecAux)    delete[] fPrecAux; fPrecAux = 0;
+  if (fDiagLU)   delete[] fDiagLU; fDiagLU = 0;
   if (fMatL)       delete fMatL; fMatL = 0;
   if (fMatU)       delete fMatU; fMatU = 0;
+  if (fMatBD)      delete fMatBD; fMatBD = 0;
+}
+
+//___________________________________________________________
+Int_t  AliMinResSolve::BuildPreconBD(Int_t hwidth)
+{
+  // build preconditioner
+  AliInfo(Form("Building Band-Diagonal preconditioner of half-width = %d",hwidth));
+  fMatBD = new AliSymBDMatrix( fMatrix->GetSize(), hwidth );
+  //
+  // fill the band-diagonal part of the matrix
+  if (fMatrix->InheritsFrom("AliMatrixSparse")) {
+    for (int ir=fMatrix->GetSize();ir--;) {
+      int jmin = TMath::Max(0,ir-hwidth);
+      AliVectorSparse& irow = *((AliMatrixSparse*)fMatrix)->GetRow(ir);
+      for(int j=irow.GetNElems();j--;) {
+       int jind = irow.GetIndex(j);
+       if (jind<jmin) break;
+       (*fMatBD)(ir,jind) = irow.GetElem(j);
+      }
+    }
+  }
+  else {
+    for (int ir=fMatrix->GetSize();ir--;) {
+      int jmin = TMath::Max(0,ir-hwidth);
+      for(int jr=jmin;jr<=ir;jr++) (*fMatBD)(ir,jr) = fMatrix->Query(ir,jr);
+    }
+  }
+  //
+  fMatBD->DecomposeLDLT();
+  //
+  return 0;
 }
 
 //___________________________________________________________
@@ -627,14 +658,14 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
    *----------------------------------------------------------------------------*/
   //
   AliInfo(Form("Building ILU%d preconditioner",lofM));
-  AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
   //
   TStopwatch sw; sw.Start();
   fMatL = new AliMatrixSparse(fSize);
   fMatU = new AliMatrixSparse(fSize);
   fMatL->SetSymmetric(kFALSE);
   fMatU->SetSymmetric(kFALSE);
-  fPrecDiag = new Double_t[fSize];
+  fDiagLU = new Double_t[fSize];
+  AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
   //
   // symbolic factorization to calculate level of fill index arrays
   if ( PreconILUKsymb(lofM)<0 ) {
@@ -647,7 +678,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
   //
   for(int i=0; i<fSize; i++ ) {            // beginning of main loop
     if ( (i%int(0.1*fSize)) == 0) {
-      printf("BuildPrecon: row %d of %d\n",i,fSize);
+      AliInfo(Form("BuildPrecon: row %d of %d",i,fSize));
       sw.Stop();
       sw.Print();
       sw.Start(kFALSE);
@@ -655,7 +686,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     /* setup array jw[], and initial i-th row */
     AliVectorSparse& rowLi = *fMatL->GetRow(i);
     AliVectorSparse& rowUi = *fMatU->GetRow(i);
-    AliVectorSparse& rowM  = *Matrix->GetRow(i);
+    AliVectorSparse& rowM  = *matrix->GetRow(i);
     //
     for(int j=rowLi.GetNElems(); j--;) {  // initialize L part
       int col = rowLi.GetIndex(j);
@@ -663,7 +694,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
       rowLi.GetElem(j) = 0.;   // do we need this ?
     }
     jw[i] = i;
-    fPrecDiag[i] = 0;      // initialize diagonal
+    fDiagLU[i] = 0;      // initialize diagonal
     //
     for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
       int col = rowUi.GetIndex(j);
@@ -672,15 +703,15 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     }
     // copy row from csmat into L,U D
     for(int j=rowM.GetNElems(); j--;) {  // L and D part 
-      if (rowM.GetElem(j)==0) continue;
+      if (TMath::Abs(rowM.GetElem(j))<DBL_MIN) continue;
       int col = rowM.GetIndex(j);         // (the original matrix stores only lower triangle)
       if( col < i )   rowLi.GetElem(jw[col]) = rowM.GetElem(j); 
-      else if(col==i) fPrecDiag[i] = rowM.GetElem(j);
+      else if(col==i) fDiagLU[i] = rowM.GetElem(j);
       else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
     }
-    if (fMatrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) {      // part of the row I on the right of diagonal is stored as 
-       double vl = fMatrix->Query(col,i);    // the lower part of the column I
-       if (vl==0) continue;                   
+    if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) {      // part of the row I on the right of diagonal is stored as 
+       double vl = matrix->Query(col,i);    // the lower part of the column I
+       if (TMath::Abs(vl)<DBL_MIN) continue;                   
        rowUi.GetElem(jw[col]) = vl;
       }
     //
@@ -688,7 +719,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     for(int j=0; j<rowLi.GetNElems(); j++) {
       int jrow = rowLi.GetIndex(j);
       // get the multiplier for row to be eliminated (jrow)
-      rowLi.GetElem(j) *= fPrecDiag[jrow];
+      rowLi.GetElem(j) *= fDiagLU[jrow];
       //
       // combine current row and row jrow
       AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
@@ -697,7 +728,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
        int jpos = jw[col];
        if( jpos == -1 ) continue;
        if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
-       else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
        else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
       }
     }
@@ -706,19 +737,19 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     jw[i] = -1;
     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
     //
-    if( fPrecDiag[i] == 0 ) {
+    if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
       delete[] jw;
       return -1;
     }
-    fPrecDiag[i] = 1.0 / fPrecDiag[i];
+    fDiagLU[i] = 1.0 / fDiagLU[i];
   }
   //
   delete[] jw;
   //
   sw.Stop();
   AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
-  AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
+  AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
   //
   return 0;
 }
@@ -739,7 +770,7 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
   fMatU = new AliMatrixSparse(fSize);
   fMatL->SetSymmetric(kFALSE);
   fMatU->SetSymmetric(kFALSE);
-  fPrecDiag = new Double_t[fSize];
+  fDiagLU = new Double_t[fSize];
   //
   // symbolic factorization to calculate level of fill index arrays
   if ( PreconILUKsymbDense(lofM)<0 ) {
@@ -761,7 +792,7 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
       rowLi.GetElem(j) = 0.;   // do we need this ?
     }
     jw[i] = i;
-    fPrecDiag[i] = 0;      // initialize diagonal
+    fDiagLU[i] = 0;      // initialize diagonal
     //
     for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
       int col = rowUi.GetIndex(j);
@@ -771,16 +802,16 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
     // copy row from csmat into L,U D
     for(int j=fSize; j--;) {  // L and D part 
       double vl = fMatrix->Query(i,j);
-      if (vl==0) continue;
+      if (TMath::Abs(vl)<DBL_MIN) continue;
       if( j < i )   rowLi.GetElem(jw[j]) = vl;
-      else if(j==i) fPrecDiag[i] = vl;
+      else if(j==i) fDiagLU[i] = vl;
       else rowUi.GetElem(jw[j]) = vl;
     }
     // eliminate previous rows
     for(int j=0; j<rowLi.GetNElems(); j++) {
       int jrow = rowLi.GetIndex(j);
       // get the multiplier for row to be eliminated (jrow)
-      rowLi.GetElem(j) *= fPrecDiag[jrow];
+      rowLi.GetElem(j) *= fDiagLU[jrow];
       //
       // combine current row and row jrow
       AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
@@ -789,7 +820,7 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
        int jpos = jw[col];
        if( jpos == -1 ) continue;
        if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
-       else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
+       else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
        else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
       }
     }
@@ -798,12 +829,12 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
     jw[i] = -1;
     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
     //
-    if( fPrecDiag[i] == 0 ) {
+    if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
       delete[] jw;
       return -1;
     }
-    fPrecDiag[i] = 1.0 / fPrecDiag[i];
+    fDiagLU[i] = 1.0 / fDiagLU[i];
   }
   //
   delete[] jw;
@@ -825,9 +856,9 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
    *----------------------------------------------------------------------------*/
   //
   TStopwatch sw;
-  printf("PreconILUKsymb>>\n");
+  AliInfo("PreconILUKsymb>>");
+  AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix; 
   sw.Start();
-  AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix;
   //
   UChar_t **ulvl=0,*levls=0;
   UShort_t *jbuf=0;
@@ -853,12 +884,12 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
   for(int i=0; i<fSize; i++) {
     int incl = 0;
     int incu = i; 
-    AliVectorSparse& row = *Matrix->GetRow(i);
+    AliVectorSparse& row = *matrix->GetRow(i);
     //
     // assign lof = 0 for matrix elements
     for(int j=0;j<row.GetNElems(); j++) {
       int col = row.GetIndex(j);
-      if (row.GetElem(j)==0) continue;  // !!!! matrix is sparse but sometimes 0 appears 
+      if (TMath::Abs(row.GetElem(j))<DBL_MIN) continue;  // !!!! matrix is sparse but sometimes 0 appears 
       if (col<i) {                      // L-part
        jbuf[incl] = col;
        levls[incl] = 0;
@@ -870,8 +901,8 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
        iw[col] = incu++;
       }
     }
-    if (Matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
-       if (fMatrix->Query(col,i)==0) continue;    // Due to the symmetry  == matrix(i,col)
+    if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
+       if (TMath::Abs(matrix->Query(col,i))<DBL_MIN) continue;    // Due to the symmetry  == matrix(i,col)
        jbuf[incu] = col;
        levls[incu] = 0;
        iw[col] = incu++;
@@ -949,7 +980,7 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
   fMatU->SortIndices();
   sw.Stop();
   sw.Print();
-  printf("PreconILUKsymb<<\n");
+  AliInfo("PreconILUKsymb<<");
   return 0;
 }
 
@@ -990,7 +1021,7 @@ Int_t  AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
     //
     // assign lof = 0 for matrix elements
     for(int j=0;j<fSize; j++) {
-      if (fMatrix->Query(i,j)==0) continue;
+      if (TMath::Abs(fMatrix->Query(i,j))<DBL_MIN) continue;
       if (j<i) {                      // L-part
        jbuf[incl] = j;
        levls[incl] = 0;
index c9312e9..243ef3e 100644 (file)
@@ -1,18 +1,25 @@
 #ifndef ALIMINRESSOLVE_H
 #define ALIMINRESSOLVE_H
 
+/**********************************************************************************************/
+/* General class for solving large system of linear equations                                 */
+/* Includes MINRES, FGMRES methods as well as a few precondiotiong methods                    */
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
+
 #include <TObject.h>
 #include <TVectorD.h>
-#include <TMath.h>
-#include "AliMatrixSq.h"
-#include "AliMatrixSparse.h"
-class AliLog;
-class TStopwatch;
+class AliMatrixSq;
+class AliMatrixSparse;
+class AliSymBDMatrix;
+
 
 class AliMinResSolve : public TObject {
   //
  public:
-  enum {kPreconDiag=1,kPreconDILU=2,kPreconILU0=100,kPreconILU10=kPreconILU0+10,kPreconsTot};
+  enum {kPreconBD=1,kPreconILU0=100,kPreconILU10=kPreconILU0+10,kPreconsTot};
   enum {kSolMinRes,kSolFGMRes,kNSolvers};
  public:
   AliMinResSolve();
@@ -39,6 +46,7 @@ class AliMinResSolve : public TObject {
   Int_t  GetPrecon()                                               const    {return fPrecon;} 
   void   ClearAux();
   //
+  Int_t  BuildPreconBD(Int_t hwidth);
   Int_t  BuildPreconILUK(Int_t lofM);
   Int_t  BuildPreconILUKDense(Int_t lofM);
   Int_t  PreconILUKsymb(Int_t lofM);
@@ -48,13 +56,23 @@ class AliMinResSolve : public TObject {
   //
   Int_t               fSize;                             // dimension of the input matrix
   Int_t               fPrecon;                           // preconditioner type
-  const AliMatrixSq*  fMatrix;                           // matrix defining the equations
-  const Double_t*     fRHS;                              // right hand side
+  AliMatrixSq*        fMatrix;                           // matrix defining the equations
+  Double_t*           fRHS;                              // right hand side
   //
-  Double_t            *fPVecY,*fPVecR1,*fPVecR2,*fPVecV,*fPVecW,*fPVecW1,*fPVecW2;// aux MinRes
-  Double_t            **fPvv,**fPvz,**fPhh;
-  Double_t            *fPrecDiag,*fPrecAux; // aux space
-  AliMatrixSparse     *fMatL,*fMatU;
+  Double_t            *fPVecY;                           // aux. space
+  Double_t            *fPVecR1;                          // aux. space
+  Double_t            *fPVecR2;                          // aux. space
+  Double_t            *fPVecV;                           // aux. space
+  Double_t            *fPVecW;                           // aux. space
+  Double_t            *fPVecW1;                          // aux. space
+  Double_t            *fPVecW2;                          // aux. space
+  Double_t            **fPvv;                            // aux. space
+  Double_t            **fPvz;                            // aux. space
+  Double_t            **fPhh;                            // aux. space
+  Double_t            *fDiagLU;                          // aux space
+  AliMatrixSparse     *fMatL;                            // aux. space
+  AliMatrixSparse     *fMatU;                            // aux. space
+  AliSymBDMatrix      *fMatBD;                           // aux. space
   //
   ClassDef(AliMinResSolve,0)
 };
diff --git a/STEER/AliSymBDMatrix.cxx b/STEER/AliSymBDMatrix.cxx
new file mode 100644 (file)
index 0000000..ed7f785
--- /dev/null
@@ -0,0 +1,271 @@
+
+/*********************************************************************************/
+/* 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"
+//
+
+using namespace std;
+
+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 (TMath::Abs(fElems[i])>DBL_MIN) 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(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));
+  //
+}
+
+//___________________________________________________________
+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 (TMath::Abs(QueryDiag(jr))>DBL_MIN) (*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 (TMath::Abs(dtmp)>DBL_MIN) 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);
+  //
+}
diff --git a/STEER/AliSymBDMatrix.h b/STEER/AliSymBDMatrix.h
new file mode 100644 (file)
index 0000000..dc7f9a5
--- /dev/null
@@ -0,0 +1,125 @@
+#ifndef ALISYMBDMATRIX_H
+#define ALISYMBDMATRIX_H
+
+/*********************************************************************************/
+/* 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 <string.h>
+#include <TObject.h>
+#include <TVectorD.h>
+#include "AliMatrixSq.h"
+
+
+class AliSymBDMatrix : public AliMatrixSq {
+  //
+ public:
+  enum {kDecomposedBit = 0x1};
+  //
+  AliSymBDMatrix();
+  AliSymBDMatrix(Int_t size, Int_t w = 0);
+  AliSymBDMatrix(const AliSymBDMatrix &mat);
+  virtual ~AliSymBDMatrix();
+  //
+  Int_t         GetBandHWidth()                                  const {return fNrows;}
+  Int_t         GetNElemsStored()                                const {return fNelems;}
+  void          Clear(Option_t* option="");
+  void          Reset();
+  //
+  Float_t       GetDensity()                                     const;
+  AliSymBDMatrix& operator=(const AliSymBDMatrix& src);
+  Double_t      operator()(Int_t rown, Int_t coln)               const;
+  Double_t&     operator()(Int_t rown, Int_t coln);
+  Double_t      operator()(Int_t rown)                           const;
+  Double_t&     operator()(Int_t rown);
+  //
+  Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymBDMatrix*)this)(r,r);}
+  Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
+  void          DecomposeLDLT();
+  void          Solve(Double_t *rhs);
+  void          Solve(const Double_t *rhs,Double_t *sol);
+  void          Solve(TVectorD &rhs)                                   {Solve(rhs.GetMatrixArray());}
+  void          Solve(const TVectorD &rhs,TVectorD &sol)               {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());}
+  //
+  void          Print(Option_t* option="")                       const;
+  void          SetDecomposed()                                        {SetBit(kDecomposedBit);}
+  Bool_t        IsDecomposed()                                   const {return TestBit(kDecomposedBit);}
+  //
+  void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
+  void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
+  void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
+  //
+  // protected:
+  virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
+  virtual             Int_t GetIndex(Int_t diagID)             const;
+  Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
+  void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
+  //
+ protected:
+  Double_t*  fElems;     //   Elements booked by constructor
+  //
+  ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
+};
+
+
+//___________________________________________________________
+inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
+{
+  // lower triangle band is actually filled
+  if (row<col) Swap(row,col);
+  col -= row;
+  if (col < -GetBandHWidth()) return -1;
+  return GetIndex(row) + col;
+}
+
+//___________________________________________________________
+inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
+{
+  // Get index of the diagonal element on row diagID
+  return (diagID+1)*fRowLwb-1;
+}
+
+//___________________________________________________________
+inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
+{
+  // query element
+  int idx = GetIndex(row,col);
+  return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
+}
+
+//___________________________________________________________
+inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
+{
+  // get element for assingment; assignment outside of the stored range has no effect
+  int idx = GetIndex(row,col);  
+  if (idx>=0) return fElems[idx];
+  fTol = 0; 
+  return fTol;
+}
+
+//___________________________________________________________
+inline Double_t AliSymBDMatrix::operator()(Int_t row) const
+{
+  // query diagonal 
+  return (const Double_t&) fElems[GetIndex(row)];
+}
+
+//___________________________________________________________
+inline Double_t& AliSymBDMatrix::operator()(Int_t row)
+{
+  // get diagonal for assingment; assignment outside of the stored range has no effect
+  return fElems[GetIndex(row)];
+}
+
+//___________________________________________________________
+inline void AliSymBDMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
+{
+  MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
+}
+
+#endif
index 50bed50..184c511 100644 (file)
@@ -1,10 +1,24 @@
+/**********************************************************************************************/
+/* Fast symmetric matrix with dynamically expandable size.                                    */
+/* Only part can be used for matrix operations. It is defined as:                             */ 
+/* fNCols: rows built by constructor (GetSizeBooked)                                          */ 
+/* fNRows: number of rows added dynamically (automatically added on assignment to row)        */ 
+/*         GetNRowAdded                                                                       */ 
+/* fNRowIndex: total size (fNCols+fNRows), GetSize                                            */ 
+/* fRowLwb   : actual size to used for given operation, by default = total size, GetSizeUsed  */ 
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
 #include <stdlib.h>
 #include <stdio.h>
 #include <iostream>
+#include <float.h>
 //
-#include "TClass.h"
-#include "TMath.h"
+#include <TClass.h>
+#include <TMath.h>
 #include "AliSymMatrix.h"
+#include "AliLog.h"
 //
 
 using namespace std;
@@ -28,7 +42,7 @@ AliSymMatrix::AliSymMatrix(Int_t size)
 {
   //
   fNrows = 0;
-  fNrowIndex = fNcols = size;
+  fNrowIndex = fNcols = fRowLwb = size;
   fElems     = new Double_t[fNcols*(fNcols+1)/2];
   fSymmetric = kTRUE;
   Reset();
@@ -42,15 +56,16 @@ AliSymMatrix::AliSymMatrix(const AliSymMatrix &src)
 {
   fNrowIndex = fNcols = src.GetSize();
   fNrows = 0;
+  fRowLwb = src.GetSizeUsed();
   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
+    if (src.GetSizeAdded()) { // transfer extra rows to main matrix
       Double_t *pnt = fElems + nmainel;
-      int ncl = src.fNcols + 1;
-      for (int ir=0;ir<src.fNrows;ir++) {
+      int ncl = src.GetSizeBooked() + 1;
+      for (int ir=0;ir<src.GetSizeAdded();ir++) {
        memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
        pnt += ncl;
        ncl++; 
@@ -76,22 +91,23 @@ AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
   //
   if (this != &src) {
     TObject::operator=(src);
-    if (fNcols!=src.fNcols && fNrows!=src.fNrows) {
+    if (GetSizeBooked()!=src.GetSizeBooked() && GetSizeAdded()!=src.GetSizeAdded()) {
       // recreate the matrix
       if (fElems) delete[] fElems;
-      for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; 
+      for (int i=0;i<GetSizeAdded();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);
+      fRowLwb = src.GetSizeUsed();
+      fElems     = new Double_t[GetSize()*(GetSize()+1)/2];
+      int nmainel = src.GetSizeBooked()*(src.GetSizeBooked()+1);
       memcpy(fElems,src.fElems,nmainel*sizeof(Double_t));
-      if (src.fNrows) { // transfer extra rows to main matrix
+      if (src.GetSizeAdded()) { // 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++) {
+       int ncl = src.GetSizeBooked() + 1;
+       for (int ir=0;ir<src.GetSizeAdded();ir++) {
          ncl += ir; 
          memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t));
          pnt += ncl*sizeof(Double_t);
@@ -100,9 +116,9 @@ AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
       //
     }
     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
+      memcpy(fElems,src.fElems,GetSizeBooked()*(GetSizeBooked()+1)/2*sizeof(Double_t));
+      int ncl = GetSizeBooked() + 1;
+      for (int ir=0;ir<GetSizeAdded();ir++) { // dynamic rows
        ncl += ir; 
        memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t));
       }
@@ -112,19 +128,29 @@ AliSymMatrix&  AliSymMatrix::operator=(const AliSymMatrix& src)
   return *this;
 }
 
+//___________________________________________________________
+AliSymMatrix& AliSymMatrix::operator+=(const AliSymMatrix& src)
+{
+  //
+  if (GetSizeUsed() != src.GetSizeUsed()) {
+    AliError("Matrix sizes are different");
+    return *this;
+  }
+  for (int i=0;i<GetSizeUsed();i++) for (int j=i;j<GetSizeUsed();j++) (*this)(j,i) += src(j,i);
+  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]; 
+    for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i]; 
     delete[] fElemsAdd;
     fElemsAdd = 0;
   }
-  fNrowIndex = 0;
-  fNcols = 0;
-  fNrows = 0;
+  fNrowIndex = fNcols = fNrows = fRowLwb = 0;
   //
 }
 
@@ -133,18 +159,18 @@ 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() );
+  for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (TMath::Abs(GetEl(i,j))>DBL_MIN) nel++;
+  return 2.*nel/( (GetSizeUsed()+1)*GetSizeUsed() );
 }
 
 //___________________________________________________________
 void AliSymMatrix::Print(Option_t* option) const
 {
-  printf("Symmetric Matrix: Size = %d (%d rows added dynamically)\n",GetSize(),fNrows);
+  printf("Symmetric Matrix: Size = %d (%d rows added dynamically), %d used\n",GetSize(),GetSizeAdded(),GetSizeUsed());
   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++) {
+  for (Int_t i=0;i<GetSizeUsed();i++) {
     printf(opt,i);
     for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j));
     printf("\n");
@@ -156,9 +182,9 @@ 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--;) {
+  for (int i=GetSizeUsed();i--;) {
     vecOut[i] = 0.0;
-    for (int j=fNrowIndex;j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
+    for (int j=GetSizeUsed();j--;) vecOut[i] += vecIn[j]*GetEl(i,j);
   }
   //
 }
@@ -174,7 +200,7 @@ AliSymMatrix* AliSymMatrix::DecomposeChol()
   // In opposite to function from the book, the matrix is modified:
   // lower triangle and diagonal are refilled.
   //
-  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+  if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) {
     delete fgBuffer; 
     try {
       fgBuffer = new AliSymMatrix(*this);
@@ -188,9 +214,9 @@ AliSymMatrix* AliSymMatrix::DecomposeChol()
   //
   AliSymMatrix& mchol = *fgBuffer;
   //
-  for (int i=0;i<fNrowIndex;i++) {
+  for (int i=0;i<GetSizeUsed();i++) {
     Double_t *rowi = mchol.GetRow(i);
-    for (int j=i;j<fNrowIndex;j++) {
+    for (int j=i;j<GetSizeUsed();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];
@@ -233,9 +259,9 @@ void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
   AliSymMatrix& mchol = *pmchol;
   //
   // Invert decomposed triangular L matrix (Lower triangle is filled)
-  for (int i=0;i<fNrowIndex;i++) { 
+  for (int i=0;i<GetSizeUsed();i++) { 
     mchol(i,i) =  1.0/mchol(i,i);
-    for (int j=i+1;j<fNrowIndex;j++) { 
+    for (int j=i+1;j<GetSizeUsed();j++) { 
       Double_t *rowj = mchol.GetRow(j);
       sum = 0.0; 
       for (int k=i;k<j;k++) if (rowj[k]) { 
@@ -246,10 +272,10 @@ void AliSymMatrix::InvertChol(AliSymMatrix* pmchol)
   }
   //
   // take product of the inverted Choleski L matrix with its transposed
-  for (int i=fNrowIndex;i--;) {
+  for (int i=GetSizeUsed();i--;) {
     for (int j=i+1;j--;) {
       sum = 0;
-      for (int k=i;k<fNrowIndex;k++) {
+      for (int k=i;k<GetSizeUsed();k++) {
        double &mik = mchol(i,k); 
        if (mik) {
          double &mjk = mchol(j,k);
@@ -284,14 +310,14 @@ Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert)
   }
   AliSymMatrix& mchol = *pmchol;
   //
-  for (i=0;i<fNrowIndex;i++) {
+  for (i=0;i<GetSizeUsed();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]) {
+  for (i=GetSizeUsed()-1;i>=0;i--) {
+    for (sum=b[i],k=i+1;k<GetSizeUsed();k++) if (b[k]) {
       double &mki=mchol(k,i); if (mki) sum -= mki*b[k];
     }
     b[i]=sum/mchol(i,i);
@@ -312,7 +338,7 @@ Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert)
 //___________________________________________________________
 Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert) 
 {
-  memcpy(bsol,brhs,GetSize()*sizeof(Double_t));
+  memcpy(bsol,brhs,GetSizeUsed()*sizeof(Double_t));
   return SolveChol(bsol,invert);  
 }
 
@@ -335,6 +361,7 @@ void AliSymMatrix::AddRows(int nrows)
     memset(pnew[fNrows],0,ncl*sizeof(Double_t));
     fNrows++;
     fNrowIndex++;
+    fRowLwb++;
   }
   delete[] fElemsAdd;
   fElemsAdd = pnew;
@@ -349,11 +376,11 @@ void AliSymMatrix::Reset()
     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];
+    fNcols = fRowLwb = fNrowIndex;
+    fElems = new Double_t[GetSize()*(GetSize()+1)/2];
     fNrows = 0;
   }
-  if (fElems) memset(fElems,0,fNcols*(fNcols+1)/2*sizeof(Double_t));
+  if (fElems) memset(fElems,0,GetSize()*(GetSize()+1)/2*sizeof(Double_t));
   //
 }
 
@@ -392,11 +419,11 @@ void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
 //___________________________________________________________
 Double_t* AliSymMatrix::GetRow(Int_t r)
 {
-  if (r>=fNrowIndex) {
-    int nn = fNrowIndex;
-    AddRows(r-fNrowIndex+1); 
+  if (r>=GetSize()) {
+    int nn = GetSize();
+    AddRows(r-GetSize()+1); 
     printf("create %d of %d\n",r, nn);
-    return &((fElemsAdd[r-fNcols])[0]);
+    return &((fElemsAdd[r-GetSizeBooked()])[0]);
   }
   else return &fElems[GetIndex(r,0)];
 }
@@ -413,8 +440,8 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
   Int_t nRank = 0;
   int iPivot;
   double vPivot = 0.;
-  double eps = 0.00000000000001;
-  int nGlo = GetSize();
+  double eps = 1e-14;
+  int nGlo = GetSizeUsed();
   bool   *bUnUsed = new bool[nGlo];
   double *rowMax,*colMax=0;
   rowMax  = new double[nGlo];
@@ -424,7 +451,7 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
     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<DBL_MIN) 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;
@@ -433,15 +460,15 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
       }
     //
     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
+      if (TMath::Abs(rowMax[i])>DBL_MIN) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
+      if (TMath::Abs(colMax[i])>DBL_MIN) colMax[i] = 1./colMax[i]; // Max elemt of column i
     }
     //
   }
   //
   for (Int_t i=nGlo; i--;) bUnUsed[i] = true;
   //  
-  if (!fgBuffer || fgBuffer->GetSize()!=GetSize()) {
+  if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) {
     delete fgBuffer; 
     try {
       fgBuffer = new AliSymMatrix(*this);
@@ -456,11 +483,11 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
   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
+       if (TMath::Abs(vl)>DBL_MIN) 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
+       if (TMath::Abs(vl)>DBL_MIN) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
       }
     }
   //
index f22d5d5..1f0eaea 100644 (file)
@@ -1,5 +1,17 @@
 #ifndef ALISYMMATRIX_H
 #define ALISYMMATRIX_H
+/**********************************************************************************************/
+/* Fast symmetric matrix with dynamically expandable size.                                    */
+/* Only part can be used for matrix operations. It is defined as:                             */ 
+/* fNCols: rows built by constructor (GetSizeBooked)                                          */ 
+/* fNRows: number of rows added dynamically (automatically added on assignment to row)        */ 
+/*         GetNRowAdded                                                                       */ 
+/* fNRowIndex: total size (fNCols+fNRows), GetSize                                            */ 
+/* fRowLwb   : actual size to used for given operation, by default = total size, GetSizeUsed  */ 
+/*                                                                                            */ 
+/* Author: ruben.shahoyan@cern.ch                                                             */
+/*                                                                                            */ 
+/**********************************************************************************************/
 
 #include <string.h>
 #include <TObject.h>
@@ -20,8 +32,12 @@ class AliSymMatrix : public AliMatrixSq {
   void          Reset();
   //
   Int_t         GetSize()                                        const {return fNrowIndex;}
+  Int_t         GetSizeUsed()                                    const {return fRowLwb;}
+  Int_t         GetSizeBooked()                                  const {return fNcols;}
+  Int_t         GetSizeAdded()                                   const {return fNrows;}
   Float_t       GetDensity()                                     const;
   AliSymMatrix& operator=(const AliSymMatrix& src);
+  AliSymMatrix& operator+=(const AliSymMatrix& src);
   Double_t      operator()(Int_t rown, Int_t coln)               const;
   Double_t&     operator()(Int_t rown, Int_t coln);
   //
@@ -32,6 +48,7 @@ class AliSymMatrix : public AliMatrixSq {
   //
   void          Print(Option_t* option="")                       const;
   void          AddRows(int nrows=1);
+  void          SetSizeUsed(Int_t sz)                                  {fRowLwb = sz;}
   //
   void          Scale(Double_t coeff);
   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;