STEER/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   fPVecY   = new double[fSize];
 555   fPVecR1  = new double[fSize];
 556   fPVecR2  = new double[fSize];
 557   fPVecV   = new double[fSize];
 558   fPVecW   = new double[fSize];
 559   fPVecW1  = new double[fSize];
 560   fPVecW2  = new double[fSize];
 561   //
 562   for (int i=fSize;i--;) fPVecY[i]=fPVecR1[i]=fPVecR2[i]=fPVecV[i]=fPVecW[i]=fPVecW1[i]=fPVecW2[i]=0.0;
 563   //
 564   return kTRUE;
 565 }
 566
 567
 568 //___________________________________________________________
 569 Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov)
 570 {
 571   // init auxiliary space for fgmres
 572   fPvv     = new double*[nkrylov+1];
 573   fPvz     = new double*[nkrylov];
 574   for (int i=0; i<=nkrylov; i++) fPvv[i] = new double[fSize];
 575   fPhh     = new double*[nkrylov];
 576   for (int i=0; i<nkrylov; i++) {
 577     fPhh[i] = new double[i+2];
 578     fPvz[i]  = new double[fSize];
 579   }
 580   //
 581   fPVecR1  = new double[nkrylov];
 582   fPVecR2  = new double[nkrylov];
 583   fPVecV   = new double[nkrylov+1];
 584   //
 585   return kTRUE;
 586 }
 587
 588
 589 //___________________________________________________________
 590 void AliMinResSolve::ClearAux()
 591 {
 592   // clear aux. space
 593   if (fPVecY)      delete[] fPVecY;  fPVecY   = 0;
 594   if (fPVecR1)     delete[] fPVecR1; fPVecR1  = 0;
 595   if (fPVecR2)     delete[] fPVecR2; fPVecR2  = 0;
 596   if (fPVecV)      delete[] fPVecV;  fPVecV   = 0;
 597   if (fPVecW)      delete[] fPVecW;  fPVecW   = 0;
 598   if (fPVecW1)     delete[] fPVecW1; fPVecW1  = 0;
 599   if (fPVecW2)     delete[] fPVecW2; fPVecW2  = 0;
 600   if (fDiagLU)   delete[] fDiagLU; fDiagLU = 0;
 601   if (fMatL)       delete fMatL; fMatL = 0;
 602   if (fMatU)       delete fMatU; fMatU = 0;
 603   if (fMatBD)      delete fMatBD; fMatBD = 0;
 604 }
 605
 606 //___________________________________________________________
 607 Int_t  AliMinResSolve::BuildPreconBD(Int_t hwidth)
 608 {
 609   // build preconditioner
 610   AliInfo(Form("Building Band-Diagonal preconditioner of half-width = %d",hwidth));
 611   fMatBD = new AliSymBDMatrix( fMatrix->GetSize(), hwidth );
 612   //
 613   // fill the band-diagonal part of the matrix
 614   if (fMatrix->InheritsFrom("AliMatrixSparse")) {
 615     for (int ir=fMatrix->GetSize();ir--;) {
 616       int jmin = TMath::Max(0,ir-hwidth);
 617       AliVectorSparse& irow = *((AliMatrixSparse*)fMatrix)->GetRow(ir);
 618       for(int j=irow.GetNElems();j--;) {
 619         int jind = irow.GetIndex(j);
 620         if (jind<jmin) break;
 621         (*fMatBD)(ir,jind) = irow.GetElem(j);
 622       }
 623     }
 624   }
 625   else {
 626     for (int ir=fMatrix->GetSize();ir--;) {
 627       int jmin = TMath::Max(0,ir-hwidth);
 628       for(int jr=jmin;jr<=ir;jr++) (*fMatBD)(ir,jr) = fMatrix->Query(ir,jr);
 629     }
 630   }
 631   //
 632   fMatBD->DecomposeLDLT();
 633   //
 634   return 0;
 635 }
 636
 637 //___________________________________________________________
 638 Int_t  AliMinResSolve::BuildPreconILUK(Int_t lofM)
 639 {
 640   /*----------------------------------------------------------------------------
 641    * ILUK preconditioner
 642    * incomplete LU factorization with level of fill dropping
 643    * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
 644    *----------------------------------------------------------------------------*/
 645   //
 646   AliInfo(Form("Building ILU%d preconditioner",lofM));
 647   //
 648   TStopwatch sw; sw.Start();
 649   fMatL = new AliMatrixSparse(fSize);
 650   fMatU = new AliMatrixSparse(fSize);
 651   fMatL->SetSymmetric(kFALSE);
 652   fMatU->SetSymmetric(kFALSE);
 653   fDiagLU = new Double_t[fSize];
 654   AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
 655   //
 656   // symbolic factorization to calculate level of fill index arrays
 657   if ( PreconILUKsymb(lofM)<0 ) {
 658     ClearAux();
 659     return -1;
 660   }
 661   //
 662   Int_t *jw = new Int_t[fSize];
 663   for(int j=fSize;j--;) jw[j] = -1;     // set indicator array jw to -1
 664   //
 665   for(int i=0; i<fSize; i++ ) {            // beginning of main loop
 666     if ( (i%int(0.1*fSize)) == 0) {
 667       AliInfo(Form("BuildPrecon: row %d of %d",i,fSize));
 668       sw.Stop();
 669       sw.Print();
 670       sw.Start(kFALSE);
 671     }
 672     /* setup array jw[], and initial i-th row */
 673     AliVectorSparse& rowLi = *fMatL->GetRow(i);
 674     AliVectorSparse& rowUi = *fMatU->GetRow(i);
 675     AliVectorSparse& rowM  = *matrix->GetRow(i);
 676     //
 677     for(int j=rowLi.GetNElems(); j--;) {  // initialize L part
 678       int col = rowLi.GetIndex(j);
 679       jw[col] = j;
 680       rowLi.GetElem(j) = 0.;   // do we need this ?
 681     }
 682     jw[i] = i;
 683     fDiagLU[i] = 0;      // initialize diagonal
 684     //
 685     for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
 686       int col = rowUi.GetIndex(j);
 687       jw[col] = j;
 688       rowUi.GetElem(j) = 0;
 689     }
 690     // copy row from csmat into L,U D
 691     for(int j=rowM.GetNElems(); j--;) {  // L and D part
 692       if (AliMatrixSq::IsZero(rowM.GetElem(j))) continue;
 693       int col = rowM.GetIndex(j);         // (the original matrix stores only lower triangle)
 694       if( col < i )   rowLi.GetElem(jw[col]) = rowM.GetElem(j);
 695       else if(col==i) fDiagLU[i] = rowM.GetElem(j);
 696       else rowUi.GetElem(jw[col]) = rowM.GetElem(j);
 697     }
 698     if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) {      // part of the row I on the right of diagonal is stored as
 699         double vl = matrix->Query(col,i);    // the lower part of the column I
 700         if (AliMatrixSq::IsZero(vl)) continue;
 701         rowUi.GetElem(jw[col]) = vl;
 702       }
 703     //
 704     // eliminate previous rows
 705     for(int j=0; j<rowLi.GetNElems(); j++) {
 706       int jrow = rowLi.GetIndex(j);
 707       // get the multiplier for row to be eliminated (jrow)
 708       rowLi.GetElem(j) *= fDiagLU[jrow];
 709       //
 710       // combine current row and row jrow
 711       AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
 712       for(int k=0; k<rowUj.GetNElems(); k++ ) {
 713         int col = rowUj.GetIndex(k);
 714         int jpos = jw[col];
 715         if( jpos == -1 ) continue;
 716         if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
 717         else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
 718         else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
 719       }
 720     }
 721     // reset double-pointer to -1 ( U-part)
 722     for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
 723     jw[i] = -1;
 724     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
 725     //
 726     if( AliMatrixSq::IsZero(fDiagLU[i]) ) {
 727       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
 728       delete[] jw;
 729       return -1;
 730     }
 731     fDiagLU[i] = 1.0 / fDiagLU[i];
 732   }
 733   //
 734   delete[] jw;
 735   //
 736   sw.Stop();
 737   AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
 738   AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
 739   //
 740   return 0;
 741 }
 742
 743 //___________________________________________________________
 744 Int_t  AliMinResSolve::BuildPreconILUKDense(Int_t lofM)
 745 {
 746   /*----------------------------------------------------------------------------
 747    * ILUK preconditioner
 748    * incomplete LU factorization with level of fill dropping
 749    * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
 750    *----------------------------------------------------------------------------*/
 751   //
 752   TStopwatch sw; sw.Start();
 753   AliInfo(Form("Building ILU%d preconditioner for dense matrix",lofM));
 754   //
 755   fMatL = new AliMatrixSparse(fSize);
 756   fMatU = new AliMatrixSparse(fSize);
 757   fMatL->SetSymmetric(kFALSE);
 758   fMatU->SetSymmetric(kFALSE);
 759   fDiagLU = new Double_t[fSize];
 760   //
 761   // symbolic factorization to calculate level of fill index arrays
 762   if ( PreconILUKsymbDense(lofM)<0 ) {
 763     ClearAux();
 764     return -1;
 765   }
 766   //
 767   Int_t *jw = new Int_t[fSize];
 768   for(int j=fSize;j--;) jw[j] = -1;     // set indicator array jw to -1
 769   //
 770   for(int i=0; i<fSize; i++ ) {            // beginning of main loop
 771     /* setup array jw[], and initial i-th row */
 772     AliVectorSparse& rowLi = *fMatL->GetRow(i);
 773     AliVectorSparse& rowUi = *fMatU->GetRow(i);
 774     //
 775     for(int j=rowLi.GetNElems(); j--;) {  // initialize L part
 776       int col = rowLi.GetIndex(j);
 777       jw[col] = j;
 778       rowLi.GetElem(j) = 0.;   // do we need this ?
 779     }
 780     jw[i] = i;
 781     fDiagLU[i] = 0;      // initialize diagonal
 782     //
 783     for(int j=rowUi.GetNElems();j--;)   {  // initialize U part
 784       int col = rowUi.GetIndex(j);
 785       jw[col] = j;
 786       rowUi.GetElem(j) = 0;
 787     }
 788     // copy row from csmat into L,U D
 789     for(int j=fSize; j--;) {  // L and D part
 790       double vl = fMatrix->Query(i,j);
 791       if (AliMatrixSq::IsZero(vl)) continue;
 792       if( j < i )   rowLi.GetElem(jw[j]) = vl;
 793       else if(j==i) fDiagLU[i] = vl;
 794       else rowUi.GetElem(jw[j]) = vl;
 795     }
 796     // eliminate previous rows
 797     for(int j=0; j<rowLi.GetNElems(); j++) {
 798       int jrow = rowLi.GetIndex(j);
 799       // get the multiplier for row to be eliminated (jrow)
 800       rowLi.GetElem(j) *= fDiagLU[jrow];
 801       //
 802       // combine current row and row jrow
 803       AliVectorSparse& rowUj = *fMatU->GetRow(jrow);
 804       for(int k=0; k<rowUj.GetNElems(); k++ ) {
 805         int col = rowUj.GetIndex(k);
 806         int jpos = jw[col];
 807         if( jpos == -1 ) continue;
 808         if( col < i )   rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
 809         else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k);
 810         else            rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k);
 811       }
 812     }
 813     // reset double-pointer to -1 ( U-part)
 814     for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1;
 815     jw[i] = -1;
 816     for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1;
 817     //
 818     if( AliMatrixSq::IsZero(fDiagLU[i])) {
 819       AliInfo(Form("Fatal error in ILIk: Zero diagonal found..."));
 820       delete[] jw;
 821       return -1;
 822     }
 823     fDiagLU[i] = 1.0 / fDiagLU[i];
 824   }
 825   //
 826   delete[] jw;
 827   //
 828   sw.Stop();
 829   AliInfo(Form("ILU%d dense preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime()));
 830   //  AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity()));
 831   //
 832   return 0;
 833 }
 834
 835 //___________________________________________________________
 836 Int_t  AliMinResSolve::PreconILUKsymb(Int_t lofM)
 837 {
 838   /*----------------------------------------------------------------------------
 839    * ILUK preconditioner
 840    * incomplete LU factorization with level of fill dropping
 841    * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
 842    *----------------------------------------------------------------------------*/
 843   //
 844   TStopwatch sw;
 845   AliInfo("PreconILUKsymb>>");
 846   AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix;
 847   sw.Start();
 848   //
 849   UChar_t **ulvl=0,*levls=0;
 850   UShort_t *jbuf=0;
 851   Int_t    *iw=0;
 852   ulvl = new UChar_t*[fSize];      // stores lev-fils for U part of ILU factorization
 853   levls = new UChar_t[fSize];
 854   jbuf = new UShort_t[fSize];
 855   iw = new Int_t[fSize];
 856   //
 857   for(int j=fSize; j--;) iw[j] = -1;           // initialize iw
 858   for(int i=0; i<fSize; i++) {
 859     int incl = 0;
 860     int incu = i;
 861     AliVectorSparse& row = *matrix->GetRow(i);
 862     //
 863     // assign lof = 0 for matrix elements
 864     for(int j=0;j<row.GetNElems(); j++) {
 865       int col = row.GetIndex(j);
 866       if (AliMatrixSq::IsZero(row.GetElem(j))) continue;  // !!!! matrix is sparse but sometimes 0 appears
 867       if (col<i) {                      // L-part
 868         jbuf[incl] = col;
 869         levls[incl] = 0;
 870         iw[col] = incl++;
 871       }
 872       else if (col>i) {                 // This works only for general matrix
 873         jbuf[incu] = col;
 874         levls[incu] = 0;
 875         iw[col] = incu++;
 876       }
 877     }
 878     if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix
 879         if (AliMatrixSq::IsZero(matrix->Query(col,i))) continue;    // Due to the symmetry  == matrix(i,col)
 880         jbuf[incu] = col;
 881         levls[incu] = 0;
 882         iw[col] = incu++;
 883       }
 884     //
 885     // symbolic k,i,j Gaussian elimination
 886     int jpiv = -1;
 887     while (++jpiv < incl) {
 888       int k = jbuf[jpiv] ;                        // select leftmost pivot
 889       int kmin = k;
 890       int jmin = jpiv;
 891       for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
 892       //
 893       // ------------------------------------  swap
 894       if(jmin!=jpiv) {
 895         jbuf[jpiv] = kmin;
 896         jbuf[jmin] = k;
 897         iw[kmin] = jpiv;
 898         iw[k] = jmin;
 899         int tj = levls[jpiv] ;
 900         levls[jpiv] = levls[jmin];
 901         levls[jmin] = tj;
 902         k = kmin;
 903       }
 904       // ------------------------------------ symbolic linear combinaiton of rows
 905       AliVectorSparse& rowU = *fMatU->GetRow(k);
 906       for(int j=0; j<rowU.GetNElems(); j++ ) {
 907         int col = rowU.GetIndex(j);
 908         int it  = ulvl[k][j]+levls[jpiv]+1;
 909         if( it > lofM ) continue;
 910         int ip = iw[col];
 911         if( ip == -1 ) {
 912           if( col < i) {
 913             jbuf[incl] = col;
 914             levls[incl] = it;
 915             iw[col] = incl++;
 916           }
 917           else if( col > i ) {
 918             jbuf[incu] = col;
 919             levls[incu] = it;
 920             iw[col] = incu++;
 921           }
 922         }
 923         else levls[ip] = TMath::Min(levls[ip], it);
 924       }
 925       //
 926     } // end - while loop
 927     //
 928     // reset iw
 929     for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
 930     for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
 931     //
 932     // copy L-part
 933     AliVectorSparse& rowLi = *fMatL->GetRow(i);
 934     rowLi.ReSize(incl);
 935     if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
 936     // copy U-part
 937     int k = incu-i;
 938     AliVectorSparse& rowUi = *fMatU->GetRow(i);
 939     rowUi.ReSize(k);
 940     if( k > 0 ) {
 941       memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
 942       ulvl[i] = new UChar_t[k];   // update matrix of levels
 943       memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
 944     }
 945   }
 946   //
 947   // free temp space and leave
 948   delete[] levls;
 949   delete[] jbuf;
 950   for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i];
 951   delete[] ulvl;
 952   delete[] iw;
 953   //
 954   fMatL->SortIndices();
 955   fMatU->SortIndices();
 956   sw.Stop();
 957   sw.Print();
 958   AliInfo("PreconILUKsymb<<");
 959   return 0;
 960 }
 961
 962
 963 //___________________________________________________________
 964 Int_t  AliMinResSolve::PreconILUKsymbDense(Int_t lofM)
 965 {
 966   /*----------------------------------------------------------------------------
 967    * ILUK preconditioner
 968    * incomplete LU factorization with level of fill dropping
 969    * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/
 970    *----------------------------------------------------------------------------*/
 971   //
 972   UChar_t **ulvl=0,*levls=0;
 973   UShort_t *jbuf=0;
 974   Int_t    *iw=0;
 975   ulvl = new UChar_t*[fSize];      // stores lev-fils for U part of ILU factorization
 976   levls = new UChar_t[fSize];
 977   jbuf = new UShort_t[fSize];
 978   iw = new Int_t[fSize];
 979   //
 980   for(int j=fSize; j--;) iw[j] = -1;           // initialize iw
 981   for(int i=0; i<fSize; i++) {
 982     int incl = 0;
 983     int incu = i;
 984     //
 985     // assign lof = 0 for matrix elements
 986     for(int j=0;j<fSize; j++) {
 987       if (AliMatrixSq::IsZero(fMatrix->Query(i,j))) continue;
 988       if (j<i) {                      // L-part
 989         jbuf[incl] = j;
 990         levls[incl] = 0;
 991         iw[j] = incl++;
 992       }
 993       else if (j>i) {                 // This works only for general matrix
 994         jbuf[incu] = j;
 995         levls[incu] = 0;
 996         iw[j] = incu++;
 997       }
 998     }
 999     //
1000     // symbolic k,i,j Gaussian elimination
1001     int jpiv = -1;
1002     while (++jpiv < incl) {
1003       int k = jbuf[jpiv] ;                        // select leftmost pivot
1004       int kmin = k;
1005       int jmin = jpiv;
1006       for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;}
1007       //
1008       // ------------------------------------  swap
1009       if(jmin!=jpiv) {
1010         jbuf[jpiv] = kmin;
1011         jbuf[jmin] = k;
1012         iw[kmin] = jpiv;
1013         iw[k] = jmin;
1014         int tj = levls[jpiv] ;
1015         levls[jpiv] = levls[jmin];
1016         levls[jmin] = tj;
1017         k = kmin;
1018       }
1019       // ------------------------------------ symbolic linear combinaiton of rows
1020       AliVectorSparse& rowU = *fMatU->GetRow(k);
1021       for(int j=0; j<rowU.GetNElems(); j++ ) {
1022         int col = rowU.GetIndex(j);
1023         int it  = ulvl[k][j]+levls[jpiv]+1;
1024         if( it > lofM ) continue;
1025         int ip = iw[col];
1026         if( ip == -1 ) {
1027           if( col < i) {
1028             jbuf[incl] = col;
1029             levls[incl] = it;
1030             iw[col] = incl++;
1031           }
1032           else if( col > i ) {
1033             jbuf[incu] = col;
1034             levls[incu] = it;
1035             iw[col] = incu++;
1036           }
1037         }
1038         else levls[ip] = TMath::Min(levls[ip], it);
1039       }
1040       //
1041     } // end - while loop
1042     //
1043     // reset iw
1044     for (int j=0;j<incl;j++) iw[jbuf[j]] = -1;
1045     for (int j=i;j<incu;j++) iw[jbuf[j]] = -1;
1046     //
1047     // copy L-part
1048     AliVectorSparse& rowLi = *fMatL->GetRow(i);
1049     rowLi.ReSize(incl);
1050     if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl);
1051     // copy U-part
1052     int k = incu-i;
1053     AliVectorSparse& rowUi = *fMatU->GetRow(i);
1054     rowUi.ReSize(k);
1055     if( k > 0 ) {
1056       memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k);
1057       ulvl[i] = new UChar_t[k];   // update matrix of levels
1058       memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) );
1059     }
1060   }
1061   //
1062   // free temp space and leave
1063   delete[] levls;
1064   delete[] jbuf;
1065   for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i];
1066   delete[] ulvl;
1067   delete[] iw;
1068   //
1069   fMatL->SortIndices();
1070   fMatU->SortIndices();
1071   return 0;
1072 }