Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEER / AliSymMatrix.cxx
CommitLineData
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 25ClassImp(AliSymMatrix)
26
27
28AliSymMatrix* AliSymMatrix::fgBuffer = 0;
29Int_t AliSymMatrix::fgCopyCnt = 0;
30//___________________________________________________________
31AliSymMatrix::AliSymMatrix()
32: fElems(0),fElemsAdd(0)
33{
68f76645 34 // default constructor
8a9ab0eb 35 fSymmetric = kTRUE;
36 fgCopyCnt++;
37}
38
39//___________________________________________________________
40AliSymMatrix::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//___________________________________________________________
54AliSymMatrix::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//___________________________________________________________
83AliSymMatrix::~AliSymMatrix()
84{
85 Clear();
86 if (--fgCopyCnt < 1 && fgBuffer) {delete fgBuffer; fgBuffer = 0;}
87}
88
89//___________________________________________________________
90AliSymMatrix& 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
132//___________________________________________________________
7c3070ec 133AliSymMatrix& 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
144//___________________________________________________________
8a9ab0eb 145void 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//___________________________________________________________
160Float_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//___________________________________________________________
169void 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 184void 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//___________________________________________________________
196AliSymMatrix* 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//___________________________________________________________
236Bool_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//___________________________________________________________
252void 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//___________________________________________________________
291Bool_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//___________________________________________________________
330Bool_t AliSymMatrix::SolveChol(TVectorD &b, Bool_t invert)
331{
332 return SolveChol((Double_t*)b.GetMatrixArray(),invert);
333}
334
335
336//___________________________________________________________
337Bool_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 344Bool_t AliSymMatrix::SolveChol(const TVectorD &brhs, TVectorD &bsol,Bool_t invert)
8a9ab0eb 345{
346 bsol = brhs;
347 return SolveChol(bsol,invert);
348}
349
350//___________________________________________________________
351void 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//___________________________________________________________
371void 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
386//___________________________________________________________
de34b538 387/*
388void 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//___________________________________________________________
419Double_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
432//___________________________________________________________
8a9ab0eb 433int 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