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