Adding Sym.Band-Diagonal matrix class + new preconditioners + coding conventions
[u/mrichter/AliRoot.git] / 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 }