Changed definition of the "zero" matrix element from DBL_MIN to 1e-64
authorshahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Mar 2010 14:54:54 +0000 (14:54 +0000)
committershahoian <shahoian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Mar 2010 14:54:54 +0000 (14:54 +0000)
STEER/AliMatrixSq.h
STEER/AliMinResSolve.cxx
STEER/AliSymBDMatrix.cxx
STEER/AliSymMatrix.cxx

index 1689439..e1eaff3 100644 (file)
@@ -56,6 +56,8 @@ class AliMatrixSq : public TMatrixDBase {
   virtual void Allocate      (Int_t ,Int_t ,Int_t , Int_t ,Int_t ,Int_t ) 
   {Error("Allocate","Dummy"); return;}
   //
+  static Bool_t       IsZero(Double_t x,Double_t thresh=1e-64) {return x>0 ? (x<thresh) : (x>-thresh);}
+  //
  protected:
   //
   void    Swap(int &r,int &c) const {int t=r;r=c;c=t;}
index ce7dc73..b8a2b8a 100644 (file)
@@ -703,7 +703,7 @@ 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 (TMath::Abs(rowM.GetElem(j))<DBL_MIN) continue;
+      if (AliMatrixSq::IsZero(rowM.GetElem(j))) 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) fDiagLU[i] = rowM.GetElem(j);
@@ -711,7 +711,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     }
     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;                   
+       if (AliMatrixSq::IsZero(vl)) continue;
        rowUi.GetElem(jw[col]) = vl;
       }
     //
@@ -737,7 +737,7 @@ Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
     jw[i] = -1;
     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
     //
-    if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
+    if( AliMatrixSq::IsZero(fDiagLU[i]) ) {
       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
       delete[] jw;
       return -1;
@@ -802,7 +802,7 @@ 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 (TMath::Abs(vl)<DBL_MIN) continue;
+      if (AliMatrixSq::IsZero(vl)) continue;
       if( j < i )   rowLi.GetElem(jw[j]) = vl;
       else if(j==i) fDiagLU[i] = vl;
       else rowUi.GetElem(jw[j]) = vl;
@@ -829,7 +829,7 @@ Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
     jw[i] = -1;
     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
     //
-    if( TMath::Abs(fDiagLU[i])<DBL_MIN ) {
+    if( AliMatrixSq::IsZero(fDiagLU[i])) {
       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
       delete[] jw;
       return -1;
@@ -889,7 +889,7 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
     // assign lof = 0 for matrix elements
     for(int j=0;j<row.GetNElems(); j++) {
       int col = row.GetIndex(j);
-      if (TMath::Abs(row.GetElem(j))<DBL_MIN) continue;  // !!!! matrix is sparse but sometimes 0 appears 
+      if (AliMatrixSq::IsZero(row.GetElem(j))) continue;  // !!!! matrix is sparse but sometimes 0 appears 
       if (col<i) {                      // L-part
        jbuf[incl] = col;
        levls[incl] = 0;
@@ -902,7 +902,7 @@ Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
       }
     }
     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)
+       if (AliMatrixSq::IsZero(matrix->Query(col,i))) continue;    // Due to the symmetry  == matrix(i,col)
        jbuf[incu] = col;
        levls[incu] = 0;
        iw[col] = incu++;
@@ -1021,7 +1021,7 @@ Int_t  AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
     //
     // assign lof = 0 for matrix elements
     for(int j=0;j<fSize; j++) {
-      if (TMath::Abs(fMatrix->Query(i,j))<DBL_MIN) continue;
+      if (AliMatrixSq::IsZero(fMatrix->Query(i,j))) continue;
       if (j<i) {                      // L-part
        jbuf[incl] = j;
        levls[incl] = 0;
index 38d3ae3..59bd02b 100644 (file)
@@ -113,7 +113,7 @@ 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++;
+  for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
   return nel/fNelems;
 }
 
@@ -225,8 +225,8 @@ void AliSymBDMatrix::DecomposeLDLT()
       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;
+       if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
+       else                        (*this)(kr,jr) = 0.0;
       }
       else if (kr<iDiag) {
        dtmp = theta/beta;
@@ -239,7 +239,7 @@ void AliSymBDMatrix::DecomposeLDLT()
   //
   for (int i=0;i<GetSize();i++) {
     dtmp = QueryDiag(i);
-    if (TMath::Abs(dtmp)>DBL_MIN) DiagElem(i) = 1./dtmp;
+    if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
   }
   //
   SetDecomposed();
index 184c511..d5a78ce 100644 (file)
@@ -159,7 +159,7 @@ Float_t AliSymMatrix::GetDensity() const
 {
   // get fraction of non-zero elements
   Int_t nel = 0;
-  for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (TMath::Abs(GetEl(i,j))>DBL_MIN) nel++;
+  for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (!IsZero(GetEl(i,j))) nel++;
   return 2.*nel/( (GetSizeUsed()+1)*GetSizeUsed() );
 }
 
@@ -206,7 +206,7 @@ AliSymMatrix* AliSymMatrix::DecomposeChol()
       fgBuffer = new AliSymMatrix(*this);
     }
     catch(bad_alloc&) {
-      printf("Failed to allocate memory for Choleski decompostions\n");
+      AliInfo("Failed to allocate memory for Choleski decompostions");
       return 0;
     }
   }
@@ -222,8 +222,9 @@ AliSymMatrix* AliSymMatrix::DecomposeChol()
       for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k];
       if (i == j) {
        if (sum <= 0.0) { // not positive-definite
-         printf("The matrix is not positive definite [%e]\n"
-                "Choleski decomposition is not possible\n",sum);
+         AliInfo(Form("The matrix is not positive definite [%e]\n"
+                      "Choleski decomposition is not possible",sum));
+         Print("l");
          return 0;
        }
        rowi[i] = TMath::Sqrt(sum);
@@ -241,7 +242,7 @@ Bool_t AliSymMatrix::InvertChol()
   //
   AliSymMatrix* mchol = DecomposeChol();
   if (!mchol) {
-    printf("Failed to invert the matrix\n");
+    AliInfo("Failed to invert the matrix");
     return kFALSE;
   }
   //
@@ -304,7 +305,7 @@ Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert)
   //
   AliSymMatrix *pmchol = DecomposeChol();
   if (!pmchol) {
-    printf("SolveChol failed\n");
+    AliInfo("SolveChol failed");
     //    Print("l");
     return kFALSE;
   }
@@ -422,7 +423,7 @@ Double_t* AliSymMatrix::GetRow(Int_t r)
   if (r>=GetSize()) {
     int nn = GetSize();
     AddRows(r-GetSize()+1); 
-    printf("create %d of %d\n",r, nn);
+    AliDebug(2,Form("create %d of %d\n",r, nn));
     return &((fElemsAdd[r-GetSizeBooked()])[0]);
   }
   else return &fElems[GetIndex(r,0)];
@@ -451,7 +452,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<DBL_MIN) continue;
+       if (IsZero(vl)) 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;
@@ -460,8 +461,8 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
       }
     //
     for (Int_t i=nGlo; 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
+      if (!IsZero(rowMax[i])) rowMax[i] = 1./rowMax[i]; // Max elemt of row i
+      if (!IsZero(colMax[i])) colMax[i] = 1./colMax[i]; // Max elemt of column i
     }
     //
   }
@@ -474,7 +475,7 @@ int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize)
       fgBuffer = new AliSymMatrix(*this);
     }
     catch(bad_alloc&) {
-      printf("Failed to allocate memory for matrix inversion buffer\n");
+      AliError("Failed to allocate memory for matrix inversion buffer");
       return 0;
     }
   }
@@ -483,11 +484,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 (TMath::Abs(vl)>DBL_MIN) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+       if (!IsZero(vl)) 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 (TMath::Abs(vl)>DBL_MIN) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
+       if (!IsZero(vl)) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix
       }
     }
   //