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