STEER/AliMinResSolve.cxx

   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
   9 #include "AliMinResSolve.h"
  10 #include "AliLog.h"
  11 #include "AliMatrixSq.h"
  12 #include "AliMatrixSparse.h"
  13 #include "AliSymBDMatrix.h"
  14 #include <TStopwatch.h>
  15 #include <float.h>
  16 #include <TMath.h>
  17
  18 using namespace std;
  19
  20 ClassImp(AliMinResSolve)
  21
  22 //______________________________________________________
  23 AliMinResSolve::AliMinResSolve() :
  24 fSize(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),
  27   fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
  28 {
  29   // default constructor
  30 }
  31
  32 //______________________________________________________
  33 AliMinResSolve::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),
  38   fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
  39 {
  40   // copy constructor
  41 }
  42
  43 //______________________________________________________
  44 AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) :
  45   fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs->GetMatrixArray()),
  46   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
  47   fPvv(0),fPvz(0),fPhh(0),
  48   fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
  49 {
  50   // copy accepting equation
  51 }
  52
  53 //______________________________________________________
  54 AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) :
  55   fSize(mat->GetSize()),fPrecon(0),fMatrix((AliMatrixSq*)mat),fRHS((double*)rhs),
  56   fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0),
  57   fPvv(0),fPvz(0),fPhh(0),
  58   fDiagLU(0),fMatL(0),fMatU(0),fMatBD(0)
  59 {
  60   // copy accepting equation
  61 }
  62
  63 //______________________________________________________
  64 AliMinResSolve::~AliMinResSolve()
  65 {
  66   // destructor
  67   ClearAux();
  68 }
  69
  70 //______________________________________________________
  71 AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src)
  72 {
  73   // assignment op.
  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 //_______________________________________________________________
  84 Int_t AliMinResSolve::BuildPrecon(Int_t prec)
  85 {
  86   // preconditioner building
  87   //  const Double_t kTiny = 1E-12;
  88   fPrecon = prec;
  89   //
  90   if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
  91     return BuildPreconBD(fPrecon - kPreconBD);     // with halfbandwidth + diagonal = fPrecon
  92   }
  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
 102 //________________________________ FGMRES METHODS ________________________________
 103 Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov)
 104 {
 105   // solve by fgmres
 106   return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov);
 107 }
 108
 109 //________________________________________________________________________________
 110 Bool_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     //
 168     if (beta < epsmac) break;                                      // success?
 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;
 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}
 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
 219       if (gam < epsmac) gam = epsmac;
 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 ________________________________
 265 Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
 266 {
 267   // solve by minres
 268   return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
 269 }
 270
 271 //________________________________________________________________________________
 272 Bool_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;
 310   double normA = 0;
 311   double condA = 0;
 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   //
 324   for (int i=fSize;i--;) fPVecY[i]  = fPVecR1[i] = fRHS[i];
 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   //
 336   if (beta1 < eps) {
 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.
 459     normA  = TMath::Sqrt( tnorm2 );
 460     ynorm  = TMath::Sqrt( ynorm2 );
 461     epsa   = normA * eps;
 462     epsx   = normA * ynorm * eps;
 463     epsr   = normA * ynorm * rtol;
 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     */
 476     condA  = gmax / gmin;
 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     //
 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));
 483
 484     if (status == 0) {
 485       if (itn    >= itnlim    ) {status = 5; AliInfo(Form("%d iterations limit exceeded",itnlim));}
 486       if (condA  >= 0.1/eps   ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",condA,0.1/eps));}
 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",
 504                timer.CpuTime(),status,itn,normA,condA,rnorm,ynorm));
 505   //
 506   return status>=0 && status<=3;
 507   //
 508 }
 509
 510 //______________________________________________________________
 511 void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
 512 {
 513   // apply precond.
 514   ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
 515 }
 516
 517 //______________________________________________________________
 518 void 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
 522   //  const Double_t kTiny = 1E-12;
 523   if (fPrecon>=kPreconBD && fPrecon<kPreconILU0) { // band diagonal decomposition
 524     fMatBD->Solve(vecRHS,vecOut);
 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);
 542       vecOut[i] *= fDiagLU[i];
 543     }
 544     //
 545   }
 546   //
 547 }
 548
 549
 550 //___________________________________________________________
 551 Bool_t AliMinResSolve::InitAuxMinRes()
 552 {
 553   // init auxiliary space for minres
 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 //___________________________________________________________
 576 Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
 577 {
 578   // init auxiliary space for fgmres
 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 //___________________________________________________________
 604 void AliMinResSolve::ClearAux()
 605 {
 606   // clear aux. space
 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;
 614   if (fDiagLU)   delete[] fDiagLU; fDiagLU = 0;
 615   if (fMatL)       delete fMatL; fMatL = 0;
 616   if (fMatU)       delete fMatU; fMatU = 0;
 617   if (fMatBD)      delete fMatBD; fMatBD = 0;
 618 }
 619
 620 //___________________________________________________________
 621 Int_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;
 649 }
 650
 651 //___________________________________________________________
 652 Int_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));
 661   //
 662   TStopwatch sw; sw.Start();
 663   fMatL = new AliMatrixSparse(fSize);
 664   fMatU = new AliMatrixSparse(fSize);
 665   fMatL->SetSymmetric(kFALSE);
 666   fMatU->SetSymmetric(kFALSE);
 667   fDiagLU = new Double_t[fSize];
 668   AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
 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
 680     if ( (i%int(0.1*fSize)) == 0) {
 681       AliInfo(Form("BuildPrecon: row %d of %d",i,fSize));
 682       sw.Stop();
 683       sw.Print();
 684       sw.Start(kFALSE);
 685     }
 686     /* setup array jw[], and initial i-th row */
 687     AliVectorSparse& rowLi = *fMatL->GetRow(i);
 688     AliVectorSparse& rowUi = *fMatU->GetRow(i);
 689     AliVectorSparse& rowM  = *matrix->GetRow(i);
 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;
 697     fDiagLU[i] = 0;      // initialize diagonal
 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
 706       if (AliMatrixSq::IsZero(rowM.GetElem(j))) continue;
 707       int col = rowM.GetIndex(j);         // (the original matrix stores only lower triangle)
 708       if( col < i )   rowLi.GetElem(jw[col]) = rowM.GetElem(j);
 709       else if(col==i) fDiagLU[i] = rowM.GetElem(j);
 710       else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
 711     }
 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
 714         if (AliMatrixSq::IsZero(vl)) continue;
 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)
 722       rowLi.GetElem(j) *= fDiagLU[jrow];
 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);
 731         else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
 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     //
 740     if( AliMatrixSq::IsZero(fDiagLU[i]) ) {
 741       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
 742       delete[] jw;
 743       return -1;
 744     }
 745     fDiagLU[i] = 1.0 / fDiagLU[i];
 746   }
 747   //
 748   delete[] jw;
 749   //
 750   sw.Stop();
 751   AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
 752   AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
 753   //
 754   return 0;
 755 }
 756
 757 //___________________________________________________________
 758 Int_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);
 773   fDiagLU = new Double_t[fSize];
 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;
 795     fDiagLU[i] = 0;      // initialize diagonal
 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
 804       double vl = fMatrix->Query(i,j);
 805       if (AliMatrixSq::IsZero(vl)) continue;
 806       if( j < i )   rowLi.GetElem(jw[j]) = vl;
 807       else if(j==i) fDiagLU[i] = vl;
 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)
 814       rowLi.GetElem(j) *= fDiagLU[jrow];
 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);
 823         else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
 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     //
 832     if( AliMatrixSq::IsZero(fDiagLU[i])) {
 833       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
 834       delete[] jw;
 835       return -1;
 836     }
 837     fDiagLU[i] = 1.0 / fDiagLU[i];
 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 //___________________________________________________________
 850 Int_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   //
 858   TStopwatch sw;
 859   AliInfo("PreconILUKsymb>>");
 860   AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
 861   sw.Start();
 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;
 887     AliVectorSparse& row = *matrix->GetRow(i);
 888     //
 889     // assign lof = 0 for matrix elements
 890     for(int j=0;j<row.GetNElems(); j++) {
 891       int col = row.GetIndex(j);
 892       if (AliMatrixSq::IsZero(row.GetElem(j))) continue;  // !!!! matrix is sparse but sometimes 0 appears
 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     }
 904     if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
 905         if (AliMatrixSq::IsZero(matrix->Query(col,i))) continue;    // Due to the symmetry  == matrix(i,col)
 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();
 981   sw.Stop();
 982   sw.Print();
 983   AliInfo("PreconILUKsymb<<");
 984   return 0;
 985 }
 986
 987
 988 //___________________________________________________________
 989 Int_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++) {
1024       if (AliMatrixSq::IsZero(fMatrix->Query(i,j))) continue;
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 }