]>
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 | ||
25 | using namespace std; | |
26 | ||
27 | ClassImp(AliSymMatrix) | |
28 | ||
29 | ||
30 | AliSymMatrix* AliSymMatrix::fgBuffer = 0; | |
31 | Int_t AliSymMatrix::fgCopyCnt = 0; | |
32 | //___________________________________________________________ | |
33 | AliSymMatrix::AliSymMatrix() | |
34 | : fElems(0),fElemsAdd(0) | |
35 | { | |
68f76645 | 36 | // default constructor |
8a9ab0eb | 37 | fSymmetric = kTRUE; |
38 | fgCopyCnt++; | |
39 | } | |
40 | ||
41 | //___________________________________________________________ | |
42 | AliSymMatrix::AliSymMatrix(Int_t size) | |
43 | : AliMatrixSq(),fElems(0),fElemsAdd(0) | |
44 | { | |
68f76645 | 45 | //constructor for matrix with defined size |
8a9ab0eb | 46 | fNrows = 0; |
7c3070ec | 47 | fNrowIndex = fNcols = fRowLwb = size; |
8a9ab0eb | 48 | fElems = new Double_t[fNcols*(fNcols+1)/2]; |
49 | fSymmetric = kTRUE; | |
50 | Reset(); | |
51 | fgCopyCnt++; | |
52 | // | |
53 | } | |
54 | ||
55 | //___________________________________________________________ | |
56 | AliSymMatrix::AliSymMatrix(const AliSymMatrix &src) | |
57 | : AliMatrixSq(src),fElems(0),fElemsAdd(0) | |
58 | { | |
68f76645 | 59 | // copy constructor |
8a9ab0eb | 60 | fNrowIndex = fNcols = src.GetSize(); |
61 | fNrows = 0; | |
7c3070ec | 62 | fRowLwb = src.GetSizeUsed(); |
8a9ab0eb | 63 | if (fNcols) { |
64 | int nmainel = fNcols*(fNcols+1)/2; | |
65 | fElems = new Double_t[nmainel]; | |
66 | nmainel = src.fNcols*(src.fNcols+1)/2; | |
67 | memcpy(fElems,src.fElems,nmainel*sizeof(Double_t)); | |
7c3070ec | 68 | if (src.GetSizeAdded()) { // transfer extra rows to main matrix |
8a9ab0eb | 69 | Double_t *pnt = fElems + nmainel; |
7c3070ec | 70 | int ncl = src.GetSizeBooked() + 1; |
71 | for (int ir=0;ir<src.GetSizeAdded();ir++) { | |
8a9ab0eb | 72 | memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t)); |
73 | pnt += ncl; | |
74 | ncl++; | |
75 | } | |
76 | } | |
77 | } | |
78 | else fElems = 0; | |
79 | fElemsAdd = 0; | |
80 | fgCopyCnt++; | |
81 | // | |
82 | } | |
83 | ||
84 | //___________________________________________________________ | |
85 | AliSymMatrix::~AliSymMatrix() | |
86 | { | |
87 | Clear(); | |
88 | if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;} | |
89 | } | |
90 | ||
91 | //___________________________________________________________ | |
92 | AliSymMatrix& AliSymMatrix::operator=(const AliSymMatrix& src) | |
93 | { | |
68f76645 | 94 | // assignment operator |
8a9ab0eb | 95 | if (this != &src) { |
96 | TObject::operator=(src); | |
7c3070ec | 97 | if (GetSizeBooked()!=src.GetSizeBooked() && GetSizeAdded()!=src.GetSizeAdded()) { |
8a9ab0eb | 98 | // recreate the matrix |
99 | if (fElems) delete[] fElems; | |
7c3070ec | 100 | for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i]; |
8a9ab0eb | 101 | delete[] fElemsAdd; |
102 | // | |
de34b538 | 103 | fNrowIndex = src.GetSize(); |
104 | fNcols = src.GetSize(); | |
8a9ab0eb | 105 | fNrows = 0; |
7c3070ec | 106 | fRowLwb = src.GetSizeUsed(); |
107 | fElems = new Double_t[GetSize()*(GetSize()+1)/2]; | |
108 | int nmainel = src.GetSizeBooked()*(src.GetSizeBooked()+1); | |
8a9ab0eb | 109 | memcpy(fElems,src.fElems,nmainel*sizeof(Double_t)); |
7c3070ec | 110 | if (src.GetSizeAdded()) { // transfer extra rows to main matrix |
c8f37c50 | 111 | Double_t *pnt = fElems + nmainel;//*sizeof(Double_t); |
7c3070ec | 112 | int ncl = src.GetSizeBooked() + 1; |
113 | for (int ir=0;ir<src.GetSizeAdded();ir++) { | |
8a9ab0eb | 114 | ncl += ir; |
115 | memcpy(pnt,src.fElemsAdd[ir],ncl*sizeof(Double_t)); | |
c8f37c50 | 116 | pnt += ncl;//*sizeof(Double_t); |
8a9ab0eb | 117 | } |
118 | } | |
119 | // | |
120 | } | |
121 | else { | |
7c3070ec | 122 | memcpy(fElems,src.fElems,GetSizeBooked()*(GetSizeBooked()+1)/2*sizeof(Double_t)); |
123 | int ncl = GetSizeBooked() + 1; | |
124 | for (int ir=0;ir<GetSizeAdded();ir++) { // dynamic rows | |
8a9ab0eb | 125 | ncl += ir; |
126 | memcpy(fElemsAdd[ir],src.fElemsAdd[ir],ncl*sizeof(Double_t)); | |
127 | } | |
128 | } | |
129 | } | |
130 | // | |
131 | return *this; | |
132 | } | |
133 | ||
7c3070ec | 134 | //___________________________________________________________ |
135 | AliSymMatrix& AliSymMatrix::operator+=(const AliSymMatrix& src) | |
136 | { | |
68f76645 | 137 | // add operator |
7c3070ec | 138 | if (GetSizeUsed() != src.GetSizeUsed()) { |
139 | AliError("Matrix sizes are different"); | |
140 | return *this; | |
141 | } | |
142 | for (int i=0;i<GetSizeUsed();i++) for (int j=i;j<GetSizeUsed();j++) (*this)(j,i) += src(j,i); | |
143 | return *this; | |
144 | } | |
145 | ||
8a9ab0eb | 146 | //___________________________________________________________ |
147 | void AliSymMatrix::Clear(Option_t*) | |
148 | { | |
68f76645 | 149 | // clear dynamic part |
8a9ab0eb | 150 | if (fElems) {delete[] fElems; fElems = 0;} |
151 | // | |
152 | if (fElemsAdd) { | |
7c3070ec | 153 | for (int i=0;i<GetSizeAdded();i++) delete[] fElemsAdd[i]; |
8a9ab0eb | 154 | delete[] fElemsAdd; |
155 | fElemsAdd = 0; | |
156 | } | |
7c3070ec | 157 | fNrowIndex = fNcols = fNrows = fRowLwb = 0; |
8a9ab0eb | 158 | // |
159 | } | |
160 | ||
161 | //___________________________________________________________ | |
162 | Float_t AliSymMatrix::GetDensity() const | |
163 | { | |
164 | // get fraction of non-zero elements | |
165 | Int_t nel = 0; | |
7c185459 | 166 | for (int i=GetSizeUsed();i--;) for (int j=i+1;j--;) if (!IsZero(GetEl(i,j))) nel++; |
7c3070ec | 167 | return 2.*nel/( (GetSizeUsed()+1)*GetSizeUsed() ); |
8a9ab0eb | 168 | } |
169 | ||
170 | //___________________________________________________________ | |
171 | void AliSymMatrix::Print(Option_t* option) const | |
172 | { | |
68f76645 | 173 | // print itself |
7c3070ec | 174 | printf("Symmetric Matrix: Size = %d (%d rows added dynamically), %d used\n",GetSize(),GetSizeAdded(),GetSizeUsed()); |
8a9ab0eb | 175 | TString opt = option; opt.ToLower(); |
176 | if (opt.IsNull()) return; | |
177 | opt = "%"; opt += 1+int(TMath::Log10(double(GetSize()))); opt+="d|"; | |
7c3070ec | 178 | for (Int_t i=0;i<GetSizeUsed();i++) { |
8a9ab0eb | 179 | printf(opt,i); |
180 | for (Int_t j=0;j<=i;j++) printf("%+.3e|",GetEl(i,j)); | |
181 | printf("\n"); | |
182 | } | |
183 | } | |
184 | ||
185 | //___________________________________________________________ | |
551c9e69 | 186 | void AliSymMatrix::MultiplyByVec(const Double_t *vecIn,Double_t *vecOut) const |
8a9ab0eb | 187 | { |
188 | // fill vecOut by matrix*vecIn | |
189 | // vector should be of the same size as the matrix | |
7c3070ec | 190 | for (int i=GetSizeUsed();i--;) { |
8a9ab0eb | 191 | vecOut[i] = 0.0; |
7c3070ec | 192 | for (int j=GetSizeUsed();j--;) vecOut[i] += vecIn[j]*GetEl(i,j); |
8a9ab0eb | 193 | } |
194 | // | |
195 | } | |
196 | ||
197 | //___________________________________________________________ | |
198 | AliSymMatrix* AliSymMatrix::DecomposeChol() | |
199 | { | |
200 | // Return a matrix with Choleski decomposition | |
de34b538 | 201 | // Adopted from Numerical Recipes in C, ch.2-9, http://www.nr.com |
202 | // consturcts Cholesky decomposition of SYMMETRIC and | |
203 | // POSITIVELY-DEFINED matrix a (a=L*Lt) | |
204 | // Only upper triangle of the matrix has to be filled. | |
205 | // In opposite to function from the book, the matrix is modified: | |
206 | // lower triangle and diagonal are refilled. | |
8a9ab0eb | 207 | // |
7c3070ec | 208 | if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) { |
8a9ab0eb | 209 | delete fgBuffer; |
3c5b4cc8 | 210 | fgBuffer = new AliSymMatrix(*this); |
8a9ab0eb | 211 | } |
212 | else (*fgBuffer) = *this; | |
213 | // | |
214 | AliSymMatrix& mchol = *fgBuffer; | |
215 | // | |
7c3070ec | 216 | for (int i=0;i<GetSizeUsed();i++) { |
de34b538 | 217 | Double_t *rowi = mchol.GetRow(i); |
7c3070ec | 218 | for (int j=i;j<GetSizeUsed();j++) { |
de34b538 | 219 | Double_t *rowj = mchol.GetRow(j); |
220 | double sum = rowj[i]; | |
221 | for (int k=i-1;k>=0;k--) if (rowi[k]&&rowj[k]) sum -= rowi[k]*rowj[k]; | |
8a9ab0eb | 222 | if (i == j) { |
223 | if (sum <= 0.0) { // not positive-definite | |
7c185459 | 224 | AliInfo(Form("The matrix is not positive definite [%e]\n" |
225 | "Choleski decomposition is not possible",sum)); | |
226 | Print("l"); | |
8a9ab0eb | 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) { | |
7c185459 | 244 | AliInfo("Failed to invert the matrix"); |
8a9ab0eb | 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) { | |
7c185459 | 307 | AliInfo("SolveChol failed"); |
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 | //___________________________________________________________ | |
339fbe23 | 346 | Bool_t AliSymMatrix::SolveChol(const TVectorD &brhs, TVectorD &bsol,Bool_t invert) |
8a9ab0eb | 347 | { |
348 | bsol = brhs; | |
349 | return SolveChol(bsol,invert); | |
350 | } | |
351 | ||
352 | //___________________________________________________________ | |
353 | void AliSymMatrix::AddRows(int nrows) | |
354 | { | |
68f76645 | 355 | // add empty rows |
8a9ab0eb | 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 | { | |
68f76645 | 423 | // get pointer on the row |
7c3070ec | 424 | if (r>=GetSize()) { |
425 | int nn = GetSize(); | |
426 | AddRows(r-GetSize()+1); | |
7c185459 | 427 | AliDebug(2,Form("create %d of %d\n",r, nn)); |
7c3070ec | 428 | return &((fElemsAdd[r-GetSizeBooked()])[0]); |
de34b538 | 429 | } |
430 | else return &fElems[GetIndex(r,0)]; | |
431 | } | |
432 | ||
433 | ||
8a9ab0eb | 434 | //___________________________________________________________ |
435 | int AliSymMatrix::SolveSpmInv(double *vecB, Bool_t stabilize) | |
436 | { | |
437 | // Solution a la MP1: gaussian eliminations | |
438 | /// Obtain solution of a system of linear equations with symmetric matrix | |
439 | /// and the inverse (using 'singular-value friendly' GAUSS pivot) | |
440 | // | |
441 | ||
442 | Int_t nRank = 0; | |
443 | int iPivot; | |
444 | double vPivot = 0.; | |
7c3070ec | 445 | double eps = 1e-14; |
446 | int nGlo = GetSizeUsed(); | |
8a9ab0eb | 447 | bool *bUnUsed = new bool[nGlo]; |
448 | double *rowMax,*colMax=0; | |
449 | rowMax = new double[nGlo]; | |
450 | // | |
451 | if (stabilize) { | |
452 | colMax = new double[nGlo]; | |
453 | for (Int_t i=nGlo; i--;) rowMax[i] = colMax[i] = 0.0; | |
454 | for (Int_t i=nGlo; i--;) for (Int_t j=i+1;j--;) { | |
de34b538 | 455 | double vl = TMath::Abs(Query(i,j)); |
7c185459 | 456 | if (IsZero(vl)) continue; |
8a9ab0eb | 457 | if (vl > rowMax[i]) rowMax[i] = vl; // Max elemt of row i |
458 | if (vl > colMax[j]) colMax[j] = vl; // Max elemt of column j | |
459 | if (i==j) continue; | |
460 | if (vl > rowMax[j]) rowMax[j] = vl; // Max elemt of row j | |
461 | if (vl > colMax[i]) colMax[i] = vl; // Max elemt of column i | |
462 | } | |
463 | // | |
464 | for (Int_t i=nGlo; i--;) { | |
7c185459 | 465 | if (!IsZero(rowMax[i])) rowMax[i] = 1./rowMax[i]; // Max elemt of row i |
466 | if (!IsZero(colMax[i])) colMax[i] = 1./colMax[i]; // Max elemt of column i | |
8a9ab0eb | 467 | } |
468 | // | |
469 | } | |
470 | // | |
471 | for (Int_t i=nGlo; i--;) bUnUsed[i] = true; | |
472 | // | |
7c3070ec | 473 | if (!fgBuffer || fgBuffer->GetSizeUsed()!=GetSizeUsed()) { |
8a9ab0eb | 474 | delete fgBuffer; |
3c5b4cc8 | 475 | fgBuffer = new AliSymMatrix(*this); |
8a9ab0eb | 476 | } |
477 | else (*fgBuffer) = *this; | |
478 | // | |
479 | if (stabilize) for (int i=0;i<nGlo; i++) { // Small loop for matrix equilibration (gives a better conditioning) | |
480 | for (int j=0;j<=i; j++) { | |
de34b538 | 481 | double vl = Query(i,j); |
7c185459 | 482 | if (!IsZero(vl)) SetEl(i,j, TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 483 | } |
484 | for (int j=i+1;j<nGlo;j++) { | |
de34b538 | 485 | double vl = Query(j,i); |
7c185459 | 486 | if (!IsZero(vl)) fgBuffer->SetEl(j,i,TMath::Sqrt(rowMax[i])*vl*TMath::Sqrt(colMax[j]) ); // Equilibrate the V matrix |
8a9ab0eb | 487 | } |
488 | } | |
489 | // | |
de34b538 | 490 | for (Int_t j=nGlo; j--;) fgBuffer->DiagElem(j) = TMath::Abs(QueryDiag(j)); // save diagonal elem absolute values |
8a9ab0eb | 491 | // |
492 | for (Int_t i=0; i<nGlo; i++) { | |
493 | vPivot = 0.0; | |
494 | iPivot = -1; | |
495 | // | |
496 | for (Int_t j=0; j<nGlo; j++) { // First look for the pivot, ie max unused diagonal element | |
497 | double vl; | |
de34b538 | 498 | if (bUnUsed[j] && (TMath::Abs(vl=QueryDiag(j))>TMath::Max(TMath::Abs(vPivot),eps*fgBuffer->QueryDiag(j)))) { |
8a9ab0eb | 499 | vPivot = vl; |
500 | iPivot = j; | |
501 | } | |
502 | } | |
503 | // | |
504 | if (iPivot >= 0) { // pivot found | |
505 | nRank++; | |
506 | bUnUsed[iPivot] = false; // This value is used | |
507 | vPivot = 1.0/vPivot; | |
508 | DiagElem(iPivot) = -vPivot; // Replace pivot by its inverse | |
509 | // | |
510 | for (Int_t j=0; j<nGlo; j++) { | |
511 | for (Int_t jj=0; jj<nGlo; jj++) { | |
512 | if (j != iPivot && jj != iPivot) {// Other elements (!!! do them first as you use old matV[k][j]'s !!!) | |
513 | double &r = j>=jj ? (*this)(j,jj) : (*fgBuffer)(jj,j); | |
de34b538 | 514 | r -= vPivot* ( j>iPivot ? Query(j,iPivot) : fgBuffer->Query(iPivot,j) ) |
515 | * ( iPivot>jj ? Query(iPivot,jj) : fgBuffer->Query(jj,iPivot)); | |
8a9ab0eb | 516 | } |
517 | } | |
518 | } | |
519 | // | |
520 | for (Int_t j=0; j<nGlo; j++) if (j != iPivot) { // Pivot row or column elements | |
521 | (*this)(j,iPivot) *= vPivot; | |
522 | (*fgBuffer)(iPivot,j) *= vPivot; | |
523 | } | |
524 | // | |
525 | } | |
526 | else { // No more pivot value (clear those elements) | |
527 | for (Int_t j=0; j<nGlo; j++) { | |
528 | if (bUnUsed[j]) { | |
529 | vecB[j] = 0.0; | |
530 | for (Int_t k=0; k<nGlo; k++) { | |
531 | (*this)(j,k) = 0.; | |
532 | if (j!=k) (*fgBuffer)(j,k) = 0; | |
533 | } | |
534 | } | |
535 | } | |
536 | break; // No more pivots anyway, stop here | |
537 | } | |
538 | } | |
539 | // | |
c8f37c50 | 540 | if (stabilize) for (Int_t i=0; i<nGlo; i++) for (Int_t j=0; j<nGlo; j++) { |
541 | double vl = TMath::Sqrt(colMax[i])*TMath::Sqrt(rowMax[j]); // Correct matrix V | |
542 | if (i>=j) (*this)(i,j) *= vl; | |
543 | else (*fgBuffer)(j,i) *= vl; | |
544 | } | |
8a9ab0eb | 545 | // |
546 | for (Int_t j=0; j<nGlo; j++) { | |
547 | rowMax[j] = 0.0; | |
548 | for (Int_t jj=0; jj<nGlo; jj++) { // Reverse matrix elements | |
549 | double vl; | |
de34b538 | 550 | if (j>=jj) vl = (*this)(j,jj) = -Query(j,jj); |
551 | else vl = (*fgBuffer)(j,jj) = -fgBuffer->Query(j,jj); | |
8a9ab0eb | 552 | rowMax[j] += vl*vecB[jj]; |
553 | } | |
554 | } | |
555 | ||
556 | for (Int_t j=0; j<nGlo; j++) { | |
557 | vecB[j] = rowMax[j]; // The final result | |
558 | } | |
559 | // | |
560 | delete [] bUnUsed; | |
561 | delete [] rowMax; | |
562 | if (stabilize) delete [] colMax; | |
563 | ||
564 | return nRank; | |
565 | } | |
5d88242b | 566 | |
567 |