EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
ratlanNr.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <string.h>
5 #include <float.h>
6 #include <complex.h>
7 #include "def.h"
8 #include "blaslapack.h"
9 #include "struct.h"
10 #include "internal_proto.h"
18 #define FILTER_VINIT 1
19 
54 int RatLanNr(double *intv, int maxit, double tol, double *vinit,
55  ratparams *rat, int *nevOut, double **lamo, double **Wo,
56  double **reso, FILE *fstats) {
57  //-------------------- to report timings/
58  double tall, tm1 = 0.0, tt;
59  tall = evsl_timer();
60  const int ifGenEv = evsldata.ifGenEv;
61  double tr0, tr1;
62  double *y, flami;
63  int i, k, kdim;
64  // handle case where fstats is NULL. Then no output. Needed for openMP.
65  int do_print = 1;
66  if (fstats == NULL){
67  do_print = 0;
68  }
69  /*-------------------- frequently used constants */
70  char cN = 'N';
71  int one = 1;
72  double done = 1.0, dzero = 0.0;
73  /*-------------------- Ntest = when to start testing convergence */
74  int Ntest = 30;
75  /*-------------------- how often to test */
76  int cycle = 20;
77  /* size of the matrix */
78  int n = evsldata.n;
79  /* max num of its */
80  maxit = min(n, maxit);
81  /*-------------------- Caveat !!!: To prevent integer overflow, we save
82  * another size_t version of n
83  * Note n itself may not overflow but something like
84  * (maxit * n) is much more likely
85  * When using size_t n, in all the multiplications
86  * integer is promoted to size_t */
87  size_t n_l = n;
88  /*-------------------- Rational filter with pole at ((a+b)/2,(b-a)/2) with
89  multiplicity pow, bar value equals 1/2 */
90  /*-------------------- a, b, used for testing only at end */
91  double bar = 0.5;
92  if (check_intv(intv, fstats) < 0) {
93  *nevOut = 0;
94  *lamo = NULL; *Wo = NULL; *reso = NULL;
95  return 0;
96  }
97  double aa = intv[0];
98  double bb = intv[1];
99  //int deg = rat->pow; // multiplicity of the pole
100  /*-----------------------------------------------------------------------*
101  * *Non-restarted* Lanczos iteration
102  *-----------------------------------------------------------------------*/
103  if (do_print) {
104  fprintf(fstats, " ** Rat-LanNr \n");
105  }
106  /*-------------------- Lanczos vectors V_m and tridiagonal matrix T_m */
107  double *V, *dT, *eT, *Z;
108  Malloc(V, n_l*(maxit+1), double);
109  if (ifGenEv) {
110  /* storage for Z = B * V */
111  Malloc(Z, n_l*(maxit+1), double);
112  } else {
113  /* Z and V are the same */
114  Z = V;
115  }
116  /*-------------------- diag. subdiag of Tridiagional matrix */
117  Malloc(dT, maxit, double);
118  Malloc(eT, maxit, double);
119  double *Rvec, *Lam, *res, *EvalT, *EvecT;
120  /*-------------------- Lam, Rvec: the converged (locked) Ritz values vecs*/
121  Malloc(Lam, maxit, double); // holds computed Ritz values
122  Malloc(res, maxit, double); // residual norms (w.r.t. ro(A))
123  Malloc(EvalT, maxit, double); // eigenvalues of tridia. matrix T
124  /*-------------------- nev = current number of converged e-pairs
125  nconv = converged eigenpairs from looking at Tk alone */
126  int nev, nconv = 0;
127  /*-------------------- u is just a pointer. wk == work space */
128  double *u, *wk, *w2, *vrand = NULL;
129  size_t wk_size = ifGenEv ? 6*n_l : 4*n_l;
130  Malloc(wk, wk_size, double);
131  w2 = wk + n;
132  /*-------------------- copy initial vector to Z(:,1) */
133 #if FILTER_VINIT
134  /* Filter the initial vector */
135  RatFiltApply(n, rat, vinit, V, wk);
136  Malloc(vrand, n, double);
137  if(ifGenEv){
138  DCOPY(&n, V, &one, Z, &one);
139  }
140 #else
141  DCOPY(&n, vinit, &one, Z, &one);
142 #endif
143  /*-------------------- normalize it */
144  double t, nt, res0;
145  if (ifGenEv) {
146  /* B norm */
147  matvec_B(Z, V);
148  t = 1.0 / sqrt(DDOT(&n, V, &one, Z, &one));
149  DSCAL(&n, &t, Z, &one);
150  } else {
151  /* 2-norm */
152  t = 1.0 / DNRM2(&n, V, &one); // add a test here.
153  }
154  /* unit B^{-1}-norm or 2-norm */
155  DSCAL(&n, &t, V, &one);
156  /*-------------------- for ortho test */
157  double wn = 0.0;
158  int nwn = 0;
159  /*-------------------- for stopping test [restricted trace]*/
160  tr0 = 0;
161  /*-------------------- lanczos vectors updated by rotating pointer*/
162  /*-------------------- pointers to Lanczos vectors */
163  double *zold, *z, *znew;
164  double *v, *vnew;
165  /*-------------------- Lanczos recurrence coefficients */
166  double alpha, nalpha, beta=0.0, nbeta;
167  int count = 0;
168  // ---------------- main Lanczos loop
169  for (k=0; k<maxit; k++) {
170  /*-------------------- quick reference to Z(:,k-1) when k>0 */
171  zold = k > 0 ? Z+(k-1)*n_l : NULL;
172  /*-------------------- a quick reference to V(:,k) */
173  v = &V[k*n_l];
174  /*-------------------- a quick reference to Z(:,k) */
175  z = &Z[k*n_l];
176  /*-------------------- next Lanczos vector V(:,k+1)*/
177  vnew = v + n;
178  /*-------------------- next Lanczos vector Z(:,k+1)*/
179  znew = z + n;
180  /*-------------------- compute w = rat * v */
181  RatFiltApply(n, rat, v, znew, wk);
182  /*------------------ znew = znew - beta*zold */
183  if (zold) {
184  nbeta = -beta;
185  DAXPY(&n, &nbeta, zold, &one, znew, &one);
186  }
187  /*-------------------- alpha = znew'*v */
188  alpha = DDOT(&n, v, &one, znew, &one);
189  /*-------------------- T(k,k) = alpha */
190  dT[k] = alpha;
191  wn += fabs(alpha);
192  /*-------------------- znew = znew - alpha*z */
193  nalpha = -alpha;
194  DAXPY(&n, &nalpha, z, &one, znew, &one);
195  /*-------------------- FULL reortho to all previous Lan vectors */
196  if (ifGenEv) {
197  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
198  CGS_DGKS2(n, k+1, NGS_MAX, Z, V, znew, wk);
199  /* -------------- NOTE: B-matvec
200  * vnew = B * znew */
201  matvec_B(znew, vnew);
202  /*-------------------- beta = (vnew, znew)^{1/2} */
203  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
204  } else {
205  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
206  /* beta = norm(vnew) */
207  CGS_DGKS(n, k+1, NGS_MAX, V, vnew, &beta, wk);
208  }
209  wn += 2.0 * beta;
210  nwn += 3;
211  /*-------------------- lucky breakdown test */
212  if (beta*nwn < orthTol*wn) {
213  if (do_print) {
214  fprintf(fstats, "it %4d: Lucky breakdown, beta = %.15e\n", k, beta);
215  }
216 #if FILTER_VINIT
217  /*------------------ generate a new init vector in znew */
218  rand_double(n, vrand);
219  /* Filter the initial vector */
220  RatFiltApply(n, rat, vrand, znew, wk);
221 #else
222  rand_double(n, znew);
223 #endif
224  if (ifGenEv) {
225  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
226  CGS_DGKS2(n, k+1, NGS_MAX, Z, V, znew, wk);
227  /* -------------- NOTE: B-matvec */
228  matvec_B(znew, vnew);
229  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
230  /*-------------------- vnew = vnew / beta */
231  t = 1.0 / beta;
232  DSCAL(&n, &t, vnew, &one);
233  /*-------------------- znew = znew / beta */
234  DSCAL(&n, &t, znew, &one);
235  beta = 0.0;
236  } else {
237  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
238  /* beta = norm(vnew) */
239  CGS_DGKS(n, k+1, NGS_MAX, V, vnew, &beta, wk);
240  /*-------------------- vnew = vnew / beta */
241  t = 1.0 / beta;
242  DSCAL(&n, &t, vnew, &one);
243  beta = 0.0;
244  }
245  } else {
246  /*-------------------- vnew = vnew / beta */
247  t = 1.0 / beta;
248  DSCAL(&n, &t, vnew, &one);
249  if (ifGenEv) {
250  /*-------------------- znew = znew / beta */
251  DSCAL(&n, &t, znew, &one);
252  }
253  }
254  /*-------------------- T(k+1,k) = beta */
255  eT[k] = beta;
256 #if 0
257  /*-------------------- Reallocate memory if maxit is smaller than # of eigs */
258  if (k == maxit-1) {
259  maxit = 1 + (int) (maxit * 1.5);
260  Realloc(V, (maxit+1)*n_l, double);
261  if (ifGenEv) {
262  Realloc(Z, (maxit+1)*n_l, double);
263  } else {
264  /* make sure Z == V since V may be changed in the Realloc above */
265  Z = V;
266  }
267  Realloc(dT, maxit, double);
268  Realloc(eT, maxit, double);
269  Realloc(Lam, maxit, double);
270  Realloc(res, maxit, double);
271  Realloc(EvalT, maxit, double);
272  }
273 #endif
274  /*---------------------- test for Ritz vectors */
275  if ( (k < Ntest || (k-Ntest) % cycle != 0) && k != maxit-1 ) {
276  continue;
277  }
278  /*---------------------- diagonalize T(1:k,1:k)
279  vals in EvalT, vecs in EvecT */
280  kdim = k+1;
281 #if 1
282  /*-------------------- THIS uses dsetv, do not need eig vector */
283  SymmTridEig(EvalT, NULL, kdim, dT, eT);
284  count = kdim;
285 #else
286  /*-------------------- THIS uses dstemr */
287  double vl = bar - DBL_EPSILON, vu = 10000.0; /* needed by SymmTridEigS */
288  SymmTridEigS(EvalT, EvecT, kdim, vl, vu, &count, dT, eT);
289 #endif
290  /*-------------------- restricted trace: used for convergence test */
291  tr1 = 0;
292  /*-------------------- get residual norms and check acceptance of Ritz
293  * values for r(A). nconv records number of
294  * eigenvalues whose residual for r(A) is smaller
295  * than tol. */
296  nconv = 0;
297  for (i=0; i<count; i++) {
298  flami = EvalT[i];
299  if (flami + DBL_EPSILON >= bar) {
300  tr1+= flami;
301  nconv++;
302  }
303  }
304 
305  if (do_print) {
306  fprintf(fstats, "k %4d: nconv %4d tr1 %21.15e\n",
307  k, nconv,tr1);
308  }
309  /* -------------------- simple test because all eigenvalues
310  are between gamB and ~1. */
311  if ( (fabs(tr1-tr0) < 2e-12) || (fabs(tr1)+fabs(tr0)<1e-10) ) {
312  break;
313  }
314  tr0 = tr1;
315  } /* end of the main loop */
316 
317  if (k >= maxit) {
318  fprintf(fstats, "The max number of iterations [%d] has been reached. The eigenvalues computed may not have converged\n", maxit);
319  }
320 
321  /*-------------------- compute eig vals and vector */
322  size_t kdim_l = kdim; /* just in case if kdim is > 65K */
323  Malloc(EvecT, kdim_l*kdim_l, double); // Eigen vectors of T
324  SymmTridEig(EvalT, EvecT, kdim, dT, eT);
325 
326  tt = evsl_timer();
327  /*-------------------- done == compute Ritz vectors */
328  Malloc(Rvec, nconv*n_l, double); // holds computed Ritz vectors
329 
330  nev = 0;
331  for (i=0; i<count; i++) {
332  flami = EvalT[i];
333  //-------------------- reject eigenvalue if rho(lam)<bar
334  if (fabs(flami) < bar) {
335  continue;
336  }
337  y = &EvecT[i*kdim_l];
338  /*-------------------- make sure to normalize */
339  /*
340  t = DNRM2(&kdim, y, &one);
341  t = 1.0 / t;
342  DSCAL(&kdim, &t, y, &one);
343  */
344  /*-------------------- compute Ritz vectors
345  * NOTE: use Z for gen e.v */
346  u = &Rvec[nev*n_l];
347  DGEMV(&cN, &n, &kdim, &done, Z, &n, y, &one, &dzero, u, &one);
348  /*-------------------- normalize u */
349  if (ifGenEv) {
350  /* B-norm, w2 = B*u */
351  matvec_B(u, w2);
352  t = sqrt(DDOT(&n, u, &one, w2, &one)); /* should be one */
353  } else {
354  /* 2-norm */
355  t = DNRM2(&n, u, &one); /* should be one */
356  }
357  /*-------------------- return code 2 --> zero eigenvector found */
358  if (t == 0.0) {
359  return 2;
360  }
361  /*-------------------- scal u */
362  t = 1.0 / t;
363  DSCAL(&n, &t, u, &one);
364  /*-------------------- scal B*u */
365  if (ifGenEv) {
366  /*------------------ w2 = B*u */
367  DSCAL(&n, &t, w2, &one);
368  }
369  /*-------------------- w = A*u */
370  matvec_A(u, wk);
371  /*-------------------- Ritz val: t = (u'*w)/(u'*u)
372  t = (u'*w)/(u'*B*u) */
373  t = DDOT(&n, wk, &one, u, &one);
374  /*-------------------- if lambda (==t) is in [a,b] */
375  if (t < aa - DBL_EPSILON || t > bb + DBL_EPSILON) {
376  continue;
377  }
378  /*-------------------- compute residual wrt A for this pair */
379  nt = -t;
380  if (ifGenEv) {
381  /*-------------------- w = w - t*B*u */
382  DAXPY(&n, &nt, w2, &one, wk, &one);
383  } else {
384  /*-------------------- w = w - t*u */
385  DAXPY(&n, &nt, u, &one, wk, &one);
386  }
387  /*-------------------- if diag scaling is present */
388  if (evsldata.ds) {
389  /* scale the residual */
390  int j;
391  for (j=0; j<n; j++) {
392  /* NOTE that we saved reciprocal of the diag scaling */
393  wk[j] /= evsldata.ds[j];
394  }
395  }
396  /*-------------------- res0 = 2-norm(wk) */
397  res0 = DNRM2(&n, wk, &one);
398  /*-------------------- accept (t, y) */
399  if (res0 < tol) {
400  Lam[nev] = t;
401  res[nev] = res0;
402  nev++;
403  }
404  }
405  tm1 = evsl_timer() - tt;
406 
407  /*-------------------- Done. output : */
408  *nevOut = nev;
409  *lamo = Lam;
410  *Wo = Rvec;
411  *reso = res;
412  /*-------------------- free arrays */
413  free(V);
414  free(dT);
415  free(eT);
416  free(EvalT);
417  free(EvecT);
418  free(wk);
419  if (vrand) {
420  free(vrand);
421  }
422  if (ifGenEv) {
423  free(Z);
424  }
425  /*-------------------- record stats */
426  tall = evsl_timer() - tall;
427 
428  evslstat.t_iter = tall;
429  evslstat.t_ritz = tm1;
430 
431  return 0;
432 }
433 
defs in EVSL
int ifGenEv
Definition: struct.h:166
#define orthTol
Definition: def.h:17
double DDOT(int *n, double *x, int *incx, double *y, int *incy)
void rand_double(int n, double *v)
Definition: vect.c:11
void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, double *w)
Classical GS reortho with Daniel, Gragg, Kaufman, Stewart test.
Definition: misc_la.c:195
#define NGS_MAX
Definition: misc_la.c:260
int n
Definition: struct.h:165
void DSCAL(int *n, double *a, double *x, int *incx)
int SymmTridEig(double *eigVal, double *eigVec, int n, const double *diag, const double *sdiag)
compute all eigenvalues and eigenvectors of a symmetric tridiagonal matrix
Definition: misc_la.c:36
This file contains function prototypes and constant definitions internally used in EVSL...
#define min(a, b)
Definition: def.h:62
#define Malloc(base, nmem, type)
Definition: def.h:22
void DGEMV(char *trans, int *m, int *n, double *alpha, double *a, int *lda, double *x, int *incx, double *beta, double *y, int *incy)
evslStat evslstat
global statistics of EVSL
Definition: evsl.c:21
int RatLanNr(double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
Rational filtering Lanczos process [NON-restarted version].
Definition: ratlanNr.c:54
void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q, double *v, double *w)
Classical GS reortho. No test. just do i_max times used in generalized ev problems.
Definition: misc_la.c:235
double t_ritz
Definition: struct.h:210
structs used in evsl
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
double DNRM2(int *n, double *x, int *incx)
Defs for blaslapack routines.
double * ds
Definition: struct.h:171
int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu, int *nevO, const double *diag, const double *sdiag)
compute eigenvalues and eigenvectors of a symmetric tridiagonal matrix in a slice ...
Definition: misc_la.c:93
evslData evsldata
global variable of EVSL
Definition: evsl.c:15
#define Realloc(base, nmem, type)
Definition: def.h:42
void RatFiltApply(int n, ratparams *rat, double *b, double *x, double *w6)
Apply rational filter R to a vetor b.
Definition: ratfilter.c:448
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
double t_iter
Definition: struct.h:200
void DCOPY(int *n, double *dx, int *incx, double *dy, int *incy)
parameters for rational filter
Definition: struct.h:112