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