]>
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; | |
7c185459 | 162 | for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (!IsZero(GetEl(i,j))) nel++; |
7c3070ec | 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&) { | |
7c185459 | 209 | AliInfo("Failed to allocate memory for Choleski decompostions"); |
8a9ab0eb | 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 | |
7c185459 | 225 | AliInfo(Form("The matrix is not positive definite [%e]\n" |
226 | "Choleski decomposition is not possible",sum)); | |
227 | Print("l"); | |
8a9ab0eb | 228 | return 0; |
229 | } | |
de34b538 | 230 | rowi[i] = TMath::Sqrt(sum); |
8a9ab0eb | 231 | // |
de34b538 | 232 | } else rowj[i] = sum/rowi[i]; |
8a9ab0eb | 233 | } |
234 | } | |
235 | return fgBuffer; | |
236 | } | |
237 | ||
238 | //___________________________________________________________ | |
239 | Bool_t AliSymMatrix::InvertChol() | |
240 | { | |
241 | // Invert matrix using Choleski decomposition | |
242 | // | |
243 | AliSymMatrix* mchol = DecomposeChol(); | |
244 | if (!mchol) { | |
7c185459 | 245 | AliInfo("Failed to invert the matrix"); |
8a9ab0eb | 246 | return kFALSE; |
247 | } | |
248 | // | |
249 | InvertChol(mchol); | |
250 | return kTRUE; | |
251 | // | |
252 | } | |
de34b538 | 253 | |
8a9ab0eb | 254 | //___________________________________________________________ |
255 | void AliSymMatrix::InvertChol(AliSymMatrix* pmchol) | |
256 | { | |
257 | // Invert matrix using Choleski decomposition, provided the Cholseki's L matrix | |
5d88242b | 258 | // |
8a9ab0eb | 259 | Double_t sum; |
260 | AliSymMatrix& mchol = *pmchol; | |
261 | // | |
262 | // Invert decomposed triangular L matrix (Lower triangle is filled) | |
7c3070ec | 263 | for (int i=0;i<GetSizeUsed();i++) { |
8a9ab0eb | 264 | mchol(i,i) = 1.0/mchol(i,i); |
7c3070ec | 265 | for (int j=i+1;j<GetSizeUsed();j++) { |
de34b538 | 266 | Double_t *rowj = mchol.GetRow(j); |
8a9ab0eb | 267 | sum = 0.0; |
de34b538 | 268 | for (int k=i;k<j;k++) if (rowj[k]) { |
269 | double &mki = mchol(k,i); if (mki) sum -= rowj[k]*mki; | |
270 | } | |
271 | rowj[i] = sum/rowj[j]; | |
8a9ab0eb | 272 | } |
273 | } | |
274 | // | |
275 | // take product of the inverted Choleski L matrix with its transposed | |
7c3070ec | 276 | for (int i=GetSizeUsed();i--;) { |
8a9ab0eb | 277 | for (int j=i+1;j--;) { |
278 | sum = 0; | |
7c3070ec | 279 | for (int k=i;k<GetSizeUsed();k++) { |
de34b538 | 280 | double &mik = mchol(i,k); |
281 | if (mik) { | |
282 | double &mjk = mchol(j,k); | |
283 | if (mjk) sum += mik*mjk; | |
284 | } | |
285 | } | |
8a9ab0eb | 286 | (*this)(j,i) = sum; |
287 | } | |
288 | } | |
289 | // | |
290 | } | |
291 | ||
292 | ||
293 | //___________________________________________________________ | |
294 | Bool_t AliSymMatrix::SolveChol(Double_t *b, Bool_t invert) | |
295 | { | |
de34b538 | 296 | // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com |
297 | // Solves the set of n linear equations A x = b, | |
298 | // where a is a positive-definite symmetric matrix. | |
299 | // a[1..n][1..n] is the output of the routine CholDecomposw. | |
300 | // Only the lower triangle of a is accessed. b[1..n] is input as the | |
301 | // right-hand side vector. The solution vector is returned in b[1..n]. | |
8a9ab0eb | 302 | // |
303 | Int_t i,k; | |
304 | Double_t sum; | |
305 | // | |
306 | AliSymMatrix *pmchol = DecomposeChol(); | |
307 | if (!pmchol) { | |
7c185459 | 308 | AliInfo("SolveChol failed"); |
5d88242b | 309 | // Print("l"); |
8a9ab0eb | 310 | return kFALSE; |
311 | } | |
312 | AliSymMatrix& mchol = *pmchol; | |
313 | // | |
7c3070ec | 314 | for (i=0;i<GetSizeUsed();i++) { |
de34b538 | 315 | Double_t *rowi = mchol.GetRow(i); |
316 | for (sum=b[i],k=i-1;k>=0;k--) if (rowi[k]&&b[k]) sum -= rowi[k]*b[k]; | |
317 | b[i]=sum/rowi[i]; | |
8a9ab0eb | 318 | } |
de34b538 | 319 | // |
7c3070ec | 320 | for (i=GetSizeUsed()-1;i>=0;i--) { |
321 | for (sum=b[i],k=i+1;k<GetSizeUsed();k++) if (b[k]) { | |
de34b538 | 322 | double &mki=mchol(k,i); if (mki) sum -= mki*b[k]; |
323 | } | |
8a9ab0eb | 324 | b[i]=sum/mchol(i,i); |
325 | } | |
326 | // | |
327 | if (invert) InvertChol(pmchol); | |
328 | return kTRUE; | |
329 | // | |
330 | } | |
331 | ||
332 | //___________________________________________________________ | |
333 | Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert) | |
334 | { | |
335 | return SolveChol((Double_t*)b.GetMatrixArray(),invert); | |
336 | } | |
337 | ||
338 | ||
339 | //___________________________________________________________ | |
340 | Bool_t AliSymMatrix::SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert) | |
341 | { | |
7c3070ec | 342 | memcpy(bsol,brhs,GetSizeUsed()*sizeof(Double_t)); |
8a9ab0eb | 343 | return SolveChol(bsol,invert); |
344 | } | |
345 | ||
346 | //___________________________________________________________ | |
347 | Bool_t AliSymMatrix::SolveChol(TVectorD &brhs, TVectorD &bsol,Bool_t invert) | |
348 | { | |
349 | bsol = brhs; | |
350 | return SolveChol(bsol,invert); | |
351 | } | |
352 | ||
353 | //___________________________________________________________ | |
354 | void AliSymMatrix::AddRows(int nrows) | |
355 | { | |
356 | if (nrows<1) return; | |
357 | Double_t **pnew = new Double_t*[nrows+fNrows]; | |
358 | for (int ir=0;ir<fNrows;ir++) pnew[ir] = fElemsAdd[ir]; // copy old extra rows | |
359 | for (int ir=0;ir<nrows;ir++) { | |
360 | int ncl = GetSize()+1; | |
361 | pnew[fNrows] = new Double_t[ncl]; | |
362 | memset(pnew[fNrows],0,ncl*sizeof(Double_t)); | |
363 | fNrows++; | |
364 | fNrowIndex++; | |
7c3070ec | 365 | fRowLwb++; |
8a9ab0eb | 366 | } |
367 | delete[] fElemsAdd; | |
368 | fElemsAdd = pnew; | |
369 | // | |
370 | } | |
371 | ||
372 | //___________________________________________________________ | |
373 | void AliSymMatrix::Reset() | |
374 | { | |
375 | // if additional rows exist, regularize it | |
376 | if (fElemsAdd) { | |
377 | delete[] fElems; | |
378 | for (int i=0;i<fNrows;i++) delete[] fElemsAdd[i]; | |
379 | delete[] fElemsAdd; fElemsAdd = 0; | |
7c3070ec | 380 | fNcols = fRowLwb = fNrowIndex; |
381 | fElems = new Double_t[GetSize()*(GetSize()+1)/2]; | |
8a9ab0eb | 382 | fNrows = 0; |
383 | } | |
7c3070ec | 384 | if (fElems) memset(fElems,0,GetSize()*(GetSize()+1)/2*sizeof(Double_t)); |
8a9ab0eb | 385 | // |
386 | } | |
387 | ||
de34b538 | 388 | //___________________________________________________________ |
389 | /* | |
390 | void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n) | |
391 | { | |
392 | // for (int i=n;i--;) { | |
393 | // (*this)(indc[i],r) += valc[i]; | |
394 | // } | |
395 | // return; | |
396 | ||
397 | double *row; | |
398 | if (r>=fNrowIndex) { | |
399 | AddRows(r-fNrowIndex+1); | |
400 | row = &((fElemsAdd[r-fNcols])[0]); | |
401 | } | |
402 | else row = &fElems[GetIndex(r,0)]; | |
403 | // | |
404 | int nadd = 0; | |
405 | for (int i=n;i--;) { | |
406 | if (indc[i]>r) continue; | |
407 | row[indc[i]] += valc[i]; | |
408 | nadd++; | |
409 | } | |
410 | if (nadd == n) return; | |
411 | // | |
412 | // add to col>row | |
413 | for (int i=n;i--;) { | |
414 | if (indc[i]>r) (*this)(indc[i],r) += valc[i]; | |
415 | } | |
416 | // | |
417 | } | |
418 | */ | |
419 | ||
420 | //___________________________________________________________ | |
421 | Double_t* AliSymMatrix::GetRow(Int_t r) | |
422 | { | |
7c3070ec | 423 | if (r>=GetSize()) { |
424 | int nn = GetSize(); | |
425 | AddRows(r-GetSize()+1); | |
7c185459 | 426 | AliDebug(2,Form("create %d of %d\n",r, nn)); |
7c3070ec | 427 | return &((fElemsAdd[r-GetSizeBooked()])[0]); |
de34b538 | 428 | } |
429 | else return &fElems[GetIndex(r,0)]; | |
430 | } | |
431 | ||
432 | ||
8a9ab0eb | 433 | //___________________________________________________________ |
434 | int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize) | |
435 | { | |
436 | // Solution a la MP1: gaussian eliminations | |
437 | /// Obtain solution of a system of linear equations with symmetric matrix | |
438 | /// and the inverse (using 'singular-value friendly' GAUSS pivot) | |
439 | // | |
440 | ||
441 | Int_t nRank = 0; | |
442 | int iPivot; | |
443 | double vPivot = 0.; | |
7c3070ec | 444 | double eps = 1e-14; |
445 | int nGlo = GetSizeUsed(); | |
8a9ab0eb | 446 | bool *bUnUsed = new bool[nGlo]; |
447 | double *rowMax,*colMax=0; | |
448 | rowMax = new double[nGlo]; | |
449 | // | |
450 | if (stabilize) { | |
451 | colMax = new double[nGlo]; | |
452 | for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0; | |
453 | for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { | |
de34b538 | 454 | double vl = TMath::Abs(Query(i,j)); |
7c185459 | 455 | if (IsZero(vl)) continue; |
8a9ab0eb | 456 | if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i |
457 | if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j | |
458 | if (i==j) continue; | |
459 | if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j | |
460 | if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i | |
461 | } | |
462 | // | |
463 | for (Int_t i=nGlo; i--;) { | |
7c185459 | 464 | if (!IsZero(rowMax[i])) rowMax[i] = 1./rowMax[i]; // Max elemt of row i |
465 | if (!IsZero(colMax[i])) colMax[i] = 1./colMax[i]; // Max elemt of column i | |
8a9ab0eb | 466 | } |
467 | // | |
468 | } | |
469 | // | |
470 | for (Int_t i=nGlo; i--;) bUnUsed[i] = true; | |
471 | // | |
7c3070ec | 472 | if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) { |
8a9ab0eb | 473 | delete fgBuffer; |
474 | try { | |
475 | fgBuffer = new AliSymMatrix(*this); | |
476 | } | |
477 | catch(bad_alloc&) { | |
7c185459 | 478 | AliError("Failed to allocate memory for matrix inversion buffer"); |
8a9ab0eb | 479 | return 0; |
480 | } | |
481 | } | |
482 | else (*fgBuffer) = *this; | |
483 | // | |
484 | if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) | |
485 | for (int j=0;j<=i; j++) { | |
de34b538 | 486 | double vl = Query(i,j); |
7c185459 | 487 | if (!IsZero(vl)) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 488 | } |
489 | for (int j=i+1;j<nGlo;j++) { | |
de34b538 | 490 | double vl = Query(j,i); |
7c185459 | 491 | if (!IsZero(vl)) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 492 | } |
493 | } | |
494 | // | |
de34b538 | 495 | for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values |
8a9ab0eb | 496 | // |
497 | for (Int_t i=0; i<nGlo; i++) { | |
498 | vPivot = 0.0; | |
499 | iPivot = -1; | |
500 | // | |
501 | for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element | |
502 | double vl; | |
de34b538 | 503 | if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) { |
8a9ab0eb | 504 | vPivot = vl; |
505 | iPivot = j; | |
506 | } | |
507 | } | |
508 | // | |
509 | if (iPivot >= 0) { // pivot found | |
510 | nRank++; | |
511 | bUnUsed[iPivot] = false; // This value is used | |
512 | vPivot = 1.0/vPivot; | |
513 | DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse | |
514 | // | |
515 | for (Int_t j=0; j<nGlo; j++) { | |
516 | for (Int_t jj=0; jj<nGlo; jj++) { | |
517 | if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!) | |
518 | double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j); | |
de34b538 | 519 | r -= vPivot* ( j>iPivot ? Query(j,iPivot) : fgBuffer->Query(iPivot,j) ) |
520 | * ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot)); | |
8a9ab0eb | 521 | } |
522 | } | |
523 | } | |
524 | // | |
525 | for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements | |
526 | (*this)(j,iPivot) *= vPivot; | |
527 | (*fgBuffer)(iPivot,j) *= vPivot; | |
528 | } | |
529 | // | |
530 | } | |
531 | else { // No more pivot value (clear those elements) | |
532 | for (Int_t j=0; j<nGlo; j++) { | |
533 | if (bUnUsed[j]) { | |
534 | vecB[j] = 0.0; | |
535 | for (Int_t k=0; k<nGlo; k++) { | |
536 | (*this)(j,k) = 0.; | |
537 | if (j!=k) (*fgBuffer)(j,k) = 0; | |
538 | } | |
539 | } | |
540 | } | |
541 | break; // No more pivots anyway, stop here | |
542 | } | |
543 | } | |
544 | // | |
545 | for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) { | |
546 | double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V | |
547 | if (i>=j) (*this)(i,j) *= vl; | |
548 | else (*fgBuffer)(j,i) *= vl; | |
549 | } | |
550 | // | |
551 | for (Int_t j=0; j<nGlo; j++) { | |
552 | rowMax[j] = 0.0; | |
553 | for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements | |
554 | double vl; | |
de34b538 | 555 | if (j>=jj) vl = (*this)(j,jj) = -Query(j,jj); |
556 | else vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj); | |
8a9ab0eb | 557 | rowMax[j] += vl*vecB[jj]; |
558 | } | |
559 | } | |
560 | ||
561 | for (Int_t j=0; j<nGlo; j++) { | |
562 | vecB[j] = rowMax[j]; // The final result | |
563 | } | |
564 | // | |
565 | delete [] bUnUsed; | |
566 | delete [] rowMax; | |
567 | if (stabilize) delete [] colMax; | |
568 | ||
569 | return nRank; | |
570 | } | |
5d88242b | 571 | |
572 |