]>
Commit | Line | Data |
---|---|---|
7c3070ec | 1 | /**********************************************************************************************/ |
2 | /* Fast symmetric matrix with dynamically expandable size. */ | |
3 | /* Only part can be used for matrix operations. It is defined as: */ | |
4 | /* fNCols: rows built by constructor (GetSizeBooked) */ | |
5 | /* fNRows: number of rows added dynamically (automatically added on assignment to row) */ | |
6 | /* GetNRowAdded */ | |
7 | /* fNRowIndex: total size (fNCols+fNRows), GetSize */ | |
8 | /* fRowLwb : actual size to used for given operation, by default = total size, GetSizeUsed */ | |
9 | /* */ | |
10 | /* Author: ruben.shahoyan@cern.ch */ | |
11 | /* */ | |
12 | /**********************************************************************************************/ | |
8a9ab0eb | 13 | #include <stdlib.h> |
14 | #include <stdio.h> | |
15 | #include <iostream> | |
7c3070ec | 16 | #include <float.h> |
68f76645 | 17 | #include <string.h> |
8a9ab0eb | 18 | // |
7c3070ec | 19 | #include <TClass.h> |
20 | #include <TMath.h> | |
8a9ab0eb | 21 | #include "AliSymMatrix.h" |
7c3070ec | 22 | #include "AliLog.h" |
8a9ab0eb | 23 | // |
24 | ||
8a9ab0eb | 25 | ClassImp(AliSymMatrix) |
26 | ||
27 | ||
28 | AliSymMatrix* AliSymMatrix::fgBuffer = 0; | |
29 | Int_t AliSymMatrix::fgCopyCnt = 0; | |
30 | //___________________________________________________________ | |
31 | AliSymMatrix::AliSymMatrix() | |
32 | : fElems(0),fElemsAdd(0) | |
33 | { | |
68f76645 | 34 | // default constructor |
8a9ab0eb | 35 | fSymmetric = kTRUE; |
36 | fgCopyCnt++; | |
37 | } | |
38 | ||
39 | //___________________________________________________________ | |
40 | AliSymMatrix::AliSymMatrix(Int_t size) | |
41 | : AliMatrixSq(),fElems(0),fElemsAdd(0) | |
42 | { | |
68f76645 | 43 | //constructor for matrix with defined size |
8a9ab0eb | 44 | fNrows = 0; |
7c3070ec | 45 | fNrowIndex = fNcols = fRowLwb = size; |
8a9ab0eb | 46 | fElems = new Double_t[fNcols*(fNcols+1)/2]; |
47 | fSymmetric = kTRUE; | |
48 | Reset(); | |
49 | fgCopyCnt++; | |
50 | // | |
51 | } | |
52 | ||
53 | //___________________________________________________________ | |
54 | AliSymMatrix::AliSymMatrix(const AliSymMatrix &src) | |
55 | : AliMatrixSq(src),fElems(0),fElemsAdd(0) | |
56 | { | |
68f76645 | 57 | // copy constructor |
8a9ab0eb | 58 | fNrowIndex = fNcols = src.GetSize(); |
59 | fNrows = 0; | |
7c3070ec | 60 | fRowLwb = src.GetSizeUsed(); |
8a9ab0eb | 61 | if (fNcols) { |
62 | int nmainel = fNcols*(fNcols+1)/2; | |
63 | fElems = new Double_t[nmainel]; | |
64 | nmainel = src.fNcols*(src.fNcols+1)/2; | |
65 | memcpy(fElems,src.fElems,nmainel*sizeof(Double_t)); | |
7c3070ec | 66 | if (src.GetSizeAdded()) { // transfer extra rows to main matrix |
8a9ab0eb | 67 | Double_t *pnt = fElems + nmainel; |
7c3070ec | 68 | int ncl = src.GetSizeBooked() + 1; |
69 | for (int ir=0;ir<src.GetSizeAdded();ir++) { | |
8a9ab0eb | 70 | memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t)); |
71 | pnt += ncl; | |
72 | ncl++; | |
73 | } | |
74 | } | |
75 | } | |
76 | else fElems = 0; | |
77 | fElemsAdd = 0; | |
78 | fgCopyCnt++; | |
79 | // | |
80 | } | |
81 | ||
82 | //___________________________________________________________ | |
83 | AliSymMatrix::~AliSymMatrix() | |
84 | { | |
85 | Clear(); | |
86 | if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;} | |
87 | } | |
88 | ||
89 | //___________________________________________________________ | |
90 | AliSymMatrix& AliSymMatrix::operator=(const AliSymMatrix& src) | |
91 | { | |
68f76645 | 92 | // assignment operator |
8a9ab0eb | 93 | if (this != &src) { |
94 | TObject::operator=(src); | |
7c3070ec | 95 | if (GetSizeBooked()!=src.GetSizeBooked() && GetSizeAdded()!=src.GetSizeAdded()) { |
8a9ab0eb | 96 | // recreate the matrix |
97 | if (fElems) delete[] fElems; | |
7c3070ec | 98 | for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i]; |
8a9ab0eb | 99 | delete[] fElemsAdd; |
100 | // | |
de34b538 | 101 | fNrowIndex = src.GetSize(); |
102 | fNcols = src.GetSize(); | |
8a9ab0eb | 103 | fNrows = 0; |
7c3070ec | 104 | fRowLwb = src.GetSizeUsed(); |
105 | fElems = new Double_t[GetSize()*(GetSize()+1)/2]; | |
106 | int nmainel = src.GetSizeBooked()*(src.GetSizeBooked()+1); | |
8a9ab0eb | 107 | memcpy(fElems,src.fElems,nmainel*sizeof(Double_t)); |
7c3070ec | 108 | if (src.GetSizeAdded()) { // transfer extra rows to main matrix |
c8f37c50 | 109 | Double_t *pnt = fElems + nmainel;//*sizeof(Double_t); |
7c3070ec | 110 | int ncl = src.GetSizeBooked() + 1; |
111 | for (int ir=0;ir<src.GetSizeAdded();ir++) { | |
8a9ab0eb | 112 | ncl += ir; |
113 | memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t)); | |
c8f37c50 | 114 | pnt += ncl;//*sizeof(Double_t); |
8a9ab0eb | 115 | } |
116 | } | |
117 | // | |
118 | } | |
119 | else { | |
7c3070ec | 120 | memcpy(fElems,src.fElems,GetSizeBooked()*(GetSizeBooked()+1)/2*sizeof(Double_t)); |
121 | int ncl = GetSizeBooked() + 1; | |
122 | for (int ir=0;ir<GetSizeAdded();ir++) { // dynamic rows | |
8a9ab0eb | 123 | ncl += ir; |
124 | memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t)); | |
125 | } | |
126 | } | |
127 | } | |
128 | // | |
129 | return *this; | |
130 | } | |
131 | ||
7c3070ec | 132 | //___________________________________________________________ |
133 | AliSymMatrix& AliSymMatrix::operator+=(const AliSymMatrix& src) | |
134 | { | |
68f76645 | 135 | // add operator |
7c3070ec | 136 | if (GetSizeUsed() != src.GetSizeUsed()) { |
137 | AliError("Matrix sizes are different"); | |
138 | return *this; | |
139 | } | |
140 | for (int i=0;i<GetSizeUsed();i++) for (int j=i;j<GetSizeUsed();j++) (*this)(j,i) += src(j,i); | |
141 | return *this; | |
142 | } | |
143 | ||
8a9ab0eb | 144 | //___________________________________________________________ |
145 | void AliSymMatrix::Clear(Option_t*) | |
146 | { | |
68f76645 | 147 | // clear dynamic part |
8a9ab0eb | 148 | if (fElems) {delete[] fElems; fElems = 0;} |
149 | // | |
150 | if (fElemsAdd) { | |
7c3070ec | 151 | for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i]; |
8a9ab0eb | 152 | delete[] fElemsAdd; |
153 | fElemsAdd = 0; | |
154 | } | |
7c3070ec | 155 | fNrowIndex = fNcols = fNrows = fRowLwb = 0; |
8a9ab0eb | 156 | // |
157 | } | |
158 | ||
159 | //___________________________________________________________ | |
160 | Float_t AliSymMatrix::GetDensity() const | |
161 | { | |
162 | // get fraction of non-zero elements | |
163 | Int_t nel = 0; | |
7c185459 | 164 | for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (!IsZero(GetEl(i,j))) nel++; |
7c3070ec | 165 | return 2.*nel/( (GetSizeUsed()+1)*GetSizeUsed() ); |
8a9ab0eb | 166 | } |
167 | ||
168 | //___________________________________________________________ | |
169 | void AliSymMatrix::Print(Option_t* option) const | |
170 | { | |
68f76645 | 171 | // print itself |
7c3070ec | 172 | printf("Symmetric Matrix: Size = %d (%d rows added dynamically), %d used\n",GetSize(),GetSizeAdded(),GetSizeUsed()); |
8a9ab0eb | 173 | TString opt = option; opt.ToLower(); |
174 | if (opt.IsNull()) return; | |
175 | opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|"; | |
7c3070ec | 176 | for (Int_t i=0;i<GetSizeUsed();i++) { |
8a9ab0eb | 177 | printf(opt,i); |
178 | for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j)); | |
179 | printf("\n"); | |
180 | } | |
181 | } | |
182 | ||
183 | //___________________________________________________________ | |
551c9e69 | 184 | void AliSymMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const |
8a9ab0eb | 185 | { |
186 | // fill vecOut by matrix*vecIn | |
187 | // vector should be of the same size as the matrix | |
7c3070ec | 188 | for (int i=GetSizeUsed();i--;) { |
8a9ab0eb | 189 | vecOut[i] = 0.0; |
7c3070ec | 190 | for (int j=GetSizeUsed();j--;) vecOut[i] += vecIn[j]*GetEl(i,j); |
8a9ab0eb | 191 | } |
192 | // | |
193 | } | |
194 | ||
195 | //___________________________________________________________ | |
196 | AliSymMatrix* AliSymMatrix::DecomposeChol() | |
197 | { | |
198 | // Return a matrix with Choleski decomposition | |
de34b538 | 199 | // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com |
200 | // consturcts Cholesky decomposition of SYMMETRIC and | |
201 | // POSITIVELY-DEFINED matrix a (a=L*Lt) | |
202 | // Only upper triangle of the matrix has to be filled. | |
203 | // In opposite to function from the book, the matrix is modified: | |
204 | // lower triangle and diagonal are refilled. | |
8a9ab0eb | 205 | // |
7c3070ec | 206 | if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) { |
8a9ab0eb | 207 | delete fgBuffer; |
3c5b4cc8 | 208 | fgBuffer = new AliSymMatrix(*this); |
8a9ab0eb | 209 | } |
210 | else (*fgBuffer) = *this; | |
211 | // | |
212 | AliSymMatrix& mchol = *fgBuffer; | |
213 | // | |
7c3070ec | 214 | for (int i=0;i<GetSizeUsed();i++) { |
de34b538 | 215 | Double_t *rowi = mchol.GetRow(i); |
7c3070ec | 216 | for (int j=i;j<GetSizeUsed();j++) { |
de34b538 | 217 | Double_t *rowj = mchol.GetRow(j); |
218 | double sum = rowj[i]; | |
219 | for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k]; | |
8a9ab0eb | 220 | if (i == j) { |
221 | if (sum <= 0.0) { // not positive-definite | |
7c185459 | 222 | AliInfo(Form("The matrix is not positive definite [%e]\n" |
223 | "Choleski decomposition is not possible",sum)); | |
224 | Print("l"); | |
8a9ab0eb | 225 | return 0; |
226 | } | |
de34b538 | 227 | rowi[i] = TMath::Sqrt(sum); |
8a9ab0eb | 228 | // |
de34b538 | 229 | } else rowj[i] = sum/rowi[i]; |
8a9ab0eb | 230 | } |
231 | } | |
232 | return fgBuffer; | |
233 | } | |
234 | ||
235 | //___________________________________________________________ | |
236 | Bool_t AliSymMatrix::InvertChol() | |
237 | { | |
238 | // Invert matrix using Choleski decomposition | |
239 | // | |
240 | AliSymMatrix* mchol = DecomposeChol(); | |
241 | if (!mchol) { | |
7c185459 | 242 | AliInfo("Failed to invert the matrix"); |
8a9ab0eb | 243 | return kFALSE; |
244 | } | |
245 | // | |
246 | InvertChol(mchol); | |
247 | return kTRUE; | |
248 | // | |
249 | } | |
de34b538 | 250 | |
8a9ab0eb | 251 | //___________________________________________________________ |
252 | void AliSymMatrix::InvertChol(AliSymMatrix* pmchol) | |
253 | { | |
254 | // Invert matrix using Choleski decomposition, provided the Cholseki's L matrix | |
5d88242b | 255 | // |
8a9ab0eb | 256 | Double_t sum; |
257 | AliSymMatrix& mchol = *pmchol; | |
258 | // | |
259 | // Invert decomposed triangular L matrix (Lower triangle is filled) | |
7c3070ec | 260 | for (int i=0;i<GetSizeUsed();i++) { |
8a9ab0eb | 261 | mchol(i,i) = 1.0/mchol(i,i); |
7c3070ec | 262 | for (int j=i+1;j<GetSizeUsed();j++) { |
de34b538 | 263 | Double_t *rowj = mchol.GetRow(j); |
8a9ab0eb | 264 | sum = 0.0; |
de34b538 | 265 | for (int k=i;k<j;k++) if (rowj[k]) { |
266 | double &mki = mchol(k,i); if (mki) sum -= rowj[k]*mki; | |
267 | } | |
268 | rowj[i] = sum/rowj[j]; | |
8a9ab0eb | 269 | } |
270 | } | |
271 | // | |
272 | // take product of the inverted Choleski L matrix with its transposed | |
7c3070ec | 273 | for (int i=GetSizeUsed();i--;) { |
8a9ab0eb | 274 | for (int j=i+1;j--;) { |
275 | sum = 0; | |
7c3070ec | 276 | for (int k=i;k<GetSizeUsed();k++) { |
de34b538 | 277 | double &mik = mchol(i,k); |
278 | if (mik) { | |
279 | double &mjk = mchol(j,k); | |
280 | if (mjk) sum += mik*mjk; | |
281 | } | |
282 | } | |
8a9ab0eb | 283 | (*this)(j,i) = sum; |
284 | } | |
285 | } | |
286 | // | |
287 | } | |
288 | ||
289 | ||
290 | //___________________________________________________________ | |
291 | Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert) | |
292 | { | |
de34b538 | 293 | // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com |
294 | // Solves the set of n linear equations A x = b, | |
295 | // where a is a positive-definite symmetric matrix. | |
296 | // a[1..n][1..n] is the output of the routine CholDecomposw. | |
297 | // Only the lower triangle of a is accessed. b[1..n] is input as the | |
298 | // right-hand side vector. The solution vector is returned in b[1..n]. | |
8a9ab0eb | 299 | // |
300 | Int_t i,k; | |
301 | Double_t sum; | |
302 | // | |
303 | AliSymMatrix *pmchol = DecomposeChol(); | |
304 | if (!pmchol) { | |
7c185459 | 305 | AliInfo("SolveChol failed"); |
5d88242b | 306 | // Print("l"); |
8a9ab0eb | 307 | return kFALSE; |
308 | } | |
309 | AliSymMatrix& mchol = *pmchol; | |
310 | // | |
7c3070ec | 311 | for (i=0;i<GetSizeUsed();i++) { |
de34b538 | 312 | Double_t *rowi = mchol.GetRow(i); |
313 | for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k]; | |
314 | b[i]=sum/rowi[i]; | |
8a9ab0eb | 315 | } |
de34b538 | 316 | // |
7c3070ec | 317 | for (i=GetSizeUsed()-1;i>=0;i--) { |
318 | for (sum=b[i],k=i+1;k<GetSizeUsed();k++) if (b[k]) { | |
de34b538 | 319 | double &mki=mchol(k,i); if (mki) sum -= mki*b[k]; |
320 | } | |
8a9ab0eb | 321 | b[i]=sum/mchol(i,i); |
322 | } | |
323 | // | |
324 | if (invert) InvertChol(pmchol); | |
325 | return kTRUE; | |
326 | // | |
327 | } | |
328 | ||
329 | //___________________________________________________________ | |
330 | Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert) | |
331 | { | |
332 | return SolveChol((Double_t*)b.GetMatrixArray(),invert); | |
333 | } | |
334 | ||
335 | ||
336 | //___________________________________________________________ | |
337 | Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert) | |
338 | { | |
7c3070ec | 339 | memcpy(bsol,brhs,GetSizeUsed()*sizeof(Double_t)); |
8a9ab0eb | 340 | return SolveChol(bsol,invert); |
341 | } | |
342 | ||
343 | //___________________________________________________________ | |
339fbe23 | 344 | Bool_t AliSymMatrix::SolveChol(const TVectorD &brhs, TVectorD &bsol,Bool_t invert) |
8a9ab0eb | 345 | { |
346 | bsol = brhs; | |
347 | return SolveChol(bsol,invert); | |
348 | } | |
349 | ||
350 | //___________________________________________________________ | |
351 | void AliSymMatrix::AddRows(int nrows) | |
352 | { | |
68f76645 | 353 | // add empty rows |
8a9ab0eb | 354 | if (nrows<1) return; |
355 | Double_t **pnew = new Double_t*[nrows+fNrows]; | |
356 | for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows | |
357 | for (int ir=0;ir<nrows;ir++) { | |
358 | int ncl = GetSize()+1; | |
359 | pnew[fNrows] = new Double_t[ncl]; | |
360 | memset(pnew[fNrows],0,ncl*sizeof(Double_t)); | |
361 | fNrows++; | |
362 | fNrowIndex++; | |
7c3070ec | 363 | fRowLwb++; |
8a9ab0eb | 364 | } |
365 | delete[] fElemsAdd; | |
366 | fElemsAdd = pnew; | |
367 | // | |
368 | } | |
369 | ||
370 | //___________________________________________________________ | |
371 | void AliSymMatrix::Reset() | |
372 | { | |
373 | // if additional rows exist, regularize it | |
374 | if (fElemsAdd) { | |
375 | delete[] fElems; | |
376 | for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; | |
377 | delete[] fElemsAdd; fElemsAdd = 0; | |
7c3070ec | 378 | fNcols = fRowLwb = fNrowIndex; |
379 | fElems = new Double_t[GetSize()*(GetSize()+1)/2]; | |
8a9ab0eb | 380 | fNrows = 0; |
381 | } | |
7c3070ec | 382 | if (fElems) memset(fElems,0,GetSize()*(GetSize()+1)/2*sizeof(Double_t)); |
8a9ab0eb | 383 | // |
384 | } | |
385 | ||
de34b538 | 386 | //___________________________________________________________ |
387 | /* | |
388 | void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n) | |
389 | { | |
390 | // for (int i=n;i--;) { | |
391 | // (*this)(indc[i],r) += valc[i]; | |
392 | // } | |
393 | // return; | |
394 | ||
395 | double *row; | |
396 | if (r>=fNrowIndex) { | |
397 | AddRows(r-fNrowIndex+1); | |
398 | row = &((fElemsAdd[r-fNcols])[0]); | |
399 | } | |
400 | else row = &fElems[GetIndex(r,0)]; | |
401 | // | |
402 | int nadd = 0; | |
403 | for (int i=n;i--;) { | |
404 | if (indc[i]>r) continue; | |
405 | row[indc[i]] += valc[i]; | |
406 | nadd++; | |
407 | } | |
408 | if (nadd == n) return; | |
409 | // | |
410 | // add to col>row | |
411 | for (int i=n;i--;) { | |
412 | if (indc[i]>r) (*this)(indc[i],r) += valc[i]; | |
413 | } | |
414 | // | |
415 | } | |
416 | */ | |
417 | ||
418 | //___________________________________________________________ | |
419 | Double_t* AliSymMatrix::GetRow(Int_t r) | |
420 | { | |
68f76645 | 421 | // get pointer on the row |
7c3070ec | 422 | if (r>=GetSize()) { |
423 | int nn = GetSize(); | |
424 | AddRows(r-GetSize()+1); | |
7c185459 | 425 | AliDebug(2,Form("create %d of %d\n",r, nn)); |
7c3070ec | 426 | return &((fElemsAdd[r-GetSizeBooked()])[0]); |
de34b538 | 427 | } |
428 | else return &fElems[GetIndex(r,0)]; | |
429 | } | |
430 | ||
431 | ||
8a9ab0eb | 432 | //___________________________________________________________ |
433 | int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize) | |
434 | { | |
435 | // Solution a la MP1: gaussian eliminations | |
436 | /// Obtain solution of a system of linear equations with symmetric matrix | |
437 | /// and the inverse (using 'singular-value friendly' GAUSS pivot) | |
438 | // | |
439 | ||
440 | Int_t nRank = 0; | |
441 | int iPivot; | |
442 | double vPivot = 0.; | |
7c3070ec | 443 | double eps = 1e-14; |
444 | int nGlo = GetSizeUsed(); | |
8a9ab0eb | 445 | bool *bUnUsed = new bool[nGlo]; |
446 | double *rowMax,*colMax=0; | |
447 | rowMax = new double[nGlo]; | |
448 | // | |
449 | if (stabilize) { | |
450 | colMax = new double[nGlo]; | |
451 | for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0; | |
452 | for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { | |
de34b538 | 453 | double vl = TMath::Abs(Query(i,j)); |
7c185459 | 454 | if (IsZero(vl)) continue; |
8a9ab0eb | 455 | if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i |
456 | if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j | |
457 | if (i==j) continue; | |
458 | if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j | |
459 | if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i | |
460 | } | |
461 | // | |
462 | for (Int_t i=nGlo; i--;) { | |
7c185459 | 463 | if (!IsZero(rowMax[i])) rowMax[i] = 1./rowMax[i]; // Max elemt of row i |
464 | if (!IsZero(colMax[i])) colMax[i] = 1./colMax[i]; // Max elemt of column i | |
8a9ab0eb | 465 | } |
466 | // | |
467 | } | |
468 | // | |
469 | for (Int_t i=nGlo; i--;) bUnUsed[i] = true; | |
470 | // | |
7c3070ec | 471 | if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) { |
8a9ab0eb | 472 | delete fgBuffer; |
3c5b4cc8 | 473 | fgBuffer = new AliSymMatrix(*this); |
8a9ab0eb | 474 | } |
475 | else (*fgBuffer) = *this; | |
476 | // | |
477 | if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) | |
478 | for (int j=0;j<=i; j++) { | |
de34b538 | 479 | double vl = Query(i,j); |
7c185459 | 480 | if (!IsZero(vl)) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 481 | } |
482 | for (int j=i+1;j<nGlo;j++) { | |
de34b538 | 483 | double vl = Query(j,i); |
7c185459 | 484 | if (!IsZero(vl)) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 485 | } |
486 | } | |
487 | // | |
de34b538 | 488 | for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values |
8a9ab0eb | 489 | // |
490 | for (Int_t i=0; i<nGlo; i++) { | |
491 | vPivot = 0.0; | |
492 | iPivot = -1; | |
493 | // | |
494 | for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element | |
495 | double vl; | |
de34b538 | 496 | if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) { |
8a9ab0eb | 497 | vPivot = vl; |
498 | iPivot = j; | |
499 | } | |
500 | } | |
501 | // | |
502 | if (iPivot >= 0) { // pivot found | |
503 | nRank++; | |
504 | bUnUsed[iPivot] = false; // This value is used | |
505 | vPivot = 1.0/vPivot; | |
506 | DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse | |
507 | // | |
508 | for (Int_t j=0; j<nGlo; j++) { | |
509 | for (Int_t jj=0; jj<nGlo; jj++) { | |
510 | if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!) | |
511 | double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j); | |
de34b538 | 512 | r -= vPivot* ( j>iPivot ? Query(j,iPivot) : fgBuffer->Query(iPivot,j) ) |
513 | * ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot)); | |
8a9ab0eb | 514 | } |
515 | } | |
516 | } | |
517 | // | |
518 | for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements | |
519 | (*this)(j,iPivot) *= vPivot; | |
520 | (*fgBuffer)(iPivot,j) *= vPivot; | |
521 | } | |
522 | // | |
523 | } | |
524 | else { // No more pivot value (clear those elements) | |
525 | for (Int_t j=0; j<nGlo; j++) { | |
526 | if (bUnUsed[j]) { | |
527 | vecB[j] = 0.0; | |
528 | for (Int_t k=0; k<nGlo; k++) { | |
529 | (*this)(j,k) = 0.; | |
530 | if (j!=k) (*fgBuffer)(j,k) = 0; | |
531 | } | |
532 | } | |
533 | } | |
534 | break; // No more pivots anyway, stop here | |
535 | } | |
536 | } | |
537 | // | |
c8f37c50 | 538 | if (stabilize) for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) { |
539 | double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V | |
540 | if (i>=j) (*this)(i,j) *= vl; | |
541 | else (*fgBuffer)(j,i) *= vl; | |
542 | } | |
8a9ab0eb | 543 | // |
544 | for (Int_t j=0; j<nGlo; j++) { | |
545 | rowMax[j] = 0.0; | |
546 | for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements | |
547 | double vl; | |
de34b538 | 548 | if (j>=jj) vl = (*this)(j,jj) = -Query(j,jj); |
549 | else vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj); | |
8a9ab0eb | 550 | rowMax[j] += vl*vecB[jj]; |
551 | } | |
552 | } | |
553 | ||
554 | for (Int_t j=0; j<nGlo; j++) { | |
555 | vecB[j] = rowMax[j]; // The final result | |
556 | } | |
557 | // | |
558 | delete [] bUnUsed; | |
559 | delete [] rowMax; | |
560 | if (stabilize) delete [] colMax; | |
561 | ||
562 | return nRank; | |
563 | } | |
5d88242b | 564 | |
565 |