Fixes for coverity.
[u/mrichter/AliRoot.git] / 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 }