STEER/AliMinResSolve.cxx

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