EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
cheblanTr.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 "def.h"
7 #include "blaslapack.h"
8 #include "struct.h"
9 #include "internal_proto.h"
10 
18 #define FILTER_VINIT 1
19 
62 int ChebLanTr(int lanm, int nev, double *intv, int maxit,
63  double tol, double *vinit, polparams *pol, int *nev2,
64  double **vals, double **W, double **resW, FILE *fstats) {
65  const int ifGenEv = evsldata.ifGenEv;
66  /*-------------------- for stats */
67  double tall=0.0;
68  //double tolP = tol;
69  double tr, last_tr;
70  tall = evsl_timer();
71  int do_print = 1;
72  /* handle case where fstats is NULL. Then no output. Needed for openMP. */
73  if (fstats == NULL) {
74  do_print = 0;
75  }
76  /* size of the matrix */
77  int n = evsldata.n;
78  size_t n_l = n;
79  /*--------------------- adjust lanm and maxit */
80  lanm = min(lanm, n);
81  int lanm1=lanm+1;
82  size_t lanm1_l = lanm1;
83  /* if use full lanczos, should not do more than n iterations */
84  if (lanm == n) {
85  maxit = min(maxit, n);
86  }
87  /*-------------------- this is needed to increment the space when we
88  discover more than nev eigenvalues in interval */
89  double nevInc = 0.2; /* add 1 + 20% each time it is needed */
90  /*-------------------- if we have at least nev/ev_frac good candidate
91  eigenvalues from p(A) == then we restart to lock them in */
92  //int evFrac = 2;
93  /*-------------------- some constants frequently used */
94  /* char cT='T'; */
95  char cN = 'N';
96  int one = 1;
97  double done=1.0, dmone=-1.0, dzero=0.0;
98  /*-------------------- Ntest = when to start testing convergence */
99  int Ntest = min(lanm, nev+50);
100  /*-------------------- how often to test */
101  int cycle = 50;
102  int i, ll, /* count, last_count,*/ jl, last_jl;
103  /*-----------------------------------------------------------------------
104  -----------------------------------------------------------------------*/
105  /* check if the given interval is valid */
106  if (check_intv(intv, fstats) < 0) {
107  *nev2 = 0;
108  *vals = NULL; *W = NULL; *resW = NULL;
109  return 0;
110  }
111  double aa = intv[0];
112  double bb = intv[1];
113  /*-------------------- Assumption: polynomial pol computed before calling cheblanTr
114  pol. approximates the delta function centered at 'gamB'
115  bar: a bar value to threshold Ritz values of p(A) */
116  //int deg = pol->deg;
117  double gamB=pol->gam, bar=pol->bar;
118  /*-------------------- gamB must be within [-1, 1] */
119  if (gamB > 1.0 || gamB < -1.0) {
120  if (fstats) {
121  fprintf(fstats, "gamB error %.15e\n", gamB);
122  }
123  return -1;
124  }
125  /*-----------------------------------------------------------------------*
126  * *thick restarted* Lanczos step
127  *-----------------------------------------------------------------------*/
128  if (do_print) {
129  fprintf(fstats, " Cheb-LanTR, dim %d cycle %d \n",lanm, cycle);
130  }
131  /*--------------------- the min number of steps to be performed for
132  * each innter loop, must be >= 1 - if not,
133  * it means that the Krylov dim is too small */
134  int min_inner_step = 5;
135  /*-------------------- it = number of Lanczos steps */
136  int it = 0;
137  /*-------------------- Lanczos vectors V_m and tridiagonal matrix T_m */
138  double *V, *T;
139  Malloc(V, n_l*lanm1, double);
140  /*-------------------- for gen eig prob, storage for Z = B * V */
141  double *Z;
142  if (ifGenEv) {
143  Malloc(Z, n_l*lanm1, double);
144  } else {
145  Z = V;
146  }
147  /*-------------------- T must be zeroed out initially */
148  Calloc(T, lanm1_l*lanm1_l, double);
149  /*-------------------- Lam, Y: the converged (locked) Ritz values/vectors
150  res: related residual norms */
151  double *Y, *Lam, *res;
152  Malloc(Y, n_l*nev, double);
153  Malloc(Lam, nev, double);
154  Malloc(res, nev, double);
155  double *BY = NULL;
156  /*-------------------- for gen eig prob, storage for B*Y */
157  if (ifGenEv) {
158  Malloc(BY, n_l*nev, double);
159  }
160  /*-------------------- lock = number of locked vectors */
161  int lock = 0;
162  /*-------------------- trlen = dim. of thick restart set */
163  int trlen = 0, prtrlen=-1;
164  /*-------------------- Ritz values and vectors of p(A) */
165  double *Rval, *Rvec, *resi, *BRvec = NULL;
166  Malloc(Rval, lanm, double);
167  Malloc(resi, lanm, double);
168  Malloc(Rvec, n_l*lanm, double);
169  if (ifGenEv) {
170  Malloc(BRvec, n_l*lanm, double);
171  }
172  /*-------------------- Eigen vectors of T */
173  double *EvecT;
174  Malloc(EvecT, lanm1_l*lanm1_l, double);
175  /*-------------------- s used by TR (the ``spike'' of 1st block in Tm)*/
176  double *s;
177  Malloc(s, lanm, double);
178  /*-------------------- alloc some work space */
179  double *work, *vrand = NULL;
180  size_t work_size = ifGenEv ? 4*n_l : 3*n_l;
181  Malloc(work, work_size, double);
182 #if FILTER_VINIT
183  /*------------------ Filter the initial vector*/
184  ChebAv(pol, vinit, V, work);
185  Malloc(vrand, n, double);
186 #else
187  /*-------------------- copy initial vector to V(:,1) */
188  DCOPY(&n, vinit, &one, V, &one);
189 #endif
190  /*-------------------- normalize it */
191  double t;
192  if (ifGenEv) {
193  /* B norm */
194  matvec_B(V, Z);
195  t = 1.0 / sqrt(DDOT(&n, V, &one, Z, &one));
196  /* z = B*v */
197  DSCAL(&n, &t, Z, &one);
198  } else {
199  /* 2-norm */
200  t = 1.0 / DNRM2(&n, V, &one);
201  }
202  /* unit B-norm or 2-norm */
203  DSCAL(&n, &t, V, &one);
204  /*-------------------- main (restarted Lan) outer loop */
205  while (it < maxit) {
206  /*-------------------- for ortho test */
207  double wn = 0.0;
208  int nwn = 0;
209  /* beta */
210  double beta = 0.0, r, res0;
211  /* start with V(:,k) */
212  int k = trlen > 0 ? trlen + 1 : 0;
213  /* ! add a test if dimension exceeds (m+1)
214  * (trlen + 1) + min_inner_step <= lanm + 1 */
215  if (k+min_inner_step > lanm1) {
216  if (fstats) {
217  fprintf(fstats, "Krylov dim too small for this problem. Try a larger dim\n");
218  }
219  exit(1);
220  }
221  /*-------------------- thick restart special step */
222  if (trlen > 0) {
223  int k1 = k-1;
224  /*------------------ a quick reference to V(:,k) */
225  double *v = &V[k1*n_l];
226  double *z = &Z[k1*n_l];
227  /*------------------ next Lanczos vector */
228  double *vnew = v + n;
229  double *znew = z + n;
230  /*------------------ znew = p[(A-cc)/dd] * v */
231  /*------------------ NOTE: z is used!!! [TODO: FIX ME] */
232  ChebAv(pol, z, znew, work);
233  /*------------------ deflation */
234  if (lock > 0) {
235  if (ifGenEv) {
236  /* orthogonalize against locked vectors first, w = w - B*Y*Y'*w */
237  CGS_DGKS2(n, lock, NGS_MAX, BY, Y, znew, work);
238  } else {
239  /* orthogonalize against locked vectors first, w = w - Y*Y'*w */
240  CGS_DGKS(n, lock, NGS_MAX, Y, vnew, NULL, work);
241  }
242  }
243  /*-------------------- restart with 'trlen' Ritz values/vectors
244  T = diag(Rval(1:trlen)) */
245  for (i=0; i<trlen; i++) {
246  T[i*lanm1_l+i] = Rval[i];
247  wn += fabs(Rval[i]);
248  }
249  /*--------------------- s(k) = V(:,k)'* znew */
250  s[k1] = DDOT(&n, v, &one, znew, &one);
251  /*--------------------- znew = znew - Z(:,1:k)*s(1:k) */
252  DGEMV(&cN, &n, &k, &dmone, Z, &n, s, &one, &done, znew, &one);
253  /*-------------------- expand T matrix to k-by-k, arrow-head shape
254  T = [T, s(1:k-1)] then T = [T; s(1:k)'] */
255  for (i=0; i<k1; i++) {
256  T[trlen*lanm1_l+i] = s[i];
257  T[i*lanm1_l+trlen] = s[i];
258  wn += 2.0 * fabs(s[i]);
259  }
260  T[trlen*lanm1_l+trlen] = s[k1];
261  wn += fabs(s[k1]);
262  if (ifGenEv) {
263  /*-------------------- vnew = B \ znew */
264  solve_B(znew, vnew);
265  /*-------------------- beta = (vnew, znew)^{1/2} */
266  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
267  } else {
268  /*-------------------- beta = norm(w) */
269  beta = DNRM2(&n, vnew, &one);
270  }
271  wn += 2.0 * beta;
272  nwn += 3*k;
273  /* beta ~ 0 */
274  if (beta*nwn < orthTol*wn) {
275  if (do_print) {
276  fprintf(fstats, "it %4d: Lucky breakdown, beta = %.15e\n", k, beta);
277  }
278 #if FILTER_VINIT
279  /* filter random vector */
280  rand_double(n, vrand);
281  ChebAv(pol, vrand, vnew, work);
282 #else
283  rand_double(n, vnew);
284 #endif
285  if (ifGenEv) {
286  /* orthogonalize against locked vectors first, v = v - Y*(B*Y)'*v */
287  CGS_DGKS2(n, lock, NGS_MAX, Y, BY, vnew, work);
288  /* vnew = vnew - V(:,1:k)*Z(:,1:k)'*vnew */
289  CGS_DGKS2(n, k, NGS_MAX, V, Z, vnew, work);
290  matvec_B(vnew, znew);
291  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
292  double ibeta = 1.0 / beta;
293  DSCAL(&n, &ibeta, vnew, &one);
294  DSCAL(&n, &ibeta, znew, &one);
295  beta = 0.0;
296  } else {
297  /* orthogonalize against locked vectors first, w = w - Y*Y'*w */
298  CGS_DGKS(n, lock, NGS_MAX, Y, vnew, NULL, work);
299  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
300  /* beta = norm(w) */
301  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
302  double ibeta = 1.0 / beta;
303  DSCAL(&n, &ibeta, vnew, &one);
304  beta = 0.0;
305  }
306  } else {
307  /*------------------- w = w / beta */
308  double ibeta = 1.0 / beta;
309  DSCAL(&n, &ibeta, vnew, &one);
310  if (ifGenEv) {
311  DSCAL(&n, &ibeta, znew, &one);
312  }
313  }
314  /*------------------- T(k+1,k) = beta; T(k,k+1) = beta; */
315  T[k1*lanm1_l+k] = beta;
316  T[k*lanm1_l+k1] = beta;
317  } /* if (trlen > 0) */
318  /*-------------------- Done with TR step. Rest of Lanczos step */
319  /*-------------------- reset Ntest at each restart. */
320  Ntest = max(20,nev-lock+10);
321  /*last_count = 0;*/ last_jl = 0; last_tr = 0.0;
322  /*-------------------- regardless of trlen, *(k+1)* is the current
323  * number of Lanczos vectors in V */
324  /*-------------------- pointer to the previous Lanczos vector */
325  double *zold = k > 0 ? Z+(k-1)*n_l : NULL;
326  /*------------------------------------------------------*/
327  /*------------------ Lanczos inner loop ----------------*/
328  /*------------------------------------------------------*/
329  while (k < lanm && it < maxit) {
330  k++;
331  /*---------------- a quick reference to V(:,k) */
332  double *v = &V[(k-1)*n_l];
333  double *z = &Z[(k-1)*n_l];
334  /*---------------- next Lanczos vector */
335  double *vnew = v + n;
336  double *znew = z + n;
337  /*------------------ znew = p[(A-cc)/dd] * v */
338  /*------------------ NOTE: z is used!!! [TODO: FIX ME] */
339  ChebAv(pol, z, znew, work);
340  it++;
341  /*-------------------- deflation: orthgonalize vs locked ones first */
342  if (lock > 0) {
343  if (ifGenEv) {
344  /* orthgonlize against locked vectors first, znew = znew - B*Y*Y'*znew */
345  CGS_DGKS2(n, lock, NGS_MAX, BY, Y, znew, work);
346  } else {
347  /*-------------------- vnew = vnew - Y*Y'*vnew */
348  CGS_DGKS(n, lock, NGS_MAX, Y, vnew, NULL, work);
349  }
350  }
351  /*-------------------- znew = znew - beta*zold */
352  if (zold) {
353  double nbeta = -beta;
354  DAXPY(&n, &nbeta, zold, &one, znew, &one);
355  }
356  /*-------------------- alpha = znew'*v */
357  double alpha = DDOT(&n, v, &one, znew, &one);
358  /*-------------------- T(k,k) = alpha */
359  T[(k-1)*lanm1_l+(k-1)] = alpha;
360  wn += fabs(alpha);
361  /*-------------------- znew = znew - alpha*z */
362  double nalpha = -alpha;
363  DAXPY(&n, &nalpha, z, &one, znew, &one);
364  /*-------------------- FULL reortho to all previous Lan vectors */
365  if (ifGenEv) {
366  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
367  CGS_DGKS2(n, k, NGS_MAX, Z, V, znew, work);
368  /*-------------------- vnew = B \ znew */
369  solve_B(znew, vnew);
370  /*-------------------- beta = (vnew, znew)^{1/2} */
371  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
372  } else {
373  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
374  /* beta = norm(w) */
375  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
376  }
377  wn += 2.0 * beta;
378  nwn += 3;
379  /*-------------------- zold = z */
380  zold = z;
381  /*-------------------- lucky breakdown test */
382  if (beta*nwn < orthTol*wn) {
383  if (do_print) {
384  fprintf(fstats, "it %4d: Lucky breakdown, beta = %.15e\n", it, beta);
385  }
386 #if FILTER_VINIT
387  /* filter random vector */
388  rand_double(n, vrand);
389  ChebAv(pol, vrand, vnew, work);
390 #else
391  /* generate a new init vector*/
392  rand_double(n, vnew);
393 #endif
394  if (ifGenEv) {
395  /* orthogonalize against locked vectors first, w = w - Y*(B*Y)'*w */
396  CGS_DGKS2(n, lock, NGS_MAX, Y, BY, vnew, work);
397  /* vnew = vnew - V(:,1:k)*Z(:,1:k)'*vnew */
398  CGS_DGKS2(n, k, NGS_MAX, V, Z, vnew, work);
399  matvec_B(vnew, znew);
400  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
401  double ibeta = 1.0 / beta;
402  DSCAL(&n, &ibeta, vnew, &one);
403  DSCAL(&n, &ibeta, znew, &one);
404  beta = 0.0;
405  } else {
406  /* orthogonalize against locked vectors first, w = w - Y*Y'*w */
407  CGS_DGKS(n, lock, NGS_MAX, Y, vnew, NULL, work);
408  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
409  /* beta = norm(w) */
410  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
411  double ibeta = 1.0 / beta;
412  DSCAL(&n, &ibeta, vnew, &one);
413  beta = 0.0;
414  }
415  } else {
416  /*---------------------- vnew = vnew / beta */
417  double ibeta = 1.0 / beta;
418  DSCAL(&n, &ibeta, vnew, &one);
419  if (ifGenEv) {
420  /*-------------------- znew = znew / beta */
421  DSCAL(&n, &ibeta, znew, &one);
422  }
423  }
424  /*-------------------- T(k,k+1) = T(k+1,k) = beta */
425  T[k*lanm1_l+(k-1)] = beta;
426  T[(k-1)*lanm1_l+k] = beta;
427  /*---------------------- Restarting test */
428  int k1 = k-trlen-Ntest;
429  if ( ((k1>=0) && (k1 % cycle == 0)) || (k == lanm) || it == maxit) {
430  /*-------------------- solve eigen-problem for T(1:k,1:k)
431  vals in Rval, vecs in EvecT */
432  /* printf("k %d, trlen %d, Ntest %d, its %d\n", k, trlen, Ntest, it); */
433  SymEigenSolver(k, T, lanm1, EvecT, lanm1, Rval);
434  /*-------------------- max dim reached-break from the inner loop */
435  if (k == lanm || it == maxit) {
436  break;
437  }
438 #if 1
439  /* we change the convergence test to be simple:
440  * we test if the sum and the number of the Ritz values that are >= bar no longer change */
441  jl = 0;
442  tr = 0.0;
443  for (i=0; i<k; i++) {
444  if (Rval[i] + DBL_EPSILON >= bar) {
445  jl++;
446  tr += Rval[i];
447  }
448  }
449  if (fabs(tr-last_tr) <= 2e-12 && jl == last_jl) {
450  if (fstats) {
451  fprintf(fstats,"break: [it %d, k %d]: last_tr %.15e, tr %.15e, last_jl %d jl %d\n",
452  it, k, last_tr, tr, last_jl, jl);
453  }
454  break;
455  }
456  if (fstats) {
457  fprintf(fstats,"[it %d, k %d] testing: last_tr %.15e, tr %.15e, last_jl %d jl %d\n",
458  it, k, last_tr, tr, last_jl, jl);
459  }
460  last_tr = tr;
461  last_jl = jl;
462 #else
463  count = 0;
464  /*-------------------- get residual norms and check acceptance of
465  Ritz values for p(A). */
466  /*-------------------- count e.vals in interval + those that convergeed */
467  jl = 0;
468  for (i=0; i<k; i++) {
469  if (Rval[i]>= bar) {
470  jl++;
471  /* for standard e.v prob, this is the 2-norm of A*y-lam*y
472  * for gen e.v prob, this is the B-norm of B^{-1}*A*y-lam*y */
473  r = fabs(beta*EvecT[i*lanm1+(k-1)]);
474  resi[i] = r;
475  if (r < tolP) {
476  count++;
477  }
478  }
479  }
480  /*-------------------- testing (partial) convergence for restart */
481  if (do_print) {
482  fprintf(fstats,"inner: testing conv k %4d, it %4d, count %3d jl %3d trlen %3d\n",
483  k, it, count, jl, trlen);
484  }
485  /*-------------------- enough good candidates 1st time -> break */
486  /* if ((count*evFrac >= nev-lock) && (prtrlen==-1)) */
487  if (count*evFrac >= nev-lock) {
488  fprintf(fstats, "count * evFrac >= nev-lock: %d * %d >= %d - %d\n", count, evFrac, nev, lock);
489  break;
490  }
491  /*-------------------- count & jl unchanged since last test --> break */
492  if ((count<=last_count) && (jl<=last_jl)) {
493  fprintf(fstats, "inner: count <= last_count && jl <= last_jl: %d <= %d && %d <= %d\n", count, last_count, jl, last_jl);
494  break;
495  }
496  last_count = count;
497  last_jl = jl;
498 #endif
499  } /* if [Restarting test] block */
500  /*-------------------- end of inner (Lanczos) loop - Next: restart*/
501  } /* while (k<mlan) loop */
502 
503 
504  //SymEigenSolver(k, T, lanm1, EvecT, lanm1, Rval);
505  //savedensemat(T, lanm1, k, k, "T.txt");
506  //savedensemat(EvecT, lanm1, k, k, "Evec.txt");
507  //save_vec(k, Rval, "eval.txt");
508  /*-------------------- TWO passes to select good candidates */
509  /* Pass-1: based on if ``p(Ritzvalue) > bar'' */
510  jl = 0;
511  for (i=0; i<k; i++) {
512  //printf("resi[%d] = %.15e\n", i, fabs(beta*EvecT[i*lanm1+(k-1)]));
513  /*-------------------- if this Ritz value is higher than ``bar'' */
514  if (Rval[i] >= bar) {
515  /* move good eigenvectors/vals to front */
516  if (i != jl) {
517  DCOPY(&k, EvecT+i*lanm1_l, &one, EvecT+jl*lanm1_l, &one);
518  Rval[jl] = Rval[i];
519  //resi[jl] = resi[i];
520  }
521  resi[jl] = fabs(beta*EvecT[i*lanm1_l+(k-1)]);
522  //printf("beta = %.15e, resi[%d] = %.15e\n", beta, i, resi[jl]);
523  jl++;
524  }
525  }
526  //exit(0);
527  //fprintf(fstats, "beta = %.1e\n", beta);
528  /*---------------------- Compute the Ritz vectors:
529  * Rvec(:,1:jl) = V(:,1:k) * EvecT(:,1:jl) */
530  DGEMM(&cN, &cN, &n, &jl, &k, &done, V, &n, EvecT, &lanm1, &dzero, Rvec, &n);
531  if (ifGenEv) {
532  DGEMM(&cN, &cN, &n, &jl, &k, &done, Z, &n, EvecT, &lanm1, &dzero, BRvec, &n);
533  }
534  /*-------------------- Pass-2: check if Ritz vals of A are in [a,b] */
535  /* number of Ritz values in [a,b] */
536  ll = 0;
537  /*-------------------- trlen = # Ritz vals that will go to TR set */
538  prtrlen = trlen;
539  trlen = 0;
540  for (i=0; i<jl; i++) {
541  double *y = Rvec + i*n_l;
542  double *By = NULL;
543  if (ifGenEv) {
544  By = BRvec + i*n_l;
545  }
546  double *w = work;
547  double *w2 = w + n;
548  /*------------------ normalize just in case. */
549  if (ifGenEv) {
550  /* B-norm, w2 = B*y */
551  matvec_B(y, w2);
552  t = sqrt(DDOT(&n, y, &one, w2, &one));
553  } else {
554  /* 2-norm */
555  t = DNRM2(&n, y, &one);
556  }
557  /*-------------------- return code 2 --> zero eigenvector found */
558  if (t == 0.0) {
559  return 2;
560  }
561  /*-------------------- scal y */
562  t = 1.0 / t;
563  DSCAL(&n, &t, y, &one);
564  /*-------------------- scal B*y */
565  if (ifGenEv) {
566  DSCAL(&n, &t, w2, &one);
567  }
568  /*-------------------- w = A*y */
569  matvec_A(y, w);
570  /*-------------------- Ritzval: t3 = (y'*w)/(y'*y) or
571  * t3 = (y'*w)/(y'*B*y) */
572  /*-------------------- Rayleigh quotient */
573  double t3 = DDOT(&n, y, &one, w, &one);
574  /*-------------------- if lambda (==t3) is in [a,b] */
575  if (t3 >= aa - DBL_EPSILON && t3 <= bb + DBL_EPSILON) {
576  ll++;
577  /*-------------------- compute residual wrt A for this pair */
578  double nt3 = -t3;
579  if (ifGenEv) {
580  /* w = w - t3*w2, w2 = B*y, (w=A*y-t3*B*y) */
581  DAXPY(&n, &nt3, w2, &one, w, &one);
582  } else {
583  /*-------------------- w = w - t3*y, (w=A*y-t3*y) */
584  DAXPY(&n, &nt3, y, &one, w, &one);
585  }
586  /*-------------------- res0 = 2-norm of w */
587  res0 = DNRM2(&n, w, &one);
588  /*-------------------- test res. of this Ritz pair against tol */
589  /* r = resi[i];*/
590  r = res0;
591  if (r < tol) {
592  //-------------------- check if need to realloc
593  if (lock >= nev) {
594  nev += 1 + (int) (nev*nevInc);
595  if (do_print) {
596  fprintf(fstats, "-- More eigval found: realloc space for %d evs\n", nev);
597  }
598  Realloc(Y, nev*n_l, double);
599  if (ifGenEv) {
600  Realloc(BY, nev*n_l, double);
601  }
602  Realloc(Lam, nev, double);
603  Realloc(res, nev, double);
604  }
605  /*-------------------- accept (t3, y) */
606  DCOPY(&n, y, &one, Y+lock*n_l, &one);
607  if (ifGenEv) {
608  DCOPY(&n, By, &one, BY+lock*n_l, &one);
609  }
610  Lam[lock] = t3;
611  res[lock] = res0;
612  lock++;
613  } else {
614  /*-------------------- restart; move Ritz pair for TR to front */
615  Rval[trlen] = Rval[i];
616  DCOPY(&n, y, &one, V+trlen*n_l, &one);
617  if (ifGenEv) {
618  DCOPY(&n, By, &one, Z+trlen*n_l, &one);
619  }
620  /* special vector for TR that is the bottom row of
621  * eigenvectors of Tm */
622  s[trlen] = beta * EvecT[i*lanm1_l+(k-1)];
623  trlen ++;
624  }
625  }
626  }/* end of 2nd pass */
627 
628  /*-------Note: jl = #evs that passed the 1st test,
629  * ll = #evs that passed the 1st and the 2nd tests.
630  * These ll evs are either locked (accepted)
631  * or put into trlen (candidates for next restarts)
632  * step when trlen = 0 last restart and ll=0 this time.
633  * this is a sort of confirmation that nothing is left.
634  * another test may be added later to make it more rigorous.
635  */
636  if (do_print) {
637  fprintf(fstats,"it %4d: k %3d, jl %3d, ll %3d, lock %3d, trlen %3d\n",
638  it, k, jl, ll, lock, trlen);
639  }
640  /*-------------------- TESTs for stopping */
641  if ((prtrlen == 0) && (ll==0)) {
642  /*-------------------- It looks like: Nothing left to compute */
643  if (do_print) {
644  fprintf(fstats, "--------------------------------------------------\n");
645  fprintf(fstats, " --> No ev.s left to be computed\n");
646  }
647  break;
648  }
649  if (it == maxit) {
650  if (do_print) {
651  fprintf(fstats, "--------------------------------------------------\n");
652  fprintf(fstats, " --> Max its reached without convergence\n");
653  }
654  break;
655  }
656  /*-------------------- prepare to restart. First zero out all T */
657  memset(T, 0, lanm1_l*lanm1_l*sizeof(double));
658  /*-------------------- move starting vector vector V(:,k+1); V(:,trlen+1) = V(:,k+1) */
659  DCOPY(&n, V+k*n_l, &one, V+trlen*n_l, &one);
660  if (ifGenEv) {
661  DCOPY(&n, Z+k*n_l, &one, Z+trlen*n_l, &one);
662  }
663  } /* outer loop (it) */
664 
665  if (do_print) {
666  fprintf(fstats, " Number of evals found = %d\n", lock);
667  fprintf(fstats, "--------------------------------------------------\n");
668  }
669 
670  /*-------------------- Done. output : */
671  *nev2 = lock;
672  *vals = Lam;
673  *W = Y;
674  *resW = res;
675  /*-------------------- free arrays */
676  free(V);
677  free(T);
678  free(Rval);
679  free(resi);
680  free(EvecT);
681  free(Rvec);
682  free(s);
683  free(work);
684  if (vrand) {
685  free(vrand);
686  }
687  if (ifGenEv) {
688  free(Z);
689  free(BY);
690  free(BRvec);
691  }
692  /*-------------------- record stats */
693  tall = evsl_timer() - tall;
694  /*-------------------- print stat */
695  evslstat.t_iter = tall;
696 
697  return 0;
698 }
699 
void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam)
interface to LAPACK SYMMETRIC EIGEN-SOLVER
Definition: misc_la.c:156
defs in EVSL
int ifGenEv
Definition: struct.h:166
int ChebLanTr(int lanm, int nev, double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats)
Chebyshev polynomial filtering Lanczos process [Thick restart version].
Definition: cheblanTr.c:62
void DGEMM(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda, double *b, int *ldb, double *beta, double *c, int *ldc)
#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
double gam
Definition: struct.h:64
#define NGS_MAX
Definition: misc_la.c:260
int n
Definition: struct.h:165
void DSCAL(int *n, double *a, double *x, int *incx)
#define Calloc(base, nmem, type)
Definition: def.h:32
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
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
structs used in evsl
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
#define max(a, b)
Definition: def.h:56
double DNRM2(int *n, double *x, int *incx)
Defs for blaslapack routines.
int ChebAv(polparams *pol, double *v, double *y, double *w)
Computes y=P(A) v, where pn is a Cheb. polynomial expansion
Definition: chebpoly.c:510
evslData evsldata
global variable of EVSL
Definition: evsl.c:15
parameters for polynomial filter
Definition: struct.h:45
#define Realloc(base, nmem, type)
Definition: def.h:42
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
double bar
Definition: struct.h:65
double t_iter
Definition: struct.h:200
void DCOPY(int *n, double *dx, int *incx, double *dy, int *incy)