Update master to aliroot
[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 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 }