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