STEER/STEER/AliSymBDMatrix.cxx

   1
   2 /*********************************************************************************/
   3 /* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal)       */
   4 /* Only lower triangle is stored in the "profile" format                         */
   5 /*                                                                               */
   6 /*                                                                               */
   7 /* Author: ruben.shahoyan@cern.ch                                                */
   8 /*                                                                               */
   9 /*********************************************************************************/
  10 #include <stdlib.h>
  11 #include <stdio.h>
  12 #include <iostream>
  13 #include <float.h>
  14 //
  15 #include "TClass.h"
  16 #include "TMath.h"
  17 #include "AliSymBDMatrix.h"
  18 //
  19
  20 ClassImp(AliSymBDMatrix)
  21
  22
  23 //___________________________________________________________
  24 AliSymBDMatrix::AliSymBDMatrix()
  25 : fElems(0)
  26 {
  27   // def. c-tor
  28   fSymmetric = kTRUE;
  29 }
  30
  31 //___________________________________________________________
  32 AliSymBDMatrix::AliSymBDMatrix(Int_t size, Int_t w)
  33   : AliMatrixSq(),fElems(0)
  34 {
  35   // c-tor for given size
  36   //
  37   fNcols = size;    // number of rows
  38   if (w<0) w = 0;
  39   if (w>=size) w = size-1;
  40   fNrows = w;
  41   fRowLwb = w+1;
  42   fSymmetric = kTRUE;
  43   //
  44   // total number of stored elements
  45   fNelems = size*(w+1) - w*(w+1)/2;
  46   //
  47   fElems = new Double_t[fNcols*fRowLwb];
  48   memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
  49   //
  50 }
  51
  52 //___________________________________________________________
  53 AliSymBDMatrix::AliSymBDMatrix(const AliSymBDMatrix &src)
  54   : AliMatrixSq(src),fElems(0)
  55 {
  56   // copy c-tor
  57   if (src.GetSize()<1) return;
  58   fNcols = src.GetSize();
  59   fNrows = src.GetBandHWidth();
  60   fRowLwb = fNrows+1;
  61   fNelems = src.GetNElemsStored();
  62   fElems = new Double_t[fNcols*fRowLwb];
  63   memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
  64   //
  65 }
  66
  67 //___________________________________________________________
  68 AliSymBDMatrix::~AliSymBDMatrix()
  69 {
  70   // d-tor
  71   Clear();
  72 }
  73
  74 //___________________________________________________________
  75 AliSymBDMatrix&  AliSymBDMatrix::operator=(const AliSymBDMatrix& src)
  76 {
  77   // assignment
  78   //
  79   if (this != &src) {
  80     TObject::operator=(src);
  81     if (fNcols!=src.fNcols) {
  82       // recreate the matrix
  83       if (fElems) delete[] fElems;
  84       fNcols = src.GetSize();
  85       fNrows = src.GetBandHWidth();
  86       fNelems = src.GetNElemsStored();
  87       fRowLwb = src.fRowLwb;
  88       fElems = new Double_t[fNcols*fRowLwb];
  89     }
  90     memcpy(fElems,src.fElems,fNcols*fRowLwb*sizeof(Double_t));
  91     fSymmetric = kTRUE;
  92   }
  93   //
  94   return *this;
  95   //
  96 }
  97
  98 //___________________________________________________________
  99 void AliSymBDMatrix::Clear(Option_t*)
 100 {
 101   // clear dynamic part
 102   if (fElems) {delete[] fElems; fElems = 0;}
 103   //
 104   fNelems = fNcols = fNrows = fRowLwb = 0;
 105   //
 106 }
 107
 108 //___________________________________________________________
 109 Float_t AliSymBDMatrix::GetDensity() const
 110 {
 111   // get fraction of non-zero elements
 112   if (!fNelems) return 0;
 113   Int_t nel = 0;
 114   for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
 115   return nel/fNelems;
 116 }
 117
 118 //___________________________________________________________
 119 void AliSymBDMatrix::Print(Option_t* option) const
 120 {
 121   // print data
 122   printf("Symmetric Band-Diagonal Matrix : Size = %d, half bandwidth = %d\n",
 123          GetSize(),GetBandHWidth());
 124   TString opt = option; opt.ToLower();
 125   if (opt.IsNull()) return;
 126   opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|";
 127   for (Int_t i=0;i<GetSize();i++) {
 128     printf(opt,i);
 129     for (Int_t j=TMath::Max(0,i-GetBandHWidth());j<=i;j++) printf("%+.3e|",GetEl(i,j));
 130     printf("\n");
 131   }
 132 }
 133
 134 //___________________________________________________________
 135 void AliSymBDMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const
 136 {
 137   // fill vecOut by matrix*vecIn
 138   // vector should be of the same size as the matrix
 139   if (IsDecomposed()) {
 140     for (int i=0;i<GetSize();i++) {
 141       double sm = 0;
 142       int jmax = TMath::Min(GetSize(),i+fRowLwb);
 143       for (int j=i+1;j<jmax;j++) sm += vecIn[j]*Query(j,i);
 144       vecOut[i] = QueryDiag(i)*(vecIn[i]+sm);
 145     }
 146     for (int i=GetSize();i--;) {
 147       double sm = 0;
 148       int jmin = TMath::Max(0,i - GetBandHWidth());
 149       int jmax = i-1;
 150       for (int j=jmin;j<jmax;j++) sm += vecOut[j]*Query(i,j);
 151       vecOut[i] += sm;
 152     }
 153   }
 154   else { // not decomposed
 155     for (int i=GetSize();i--;) {
 156       vecOut[i] = 0.0;
 157       int jmin = TMath::Max(0,i - GetBandHWidth());
 158       int jmax = TMath::Min(GetSize(),i + fRowLwb);
 159       for (int j=jmin;j<jmax;j++) vecOut[i] += vecIn[j]*Query(i,j);
 160     }
 161   }
 162   //
 163 }
 164
 165 //___________________________________________________________
 166 void AliSymBDMatrix::Reset()
 167 {
 168   // set all elems to 0
 169   if (fElems) memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
 170   SetDecomposed(kFALSE);
 171   //
 172 }
 173
 174 //___________________________________________________________
 175 void AliSymBDMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
 176 {
 177   // add list of elements to row r
 178   for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
 179 }
 180
 181 //___________________________________________________________
 182 void AliSymBDMatrix::DecomposeLDLT()
 183 {
 184   // decomposition to L Diag L^T
 185   if (IsDecomposed()) return;
 186   //
 187   Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
 188   //
 189   Double_t dtmp,gamma=0.0,xi=0.0;
 190   int iDiag;
 191   //
 192   // find max diag and number of non-0 diag.elements
 193   for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
 194     if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
 195     if (gamma < dtmp) gamma = dtmp;
 196   }
 197   //
 198   // find max. off-diag element
 199   for (int ir=1;ir<iDiag;ir++) {
 200     for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
 201       if (ic<0) continue;
 202       dtmp = TMath::Abs(Query(ir,ic));
 203       if (xi<dtmp) xi = dtmp;
 204     }
 205   }
 206   double delta = eps*TMath::Max(1.0,xi+gamma);
 207   //
 208   double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
 209   double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn)));
 210   //
 211   for (int kr=1;kr<GetSize();kr++) {
 212     int colKmin = TMath::Max(0,kr - GetBandHWidth());
 213     double theta = 0.0;
 214     //
 215     for (int jr=colKmin;jr<=kr;jr++) {
 216       int colJmin = TMath::Max(0,jr - GetBandHWidth());
 217       //
 218       dtmp = 0.0;
 219       for (int i=TMath::Max(colKmin,colJmin);i<jr;i++)
 220         dtmp += Query(kr,i)*QueryDiag(i)*Query(jr,i);
 221       dtmp = (*this)(kr,jr) -= dtmp;
 222       //
 223       theta = TMath::Max(theta, TMath::Abs(dtmp));
 224       //
 225       if (jr!=kr) {
 226         if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
 227         else                        (*this)(kr,jr) = 0.0;
 228       }
 229       else if (kr<iDiag) {
 230         dtmp = theta/beta;
 231         dtmp *= dtmp;
 232         dtmp = TMath::Max(dtmp, delta);
 233         (*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp);
 234       }
 235     } // jr
 236   } // kr
 237   //
 238   for (int i=0;i<GetSize();i++) {
 239     dtmp = QueryDiag(i);
 240     if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
 241   }
 242   //
 243   SetDecomposed();
 244 }
 245
 246 //___________________________________________________________
 247 void AliSymBDMatrix::Solve(Double_t *rhs)
 248 {
 249   // solve matrix equation
 250   //
 251   if (!IsDecomposed()) DecomposeLDLT();
 252   //
 253   for (int kr=0;kr<GetSize();kr++)
 254     for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
 255   //
 256   for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
 257   //
 258   for (int kr=GetSize();kr--;)
 259     for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
 260   //
 261 }
 262
 263 //___________________________________________________________
 264 void AliSymBDMatrix::Solve(const Double_t *rhs,Double_t *sol)
 265 {
 266   // solve matrix equation
 267   memcpy(sol,rhs,GetSize()*sizeof(Double_t));
 268   Solve(sol);
 269   //
 270 }