Removing try/catch-es
[u/mrichter/AliRoot.git] / 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 }