Corrected protection.
[u/mrichter/AliRoot.git] / STEER / AliSymBDMatrix.cxx
CommitLineData
7c3070ec 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
20using namespace std;
21
22ClassImp(AliSymBDMatrix)
23
24
25//___________________________________________________________
26AliSymBDMatrix::AliSymBDMatrix()
27: fElems(0)
28{
29 // def. c-tor
30 fSymmetric = kTRUE;
31}
32
33//___________________________________________________________
34AliSymBDMatrix::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//___________________________________________________________
55AliSymBDMatrix::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//___________________________________________________________
70AliSymBDMatrix::~AliSymBDMatrix()
71{
72 // d-tor
73 Clear();
74}
75
76//___________________________________________________________
77AliSymBDMatrix& 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//___________________________________________________________
101void 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//___________________________________________________________
111Float_t AliSymBDMatrix::GetDensity() const
112{
113 // get fraction of non-zero elements
114 if (!fNelems) return 0;
115 Int_t nel = 0;
7c185459 116 for (int i=fNelems;i--;) if (!IsZero(fElems[i])) nel++;
7c3070ec 117 return nel/fNelems;
118}
119
120//___________________________________________________________
121void 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//___________________________________________________________
551c9e69 137void AliSymBDMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const
7c3070ec 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//___________________________________________________________
168void AliSymBDMatrix::Reset()
169{
170 // set all elems to 0
171 if (fElems) memset(fElems,0,fNcols*fRowLwb*sizeof(Double_t));
c130d912 172 SetDecomposed(kFALSE);
7c3070ec 173 //
174}
175
176//___________________________________________________________
177void AliSymBDMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
178{
179 // add list of elements to row r
180 for (int i=0;i<n;i++) (*this)(r,indc[i]) = valc[i];
181}
182
183//___________________________________________________________
184void AliSymBDMatrix::DecomposeLDLT()
185{
186 // decomposition to L Diag L^T
187 if (IsDecomposed()) return;
188 //
189 Double_t eps = std::numeric_limits<double>::epsilon()*std::numeric_limits<double>::epsilon();
190 //
191 Double_t dtmp,gamma=0.0,xi=0.0;
192 int iDiag;
193 //
194 // find max diag and number of non-0 diag.elements
195 for (dtmp=0.0,iDiag=0;iDiag<GetSize();iDiag++) {
196 if ( (dtmp=QueryDiag(iDiag)) <=0.0) break;
197 if (gamma < dtmp) gamma = dtmp;
198 }
199 //
200 // find max. off-diag element
201 for (int ir=1;ir<iDiag;ir++) {
202 for (int ic=ir-GetBandHWidth();ic<ir;ic++) {
203 if (ic<0) continue;
204 dtmp = TMath::Abs(Query(ir,ic));
205 if (xi<dtmp) xi = dtmp;
206 }
207 }
208 double delta = eps*TMath::Max(1.0,xi+gamma);
209 //
210 double sn = GetSize()>1 ? 1.0/TMath::Sqrt( GetSize()*GetSize() - 1.0) : 1.0;
211 double beta = TMath::Sqrt(TMath::Max(eps,TMath::Max(gamma,xi*sn)));
212 //
213 for (int kr=1;kr<GetSize();kr++) {
214 int colKmin = TMath::Max(0,kr - GetBandHWidth());
215 double theta = 0.0;
216 //
217 for (int jr=colKmin;jr<=kr;jr++) {
218 int colJmin = TMath::Max(0,jr - GetBandHWidth());
219 //
220 dtmp = 0.0;
221 for (int i=TMath::Max(colKmin,colJmin);i<jr;i++)
222 dtmp += Query(kr,i)*QueryDiag(i)*Query(jr,i);
223 dtmp = (*this)(kr,jr) -= dtmp;
224 //
225 theta = TMath::Max(theta, TMath::Abs(dtmp));
226 //
227 if (jr!=kr) {
7c185459 228 if (!IsZero(QueryDiag(jr))) (*this)(kr,jr) /= QueryDiag(jr);
229 else (*this)(kr,jr) = 0.0;
7c3070ec 230 }
231 else if (kr<iDiag) {
232 dtmp = theta/beta;
233 dtmp *= dtmp;
234 dtmp = TMath::Max(dtmp, delta);
235 (*this)(kr,jr) = TMath::Max( TMath::Abs(Query(kr,jr)), dtmp);
236 }
237 } // jr
238 } // kr
239 //
240 for (int i=0;i<GetSize();i++) {
241 dtmp = QueryDiag(i);
7c185459 242 if (!IsZero(dtmp)) DiagElem(i) = 1./dtmp;
7c3070ec 243 }
244 //
245 SetDecomposed();
246}
247
248//___________________________________________________________
249void AliSymBDMatrix::Solve(Double_t *rhs)
250{
251 // solve matrix equation
252 //
253 if (!IsDecomposed()) DecomposeLDLT();
254 //
255 for (int kr=0;kr<GetSize();kr++)
256 for (int jr=TMath::Max(0,kr-GetBandHWidth());jr<kr;jr++) rhs[kr] -= Query(kr,jr)*rhs[jr];
257 //
258 for (int kr=GetSize();kr--;) rhs[kr] *= QueryDiag(kr);
259 //
260 for (int kr=GetSize();kr--;)
261 for (int jr=TMath::Max(0,kr - GetBandHWidth());jr<kr;jr++) rhs[jr] -= Query(kr,jr)*rhs[kr];
262 //
263}
264
265//___________________________________________________________
266void AliSymBDMatrix::Solve(const Double_t *rhs,Double_t *sol)
267{
268 // solve matrix equation
269 memcpy(sol,rhs,GetSize()*sizeof(Double_t));
270 Solve(sol);
271 //
272}