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