AliRoot Core  edcc906 (edcc906)
AliMinResSolve.cxx
Go to the documentation of this file.
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 //______________________________________________________
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 //______________________________________________________
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 //______________________________________________________
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 //______________________________________________________
63 {
64  // destructor
65  ClearAux();
66 }
67 
68 //______________________________________________________
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  //
93  if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(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 //___________________________________________________________
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 //___________________________________________________________
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  //
631  //
632  return 0;
633 }
634 
635 //___________________________________________________________
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];
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 //___________________________________________________________
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 //___________________________________________________________
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>>");
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 //___________________________________________________________
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 }
Int_t BuildPreconILUKDense(Int_t lofM)
Double_t ** fPhh
Double_t * fDiagLU
Bool_t IsSymmetric() const
Definition: AliMatrixSq.h:43
Double_t * fPVecR2
Double_t ** fPvv
Bool_t InitAuxMinRes()
static Bool_t IsZero(Double_t x, Double_t thresh=1e-64)
Definition: AliMatrixSq.h:63
TStopwatch timer
Definition: kNNSpeed.C:15
Double_t * fPVecW
Float_t GetDensity() const
Double_t * fRHS
TMatrixD mat
Definition: AnalyzeLaser.C:9
UShort_t * GetIndices() const
virtual void MultiplyByVec(const Double_t *vecIn, Double_t *vecOut) const
Definition: AliMatrixSq.cxx:33
Double_t * fPVecY
AliSymBDMatrix * fMatBD
Int_t BuildPrecon(Int_t val=0)
Double_t * fPVecW2
Double_t * fPVecW1
Double_t * fPVecR1
AliMinResSolve & operator=(const AliMinResSolve &rhs)
Bool_t SolveMinRes(Double_t *VecSol, Int_t precon=0, int itnlim=2000, double rtol=1e-12)
Int_t BuildPreconBD(Int_t hwidth)
AliVectorSparse * GetRow(Int_t ir) const
virtual Int_t GetSize() const
Definition: AliMatrixSq.h:22
#define AliInfo(message)
Definition: AliLog.h:484
Int_t PreconILUKsymbDense(Int_t lofM)
Double_t & GetElem(Int_t i) const
Bool_t SolveFGMRES(Double_t *VecSol, Int_t precon=0, int itnlim=2000, double rtol=1e-12, int nkrylov=60)
void SetSymmetric(Bool_t v=kTRUE)
Definition: AliMatrixSq.h:44
Double_t * fPVecV
Double_t ** fPvz
UShort_t & GetIndex(Int_t i)
#define AliDebug(logLevel, message)
Definition: AliLog.h:300
AliMatrixSparse * fMatU
void Solve(Double_t *rhs)
Int_t BuildPreconILUK(Int_t lofM)
AliMatrixSq * fMatrix
void ReSize(Int_t sz, Bool_t copy=kFALSE)
const Double_t kmin
class TVectorT< Double_t > TVectorD
Int_t PreconILUKsymb(Int_t lofM)
AliMatrixSparse * fMatL
void SortIndices(Bool_t valuesToo=kFALSE)
void ApplyPrecon(const TVectorD &vecRHS, TVectorD &vecOut) const
Int_t GetNElems() const
virtual Double_t Query(Int_t rown, Int_t coln) const
Definition: AliMatrixSq.h:27
Bool_t InitAuxFGMRES(int nkrylov)