]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliMinResSolve.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Make it find libraries independently of location
[u/mrichter/AliRoot.git] / STEER / AliMinResSolve.cxx
CommitLineData
7c3070ec 1/**********************************************************************************************/
2/* General class for solving large system of linear equations */
3/* Includes MINRES, FGMRES methods as well as a few precondiotiong methods */
4/* */
5/* Author: ruben.shahoyan@cern.ch */
6/* */
7/**********************************************************************************************/
8
8a9ab0eb 9#include "AliMinResSolve.h"
10#include "AliLog.h"
7c3070ec 11#include "AliMatrixSq.h"
12#include "AliMatrixSparse.h"
13#include "AliSymBDMatrix.h"
8a9ab0eb 14#include <TStopwatch.h>
7c3070ec 15#include <float.h>
16#include <TMath.h>
8a9ab0eb 17
18using namespace std;
19
20ClassImp(AliMinResSolve)
21
22//______________________________________________________
23AliMinResSolve::AliMinResSolve() :
24fSize(0),fPrecon(0),fMatrix(0),fRHS(0),
25 fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
26 fPvv(0),fPvz(0),fPhh(0),
7c3070ec 27 fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
28{
29 // default constructor
30}
8a9ab0eb 31
32//______________________________________________________
33AliMinResSolve::AliMinResSolve(const AliMinResSolve& src) :
34 TObject(src),
35 fSize(src.fSize),fPrecon(src.fPrecon),fMatrix(src.fMatrix),fRHS(src.fRHS),
36 fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
37 fPvv(0),fPvz(0),fPhh(0),
7c3070ec 38 fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
39{
40 // copy constructor
41}
8a9ab0eb 42
43//______________________________________________________
44AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) :
7c3070ec 45 fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs->GetMatrixArray()),
8a9ab0eb 46 fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
47 fPvv(0),fPvz(0),fPhh(0),
7c3070ec 48 fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
49{
50 // copy accepting equation
51}
8a9ab0eb 52
53//______________________________________________________
54AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) :
7c3070ec 55 fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs),
8a9ab0eb 56 fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
57 fPvv(0),fPvz(0),fPhh(0),
7c3070ec 58 fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
59{
60 // copy accepting equation
61}
8a9ab0eb 62
63//______________________________________________________
64AliMinResSolve::~AliMinResSolve()
65{
7c3070ec 66 // destructor
8a9ab0eb 67 ClearAux();
68}
69
70//______________________________________________________
71AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
72{
7c3070ec 73 // assignment op.
8a9ab0eb 74 if (this != &src) {
75 fSize = src.fSize;
76 fPrecon = src.fPrecon;
77 fMatrix = src.fMatrix;
78 fRHS = src.fRHS;
79 }
80 return *this;
81}
82
83//_______________________________________________________________
84Int_t AliMinResSolve::BuildPrecon(Int_t prec)
85{
7c3070ec 86 // preconditioner building
87 // const Double_t kTiny = 1E-12;
8a9ab0eb 88 fPrecon = prec;
89 //
7c3070ec 90 if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
91 return BuildPreconBD(fPrecon - kPreconBD); // with halfbandwidth + diagonal = fPrecon
92 }
8a9ab0eb 93 //
94 if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
95 if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(fPrecon-kPreconILU0);
96 else return BuildPreconILUKDense(fPrecon-kPreconILU0);
97 }
98 //
99 return -1;
100}
101
8a9ab0eb 102//________________________________ FGMRES METHODS ________________________________
103Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
104{
7c3070ec 105 // solve by fgmres
8a9ab0eb 106 return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov);
107}
108
109//________________________________________________________________________________
110Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
111{
112 // Adapted from Y.Saad fgmrs.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
113 /*----------------------------------------------------------------------
114 | *** Preconditioned FGMRES ***
115 +-----------------------------------------------------------------------
116 | This is a simple version of the ARMS preconditioned FGMRES algorithm.
117 +-----------------------------------------------------------------------
118 | Y. S. Dec. 2000. -- Apr. 2008
119 +-----------------------------------------------------------------------
120 | VecSol = real vector of length n containing an initial guess to the
121 | precon = precondtioner id (0 = no precon)
122 | itnlim = max n of iterations
123 | rtol = tolerance for stopping iteration
124 | nkrylov = N of Krylov vectors to store
125 +---------------------------------------------------------------------*/
126 int l;
127 double status = kTRUE;
128 double t,beta,eps1=0;
129 const double epsmac = 2.22E-16;
130 //
131 AliInfo(Form("Solution by FGMRes: Preconditioner#%d Max.iter.: %d Tol.: %e NKrylov: %d",
132 precon,itnlim,rtol,nkrylov));
133 //
134 int its = 0;
135 if (nkrylov<10) {
136 AliInfo(Form("Changing N Krylov vectors from %d to %d",nkrylov,10));
137 nkrylov = 10;
138 }
139 //
140 if (precon>0) {
141 if (precon>=kPreconsTot) {
142 AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon));
143 }
144 else {
145 if (BuildPrecon(precon)<0) {
146 ClearAux();
147 AliInfo("FGMRES failed to build the preconditioner");
148 return kFALSE;
149 }
150 }
151 }
152 //
153 if (!InitAuxFGMRES(nkrylov)) return kFALSE;
154 //
155 for (l=fSize;l--;) VecSol[l] = 0;
156 //
157 //-------------------- outer loop starts here
158 TStopwatch timer; timer.Start();
159 while(1) {
160 //
161 //-------------------- compute initial residual vector
162 fMatrix->MultiplyByVec(VecSol,fPvv[0]);
163 for (l=fSize;l--;) fPvv[0][l] = fRHS[l] - fPvv[0][l]; // fPvv[0]= initial residual
164 beta = 0;
165 for (l=fSize;l--;) beta += fPvv[0][l]*fPvv[0][l];
166 beta = TMath::Sqrt(beta);
167 //
7c3070ec 168 if (beta < epsmac) break; // success?
8a9ab0eb 169 t = 1.0 / beta;
170 //-------------------- normalize: fPvv[0] = fPvv[0] / beta
171 for (l=fSize;l--;) fPvv[0][l] *= t;
172 if (its == 0) eps1 = rtol*beta;
173 //
174 // ** initialize 1-st term of rhs of hessenberg system
175 fPVecV[0] = beta;
176 int i = -1;
177 do {
178 i++;
179 its++;
180 int i1 = i+1;
181 //
182 // (Right) Preconditioning Operation z_{j} = M^{-1} v_{j}
183 //
184 if (precon>0) ApplyPrecon( fPvv[i], fPvz[i]);
185 else for (l=fSize;l--;) fPvz[i][l] = fPvv[i][l];
186 //
187 //-------------------- matvec operation w = A z_{j} = A M^{-1} v_{j}
188 fMatrix->MultiplyByVec(fPvz[i],fPvv[i1]);
189 //
190 // modified gram - schmidt...
191 // h_{i,j} = (w,v_{i})
192 // w = w - h_{i,j} v_{i}
193 //
194 for (int j=0; j<=i; j++) {
195 for (t=0, l=fSize;l--;) t+= fPvv[j][l]*fPvv[i1][l];
196 fPhh[i][j] = t;
197 for (l=fSize;l--;) fPvv[i1][l] -= t*fPvv[j][l];
198 }
199 // -------------------- h_{j+1,j} = ||w||_{2}
200 for (t=0,l=fSize;l--;) t += fPvv[i1][l]*fPvv[i1][l]; t = TMath::Sqrt(t);
201 fPhh[i][i1] = t;
7c3070ec 202 if (t > 0) for (t=1./t, l=0; l<fSize; l++) fPvv[i1][l] *= t; // v_{j+1} = w / h_{j+1,j}
8a9ab0eb 203 //
204 // done with modified gram schimdt and arnoldi step
205 // now update factorization of fPhh
206 //
207 // perform previous transformations on i-th column of h
208 //
209 for (l=1; l<=i; l++) {
210 int l1 = l-1;
211 t = fPhh[i][l1];
212 fPhh[i][l1] = fPVecR1[l1]*t + fPVecR2[l1]*fPhh[i][l];
213 fPhh[i][l] = -fPVecR2[l1]*t + fPVecR1[l1]*fPhh[i][l];
214 }
215 double gam = TMath::Sqrt( fPhh[i][i]*fPhh[i][i] + fPhh[i][i1]*fPhh[i][i1]);
216 //
217 // if gamma is zero then any small value will do...
218 // will affect only residual estimate
7c3070ec 219 if (gam < epsmac) gam = epsmac;
8a9ab0eb 220 // get next plane rotation
221 fPVecR1[i] = fPhh[i][i]/gam;
222 fPVecR2[i] = fPhh[i][i1]/gam;
223 fPVecV[i1] =-fPVecR2[i]*fPVecV[i];
224 fPVecV[i] *= fPVecR1[i];
225 //
226 // determine residual norm and test for convergence
227 fPhh[i][i] = fPVecR1[i]*fPhh[i][i] + fPVecR2[i]*fPhh[i][i1];
228 beta = TMath::Abs(fPVecV[i1]);
229 //
230 } while ( (i < nkrylov-1) && (beta > eps1) && (its < itnlim) );
231 //
232 // now compute solution. 1st, solve upper triangular system
233 fPVecV[i] = fPVecV[i]/fPhh[i][i];
234 for (int j=1; j<=i; j++) {
235 int k=i-j;
236 for (t=fPVecV[k],l=k+1; l<=i; l++) t -= fPhh[l][k]*fPVecV[l];
237 fPVecV[k] = t/fPhh[k][k];
238 }
239 // -------------------- linear combination of v[i]'s to get sol.
240 for (int j=0; j<=i; j++) for (t=fPVecV[j],l=0; l<fSize; l++) VecSol[l] += t*fPvz[j][l];
241 //
242 // -------------------- restart outer loop if needed
243 //
244 if (beta <= eps1) {
245 timer.Stop();
246 AliInfo(Form("FGMRES converged in %d iterations, CPU time: %.1f sec",its,timer.CpuTime()));
247 break; // success
248 }
249 //
250 if (its >= itnlim) {
251 timer.Stop();
252 AliInfo(Form("%d iterations limit exceeded, CPU time: %.1f sec",itnlim,timer.CpuTime()));
253 status = kFALSE;
254 break;
255 }
256 }
257 //
258 ClearAux();
259 return status;
260 //
261}
262
263
264//________________________________ MINRES METHODS ________________________________
265Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
266{
7c3070ec 267 // solve by minres
8a9ab0eb 268 return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
269}
270
271//________________________________________________________________________________
272Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double rtol)
273{
274 /*
275 Adapted from author's Fortran code:
276 Michael A. Saunders na.msaunders@na-net.ornl.gov
277
278 MINRES is an implementation of the algorithm described in the following reference:
279 C. C. Paige and M. A. Saunders (1975),
280 Solution of sparse indefinite systems of linear equations,
281 SIAM J. Numer. Anal. 12(4), pp. 617-629.
282
283 */
284 if (!fMatrix->IsSymmetric()) {
285 AliInfo("MinRes cannot solve asymmetric matrices, use FGMRes instead");
286 return kFALSE;
287 }
288 //
289 ClearAux();
290 const double eps = 2.22E-16;
291 double beta1;
292 //
293 if (precon>0) {
294 if (precon>=kPreconsTot) {
295 AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon));
296 }
297 else {
298 if (BuildPrecon(precon)<0) {
299 ClearAux();
300 AliInfo("MinRes failed to build the preconditioner");
301 return kFALSE;
302 }
303 }
304 }
305 AliInfo(Form("Solution by MinRes: Preconditioner#%d Max.iter.: %d Tol.: %e",precon,itnlim,rtol));
306 //
307 // ------------------------ initialization ---------------------->>>>
308 memset(VecSol,0,fSize*sizeof(double));
309 int status=0, itn=0;
7c3070ec 310 double normA = 0;
311 double condA = 0;
8a9ab0eb 312 double ynorm = 0;
313 double rnorm = 0;
314 double gam,gmax=1,gmin=1,gbar,oldeps,epsa,epsx,epsr,diag, delta,phi,denom,z;
315 //
316 if (!InitAuxMinRes()) return kFALSE;
317 //
318 memset(VecSol,0,fSize*sizeof(double));
319 //
320 // ------------ init aux -------------------------<<<<
321 // Set up y and v for the first Lanczos vector v1.
322 // y = beta1 P' v1, where P = C**(-1). v is really P' v1.
323 //
de34b538 324 for (int i=fSize;i--;) fPVecY[i] = fPVecR1[i] = fRHS[i];
8a9ab0eb 325 //
326 if ( precon>0 ) ApplyPrecon( fRHS, fPVecY);
327 beta1 = 0; for (int i=fSize;i--;) beta1 += fRHS[i]*fPVecY[i]; //
328 //
329 if (beta1 < 0) {
330 AliInfo(Form("Preconditioner is indefinite (init) (%e).",beta1));
331 ClearAux();
332 status = 7;
333 return kFALSE;
334 }
335 //
7c3070ec 336 if (beta1 < eps) {
8a9ab0eb 337 AliInfo(Form("RHS is zero or is the nullspace of the Preconditioner: Solution is {0}"));
338 ClearAux();
339 return kTRUE;
340 }
341 //
342 beta1 = TMath::Sqrt( beta1 ); // Normalize y to get v1 later.
343 //
344 // See if Msolve is symmetric. //RS: Skept
345 // See if Aprod is symmetric. //RS: Skept
346 //
347 double oldb = 0;
348 double beta = beta1;
349 double dbar = 0;
350 double epsln = 0;
351 double qrnorm = beta1;
352 double phibar = beta1;
353 double rhs1 = beta1;
354 double rhs2 = 0;
355 double tnorm2 = 0;
356 double ynorm2 = 0;
357 double cs = -1;
358 double sn = 0;
359 for (int i=fSize;i--;) fPVecR2[i] = fPVecR1[i];
360 //
361 TStopwatch timer; timer.Start();
362 while(status==0) { //----------------- Main iteration loop ---------------------->>>>
363 //
364 itn++;
365 /*-----------------------------------------------------------------
366 Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
367 The general iteration is similar to the case k = 1 with v0 = 0:
368 p1 = Operator * v1 - beta1 * v0,
369 alpha1 = v1'p1,
370 q2 = p2 - alpha1 * v1,
371 beta2^2 = q2'q2,
372 v2 = (1/beta2) q2.
373 Again, y = betak P vk, where P = C**(-1).
374 .... more description needed.
375 -----------------------------------------------------------------*/
376 //
377 double s = 1./beta; // Normalize previous vector (in y).
378 for (int i=fSize;i--;) fPVecV[i] = s*fPVecY[i]; // v = vk if P = I
379 //
380 fMatrix->MultiplyByVec(fPVecV,fPVecY); // APROD (VecV, VecY);
381 //
382 if (itn>=2) {
383 double btrat = beta/oldb;
384 for (int i=fSize;i--;) fPVecY[i] -= btrat*fPVecR1[i];
385 }
386 double alfa = 0; for (int i=fSize;i--;) alfa += fPVecV[i]*fPVecY[i]; // alphak
387 //
388 double alf2bt = alfa/beta;
389 for (int i=fSize;i--;) {
390 fPVecY[i] -= alf2bt*fPVecR2[i];
391 fPVecR1[i] = fPVecR2[i];
392 fPVecR2[i] = fPVecY[i];
393 }
394 //
395 if ( precon>0 ) ApplyPrecon(fPVecR2, fPVecY);
396 //
397 oldb = beta; // oldb = betak
398 beta = 0; for (int i=fSize;i--;) beta += fPVecR2[i]*fPVecY[i]; // beta = betak+1^2
399 //
400 if (beta<0) {
401 AliInfo(Form("Preconditioner is indefinite (%e)",beta));
402 status = 7;
403 break;
404 }
405 //
406 beta = TMath::Sqrt(beta); // beta = betak+1
407 tnorm2 += alfa*alfa + oldb*oldb + beta*beta;
408 //
409 if (itn == 1) { // Initialize a few things.
410 if (beta/beta1 <= 10.0*eps) {
411 status = 0; //-1 //????? beta2 = 0 or ~ 0, terminate later.
412 AliInfo("RHS is eigenvector");
413 }
414 // !tnorm2 = alfa**2
415 gmax = TMath::Abs(alfa); // alpha1
416 gmin = gmax; // alpha1
417 }
418 //
419 /*
420 Apply previous rotation Qk-1 to get
421 [deltak epslnk+1] = [cs sn][dbark 0 ]
422 [gbar k dbar k+1] [sn -cs][alfak betak+1].
423 */
424 //
425 oldeps = epsln;
426 delta = cs * dbar + sn * alfa; // delta1 = 0 deltak
427 gbar = sn * dbar - cs * alfa; // gbar 1 = alfa1 gbar k
428 epsln = sn * beta; // epsln2 = 0 epslnk+1
429 dbar = - cs * beta; // dbar 2 = beta2 dbar k+1
430 //
431 // Compute the next plane rotation Qk
432 //
433 gam = TMath::Sqrt( gbar*gbar + beta*beta ); // gammak
434 cs = gbar / gam; // ck
435 sn = beta / gam; // sk
436 phi = cs * phibar; // phik
437 phibar = sn * phibar; // phibark+1
438 //
439 // Update x.
440 denom = 1./gam;
441 //
442 for (int i=fSize;i--;) {
443 fPVecW1[i] = fPVecW2[i];
444 fPVecW2[i] = fPVecW[i];
445 fPVecW[i] = denom*(fPVecV[i]-oldeps*fPVecW1[i]-delta*fPVecW2[i]);
446 VecSol[i] += phi*fPVecW[i];
447 }
448 //
449 // Go round again.
450 //
451 gmax = TMath::Max( gmax, gam );
452 gmin = TMath::Min( gmin, gam );
453 z = rhs1 / gam;
454 ynorm2 += z*z;
455 rhs1 = rhs2 - delta * z;
456 rhs2 = - epsln * z;
457 //
458 // Estimate various norms and test for convergence.
7c3070ec 459 normA = TMath::Sqrt( tnorm2 );
8a9ab0eb 460 ynorm = TMath::Sqrt( ynorm2 );
7c3070ec 461 epsa = normA * eps;
462 epsx = normA * ynorm * eps;
463 epsr = normA * ynorm * rtol;
8a9ab0eb 464 diag = gbar;
465 if (diag == 0) diag = epsa;
466 //
467 qrnorm = phibar;
468 rnorm = qrnorm;
469 /*
470 Estimate cond(A).
471 In this version we look at the diagonals of R in the
472 factorization of the lower Hessenberg matrix, Q * H = R,
473 where H is the tridiagonal matrix from Lanczos with one
474 extra row, beta(k+1) e_k^T.
475 */
7c3070ec 476 condA = gmax / gmin;
8a9ab0eb 477 //
478 // See if any of the stopping criteria are satisfied.
479 // In rare cases, istop is already -1 from above (Abar = const*I).
480 //
f0465040 481 AliDebug(2,Form("#%5d |qnrm: %+.2e Anrm:%+.2e Cnd:%+.2e Rnrm:%+.2e Ynrm:%+.2e EpsR:%+.2e EpsX:%+.2e Beta1:%+.2e",
482 itn,qrnorm, normA,condA,rnorm,ynorm,epsr,epsx,beta1));
7c3070ec 483
8a9ab0eb 484 if (status == 0) {
485 if (itn >= itnlim ) {status = 5; AliInfo(Form("%d iterations limit exceeded",itnlim));}
7c3070ec 486 if (condA >= 0.1/eps ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",condA,0.1/eps));}
8a9ab0eb 487 if (epsx >= beta1 ) {status = 3; AliInfo(Form("Approximate convergence"));}
488 if (qrnorm <= epsx ) {status = 2; AliInfo(Form("Converged within machine precision"));}
489 if (qrnorm <= epsr ) {status = 1; AliInfo(Form("Converged"));}
490 }
491 //
492 } //----------------- Main iteration loop ----------------------<<<
493 //
494 ClearAux();
495 //
496 timer.Stop();
497 AliInfo(Form("Exit from MinRes: CPU time: %.2f sec\n"
498 "Status : %2d\n"
499 "Iterations: %4d\n"
500 "Norm : %+e\n"
501 "Condition : %+e\n"
502 "Res.Norm : %+e\n"
503 "Sol.Norm : %+e",
7c3070ec 504 timer.CpuTime(),status,itn,normA,condA,rnorm,ynorm));
8a9ab0eb 505 //
506 return status>=0 && status<=3;
507 //
508}
509
510//______________________________________________________________
511void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
512{
7c3070ec 513 // apply precond.
8a9ab0eb 514 ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
515}
516
517//______________________________________________________________
518void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const
519{
520 // Application of the preconditioner matrix:
521 // implicitly defines the matrix solving the M*VecOut = VecRHS
7c3070ec 522 // const Double_t kTiny = 1E-12;
523 if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
524 fMatBD->Solve(vecRHS,vecOut);
8a9ab0eb 525 // return;
526 }
527 //
528 else if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) {
529 //
530 for(int i=0; i<fSize; i++) { // Block L solve
531 vecOut[i] = vecRHS[i];
532 AliVectorSparse &rowLi = *fMatL->GetRow(i);
533 int n = rowLi.GetNElems();
534 for(int j=0;j<n;j++) vecOut[i] -= vecOut[ rowLi.GetIndex(j) ] * rowLi.GetElem(j);
535
536 }
537 //
538 for(int i=fSize; i--; ) { // Block -- U solve
539 AliVectorSparse &rowUi = *fMatU->GetRow(i);
540 int n = rowUi.GetNElems();
541 for(int j=0;j<n;j++ ) vecOut[i] -= vecOut[ rowUi.GetIndex(j) ] * rowUi.GetElem(j);
7c3070ec 542 vecOut[i] *= fDiagLU[i];
8a9ab0eb 543 }
544 //
545 }
546 //
547}
548
549
550//___________________________________________________________
551Bool_t AliMinResSolve::InitAuxMinRes()
552{
7c3070ec 553 // init auxiliary space for minres
8a9ab0eb 554 try {
555 fPVecY = new double[fSize];
556 fPVecR1 = new double[fSize];
557 fPVecR2 = new double[fSize];
558 fPVecV = new double[fSize];
559 fPVecW = new double[fSize];
560 fPVecW1 = new double[fSize];
561 fPVecW2 = new double[fSize];
562 }
563 catch(bad_alloc&) {
564 AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize));
565 ClearAux();
566 return kFALSE;
567 }
568 //
569 for (int i=fSize;i--;) fPVecY[i]=fPVecR1[i]=fPVecR2[i]=fPVecV[i]=fPVecW[i]=fPVecW1[i]=fPVecW2[i]=0.0;
570 //
571 return kTRUE;
572}
573
574
575//___________________________________________________________
576Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
577{
7c3070ec 578 // init auxiliary space for fgmres
8a9ab0eb 579 try {
580 fPvv = new double*[nkrylov+1];
581 fPvz = new double*[nkrylov];
582 for (int i=0; i<=nkrylov; i++) fPvv[i] = new double[fSize];
583 fPhh = new double*[nkrylov];
584 for (int i=0; i<nkrylov; i++) {
585 fPhh[i] = new double[i+2];
586 fPvz[i] = new double[fSize];
587 }
588 //
589 fPVecR1 = new double[nkrylov];
590 fPVecR2 = new double[nkrylov];
591 fPVecV = new double[nkrylov+1];
592 }
593 catch(bad_alloc&) {
594 AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize));
595 ClearAux();
596 return kFALSE;
597 }
598 //
599 return kTRUE;
600}
601
602
603//___________________________________________________________
604void AliMinResSolve::ClearAux()
605{
7c3070ec 606 // clear aux. space
8a9ab0eb 607 if (fPVecY) delete[] fPVecY; fPVecY = 0;
608 if (fPVecR1) delete[] fPVecR1; fPVecR1 = 0;
609 if (fPVecR2) delete[] fPVecR2; fPVecR2 = 0;
610 if (fPVecV) delete[] fPVecV; fPVecV = 0;
611 if (fPVecW) delete[] fPVecW; fPVecW = 0;
612 if (fPVecW1) delete[] fPVecW1; fPVecW1 = 0;
613 if (fPVecW2) delete[] fPVecW2; fPVecW2 = 0;
7c3070ec 614 if (fDiagLU) delete[] fDiagLU; fDiagLU = 0;
8a9ab0eb 615 if (fMatL) delete fMatL; fMatL = 0;
616 if (fMatU) delete fMatU; fMatU = 0;
7c3070ec 617 if (fMatBD) delete fMatBD; fMatBD = 0;
618}
619
620//___________________________________________________________
621Int_t AliMinResSolve::BuildPreconBD(Int_t hwidth)
622{
623 // build preconditioner
624 AliInfo(Form("Building Band-Diagonal preconditioner of half-width = %d",hwidth));
625 fMatBD = new AliSymBDMatrix( fMatrix->GetSize(), hwidth );
626 //
627 // fill the band-diagonal part of the matrix
628 if (fMatrix->InheritsFrom("AliMatrixSparse")) {
629 for (int ir=fMatrix->GetSize();ir--;) {
630 int jmin = TMath::Max(0,ir-hwidth);
631 AliVectorSparse& irow = *((AliMatrixSparse*)fMatrix)->GetRow(ir);
632 for(int j=irow.GetNElems();j--;) {
633 int jind = irow.GetIndex(j);
634 if (jind<jmin) break;
635 (*fMatBD)(ir,jind) = irow.GetElem(j);
636 }
637 }
638 }
639 else {
640 for (int ir=fMatrix->GetSize();ir--;) {
641 int jmin = TMath::Max(0,ir-hwidth);
642 for(int jr=jmin;jr<=ir;jr++) (*fMatBD)(ir,jr) = fMatrix->Query(ir,jr);
643 }
644 }
645 //
646 fMatBD->DecomposeLDLT();
647 //
648 return 0;
8a9ab0eb 649}
650
651//___________________________________________________________
652Int_t AliMinResSolve::BuildPreconILUK(Int_t lofM)
653{
654 /*----------------------------------------------------------------------------
655 * ILUK preconditioner
656 * incomplete LU factorization with level of fill dropping
657 * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
658 *----------------------------------------------------------------------------*/
659 //
660 AliInfo(Form("Building ILU%d preconditioner",lofM));
8a9ab0eb 661 //
662 TStopwatch sw; sw.Start();
663 fMatL = new AliMatrixSparse(fSize);
664 fMatU = new AliMatrixSparse(fSize);
665 fMatL->SetSymmetric(kFALSE);
666 fMatU->SetSymmetric(kFALSE);
7c3070ec 667 fDiagLU = new Double_t[fSize];
668 AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
8a9ab0eb 669 //
670 // symbolic factorization to calculate level of fill index arrays
671 if ( PreconILUKsymb(lofM)<0 ) {
672 ClearAux();
673 return -1;
674 }
675 //
676 Int_t *jw = new Int_t[fSize];
677 for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1
678 //
679 for(int i=0; i<fSize; i++ ) { // beginning of main loop
de34b538 680 if ( (i%int(0.1*fSize)) == 0) {
7c3070ec 681 AliInfo(Form("BuildPrecon: row %d of %d",i,fSize));
de34b538 682 sw.Stop();
683 sw.Print();
684 sw.Start(kFALSE);
685 }
8a9ab0eb 686 /* setup array jw[], and initial i-th row */
687 AliVectorSparse& rowLi = *fMatL->GetRow(i);
688 AliVectorSparse& rowUi = *fMatU->GetRow(i);
7c3070ec 689 AliVectorSparse& rowM = *matrix->GetRow(i);
8a9ab0eb 690 //
691 for(int j=rowLi.GetNElems(); j--;) { // initialize L part
692 int col = rowLi.GetIndex(j);
693 jw[col] = j;
694 rowLi.GetElem(j) = 0.; // do we need this ?
695 }
696 jw[i] = i;
7c3070ec 697 fDiagLU[i] = 0; // initialize diagonal
8a9ab0eb 698 //
699 for(int j=rowUi.GetNElems();j--;) { // initialize U part
700 int col = rowUi.GetIndex(j);
701 jw[col] = j;
702 rowUi.GetElem(j) = 0;
703 }
704 // copy row from csmat into L,U D
705 for(int j=rowM.GetNElems(); j--;) { // L and D part
7c185459 706 if (AliMatrixSq::IsZero(rowM.GetElem(j))) continue;
8a9ab0eb 707 int col = rowM.GetIndex(j); // (the original matrix stores only lower triangle)
708 if( col < i ) rowLi.GetElem(jw[col]) = rowM.GetElem(j);
7c3070ec 709 else if(col==i) fDiagLU[i] = rowM.GetElem(j);
8a9ab0eb 710 else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
711 }
7c3070ec 712 if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // part of the row I on the right of diagonal is stored as
713 double vl = matrix->Query(col,i); // the lower part of the column I
7c185459 714 if (AliMatrixSq::IsZero(vl)) continue;
8a9ab0eb 715 rowUi.GetElem(jw[col]) = vl;
716 }
717 //
718 // eliminate previous rows
719 for(int j=0; j<rowLi.GetNElems(); j++) {
720 int jrow = rowLi.GetIndex(j);
721 // get the multiplier for row to be eliminated (jrow)
7c3070ec 722 rowLi.GetElem(j) *= fDiagLU[jrow];
8a9ab0eb 723 //
724 // combine current row and row jrow
725 AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
726 for(int k=0; k<rowUj.GetNElems(); k++ ) {
727 int col = rowUj.GetIndex(k);
728 int jpos = jw[col];
729 if( jpos == -1 ) continue;
730 if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
7c3070ec 731 else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
8a9ab0eb 732 else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
733 }
734 }
735 // reset double-pointer to -1 ( U-part)
736 for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
737 jw[i] = -1;
738 for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
739 //
7c185459 740 if( AliMatrixSq::IsZero(fDiagLU[i]) ) {
8a9ab0eb 741 AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
742 delete[] jw;
743 return -1;
744 }
7c3070ec 745 fDiagLU[i] = 1.0 / fDiagLU[i];
8a9ab0eb 746 }
747 //
748 delete[] jw;
749 //
750 sw.Stop();
751 AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
7c3070ec 752 AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
8a9ab0eb 753 //
754 return 0;
755}
756
757//___________________________________________________________
758Int_t AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
759{
760 /*----------------------------------------------------------------------------
761 * ILUK preconditioner
762 * incomplete LU factorization with level of fill dropping
763 * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
764 *----------------------------------------------------------------------------*/
765 //
766 TStopwatch sw; sw.Start();
767 AliInfo(Form("Building ILU%d preconditioner for dense matrix",lofM));
768 //
769 fMatL = new AliMatrixSparse(fSize);
770 fMatU = new AliMatrixSparse(fSize);
771 fMatL->SetSymmetric(kFALSE);
772 fMatU->SetSymmetric(kFALSE);
7c3070ec 773 fDiagLU = new Double_t[fSize];
8a9ab0eb 774 //
775 // symbolic factorization to calculate level of fill index arrays
776 if ( PreconILUKsymbDense(lofM)<0 ) {
777 ClearAux();
778 return -1;
779 }
780 //
781 Int_t *jw = new Int_t[fSize];
782 for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1
783 //
784 for(int i=0; i<fSize; i++ ) { // beginning of main loop
785 /* setup array jw[], and initial i-th row */
786 AliVectorSparse& rowLi = *fMatL->GetRow(i);
787 AliVectorSparse& rowUi = *fMatU->GetRow(i);
788 //
789 for(int j=rowLi.GetNElems(); j--;) { // initialize L part
790 int col = rowLi.GetIndex(j);
791 jw[col] = j;
792 rowLi.GetElem(j) = 0.; // do we need this ?
793 }
794 jw[i] = i;
7c3070ec 795 fDiagLU[i] = 0; // initialize diagonal
8a9ab0eb 796 //
797 for(int j=rowUi.GetNElems();j--;) { // initialize U part
798 int col = rowUi.GetIndex(j);
799 jw[col] = j;
800 rowUi.GetElem(j) = 0;
801 }
802 // copy row from csmat into L,U D
803 for(int j=fSize; j--;) { // L and D part
de34b538 804 double vl = fMatrix->Query(i,j);
7c185459 805 if (AliMatrixSq::IsZero(vl)) continue;
8a9ab0eb 806 if( j < i ) rowLi.GetElem(jw[j]) = vl;
7c3070ec 807 else if(j==i) fDiagLU[i] = vl;
8a9ab0eb 808 else rowUi.GetElem(jw[j]) = vl;
809 }
810 // eliminate previous rows
811 for(int j=0; j<rowLi.GetNElems(); j++) {
812 int jrow = rowLi.GetIndex(j);
813 // get the multiplier for row to be eliminated (jrow)
7c3070ec 814 rowLi.GetElem(j) *= fDiagLU[jrow];
8a9ab0eb 815 //
816 // combine current row and row jrow
817 AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
818 for(int k=0; k<rowUj.GetNElems(); k++ ) {
819 int col = rowUj.GetIndex(k);
820 int jpos = jw[col];
821 if( jpos == -1 ) continue;
822 if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
7c3070ec 823 else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
8a9ab0eb 824 else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
825 }
826 }
827 // reset double-pointer to -1 ( U-part)
828 for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
829 jw[i] = -1;
830 for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
831 //
7c185459 832 if( AliMatrixSq::IsZero(fDiagLU[i])) {
8a9ab0eb 833 AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
834 delete[] jw;
835 return -1;
836 }
7c3070ec 837 fDiagLU[i] = 1.0 / fDiagLU[i];
8a9ab0eb 838 }
839 //
840 delete[] jw;
841 //
842 sw.Stop();
843 AliInfo(Form("ILU%d dense preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
844 // AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
845 //
846 return 0;
847}
848
849//___________________________________________________________
850Int_t AliMinResSolve::PreconILUKsymb(Int_t lofM)
851{
852 /*----------------------------------------------------------------------------
853 * ILUK preconditioner
854 * incomplete LU factorization with level of fill dropping
855 * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
856 *----------------------------------------------------------------------------*/
857 //
de34b538 858 TStopwatch sw;
7c3070ec 859 AliInfo("PreconILUKsymb>>");
860 AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
de34b538 861 sw.Start();
8a9ab0eb 862 //
863 UChar_t **ulvl=0,*levls=0;
864 UShort_t *jbuf=0;
865 Int_t *iw=0;
866 try {
867 ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization
868 levls = new UChar_t[fSize];
869 jbuf = new UShort_t[fSize];
870 iw = new Int_t[fSize];
871 }
872 //
873 catch(bad_alloc&) {
874 AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb");
875 if (ulvl) delete[] ulvl;
876 if (levls) delete[] levls;
877 if (jbuf) delete[] jbuf;
878 if (iw) delete[] iw;
879 ClearAux();
880 return -1;
881 }
882 //
883 for(int j=fSize; j--;) iw[j] = -1; // initialize iw
884 for(int i=0; i<fSize; i++) {
885 int incl = 0;
886 int incu = i;
7c3070ec 887 AliVectorSparse& row = *matrix->GetRow(i);
8a9ab0eb 888 //
889 // assign lof = 0 for matrix elements
890 for(int j=0;j<row.GetNElems(); j++) {
891 int col = row.GetIndex(j);
7c185459 892 if (AliMatrixSq::IsZero(row.GetElem(j))) continue; // !!!! matrix is sparse but sometimes 0 appears
8a9ab0eb 893 if (col<i) { // L-part
894 jbuf[incl] = col;
895 levls[incl] = 0;
896 iw[col] = incl++;
897 }
898 else if (col>i) { // This works only for general matrix
899 jbuf[incu] = col;
900 levls[incu] = 0;
901 iw[col] = incu++;
902 }
903 }
7c3070ec 904 if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
7c185459 905 if (AliMatrixSq::IsZero(matrix->Query(col,i))) continue; // Due to the symmetry == matrix(i,col)
8a9ab0eb 906 jbuf[incu] = col;
907 levls[incu] = 0;
908 iw[col] = incu++;
909 }
910 //
911 // symbolic k,i,j Gaussian elimination
912 int jpiv = -1;
913 while (++jpiv < incl) {
914 int k = jbuf[jpiv] ; // select leftmost pivot
915 int kmin = k;
916 int jmin = jpiv;
917 for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
918 //
919 // ------------------------------------ swap
920 if(jmin!=jpiv) {
921 jbuf[jpiv] = kmin;
922 jbuf[jmin] = k;
923 iw[kmin] = jpiv;
924 iw[k] = jmin;
925 int tj = levls[jpiv] ;
926 levls[jpiv] = levls[jmin];
927 levls[jmin] = tj;
928 k = kmin;
929 }
930 // ------------------------------------ symbolic linear combinaiton of rows
931 AliVectorSparse& rowU = *fMatU->GetRow(k);
932 for(int j=0; j<rowU.GetNElems(); j++ ) {
933 int col = rowU.GetIndex(j);
934 int it = ulvl[k][j]+levls[jpiv]+1;
935 if( it > lofM ) continue;
936 int ip = iw[col];
937 if( ip == -1 ) {
938 if( col < i) {
939 jbuf[incl] = col;
940 levls[incl] = it;
941 iw[col] = incl++;
942 }
943 else if( col > i ) {
944 jbuf[incu] = col;
945 levls[incu] = it;
946 iw[col] = incu++;
947 }
948 }
949 else levls[ip] = TMath::Min(levls[ip], it);
950 }
951 //
952 } // end - while loop
953 //
954 // reset iw
955 for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
956 for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
957 //
958 // copy L-part
959 AliVectorSparse& rowLi = *fMatL->GetRow(i);
960 rowLi.ReSize(incl);
961 if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
962 // copy U-part
963 int k = incu-i;
964 AliVectorSparse& rowUi = *fMatU->GetRow(i);
965 rowUi.ReSize(k);
966 if( k > 0 ) {
967 memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
968 ulvl[i] = new UChar_t[k]; // update matrix of levels
969 memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
970 }
971 }
972 //
973 // free temp space and leave
974 delete[] levls;
975 delete[] jbuf;
976 for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i];
977 delete[] ulvl;
978 //
979 fMatL->SortIndices();
980 fMatU->SortIndices();
de34b538 981 sw.Stop();
982 sw.Print();
7c3070ec 983 AliInfo("PreconILUKsymb<<");
8a9ab0eb 984 return 0;
985}
986
987
988//___________________________________________________________
989Int_t AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
990{
991 /*----------------------------------------------------------------------------
992 * ILUK preconditioner
993 * incomplete LU factorization with level of fill dropping
994 * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
995 *----------------------------------------------------------------------------*/
996 //
997 UChar_t **ulvl=0,*levls=0;
998 UShort_t *jbuf=0;
999 Int_t *iw=0;
1000 try {
1001 ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization
1002 levls = new UChar_t[fSize];
1003 jbuf = new UShort_t[fSize];
1004 iw = new Int_t[fSize];
1005 }
1006 //
1007 catch(bad_alloc&) {
1008 AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb");
1009 if (ulvl) delete[] ulvl;
1010 if (levls) delete[] levls;
1011 if (jbuf) delete[] jbuf;
1012 if (iw) delete[] iw;
1013 ClearAux();
1014 return -1;
1015 }
1016 //
1017 for(int j=fSize; j--;) iw[j] = -1; // initialize iw
1018 for(int i=0; i<fSize; i++) {
1019 int incl = 0;
1020 int incu = i;
1021 //
1022 // assign lof = 0 for matrix elements
1023 for(int j=0;j<fSize; j++) {
7c185459 1024 if (AliMatrixSq::IsZero(fMatrix->Query(i,j))) continue;
8a9ab0eb 1025 if (j<i) { // L-part
1026 jbuf[incl] = j;
1027 levls[incl] = 0;
1028 iw[j] = incl++;
1029 }
1030 else if (j>i) { // This works only for general matrix
1031 jbuf[incu] = j;
1032 levls[incu] = 0;
1033 iw[j] = incu++;
1034 }
1035 }
1036 //
1037 // symbolic k,i,j Gaussian elimination
1038 int jpiv = -1;
1039 while (++jpiv < incl) {
1040 int k = jbuf[jpiv] ; // select leftmost pivot
1041 int kmin = k;
1042 int jmin = jpiv;
1043 for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
1044 //
1045 // ------------------------------------ swap
1046 if(jmin!=jpiv) {
1047 jbuf[jpiv] = kmin;
1048 jbuf[jmin] = k;
1049 iw[kmin] = jpiv;
1050 iw[k] = jmin;
1051 int tj = levls[jpiv] ;
1052 levls[jpiv] = levls[jmin];
1053 levls[jmin] = tj;
1054 k = kmin;
1055 }
1056 // ------------------------------------ symbolic linear combinaiton of rows
1057 AliVectorSparse& rowU = *fMatU->GetRow(k);
1058 for(int j=0; j<rowU.GetNElems(); j++ ) {
1059 int col = rowU.GetIndex(j);
1060 int it = ulvl[k][j]+levls[jpiv]+1;
1061 if( it > lofM ) continue;
1062 int ip = iw[col];
1063 if( ip == -1 ) {
1064 if( col < i) {
1065 jbuf[incl] = col;
1066 levls[incl] = it;
1067 iw[col] = incl++;
1068 }
1069 else if( col > i ) {
1070 jbuf[incu] = col;
1071 levls[incu] = it;
1072 iw[col] = incu++;
1073 }
1074 }
1075 else levls[ip] = TMath::Min(levls[ip], it);
1076 }
1077 //
1078 } // end - while loop
1079 //
1080 // reset iw
1081 for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
1082 for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
1083 //
1084 // copy L-part
1085 AliVectorSparse& rowLi = *fMatL->GetRow(i);
1086 rowLi.ReSize(incl);
1087 if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
1088 // copy U-part
1089 int k = incu-i;
1090 AliVectorSparse& rowUi = *fMatU->GetRow(i);
1091 rowUi.ReSize(k);
1092 if( k > 0 ) {
1093 memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
1094 ulvl[i] = new UChar_t[k]; // update matrix of levels
1095 memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
1096 }
1097 }
1098 //
1099 // free temp space and leave
1100 delete[] levls;
1101 delete[] jbuf;
1102 for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i];
1103 delete[] ulvl;
1104 //
1105 fMatL->SortIndices();
1106 fMatU->SortIndices();
1107 return 0;
1108}