EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
lanTrbounds.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 
11 #define COMP_RES 0
12 
38 int LanTrbounds(int lanm, int maxit, double tol, double *vinit,
39  int bndtype,
40  double *lammin, double *lammax, FILE *fstats) {
41  const int ifGenEv = evsldata.ifGenEv;
42  double lmin=0.0, lmax=0.0, t, t1, t2;
43  int do_print = 1;
44  /* handle case where fstats is NULL. Then no output. Needed for openMP. */
45  if (fstats == NULL) {
46  do_print = 0;
47  }
48  /* size of the matrix */
49  int n = evsldata.n;
50  size_t n_l = n;
51  /*--------------------- adjust lanm and maxit */
52  lanm = min(lanm, n);
53  int lanm1=lanm+1;
54  /* if use full lanczos, should not do more than n iterations */
55  if (lanm == n) {
56  maxit = min(maxit, n);
57  }
58  size_t lanm1_l = lanm1;
59  /*-------------------- some constants frequently used */
60  /* char cT='T'; */
61  char cN = 'N';
62  int one = 1;
63  double done=1.0, dmone=-1.0, dzero=0.0;
64  int i;
65  /*-----------------------------------------------------------------------*
66  * *thick restarted* Lanczos step
67  *-----------------------------------------------------------------------*/
68  if (do_print) {
69  fprintf(fstats, " LanTR for bounds: dim %d, maxits %d\n", lanm, maxit);
70  }
71  /*--------------------- the min number of steps to be performed for
72  * each innter loop, must be >= 1 - if not,
73  * it means that the Krylov dim is too small */
74  int min_inner_step = 5;
75  /*-------------------- it = number of Lanczos steps */
76  int it = 0;
77  /*-------------------- Lanczos vectors V_m and tridiagonal matrix T_m */
78  double *V, *T;
79  Malloc(V, n_l*lanm1, double);
80  /*-------------------- for gen eig prob, storage for Z = B * V */
81  double *Z;
82  if (ifGenEv) {
83  Malloc(Z, n_l*lanm1, double);
84  } else {
85  Z = V;
86  }
87  /*-------------------- T must be zeroed out initially */
88  Calloc(T, lanm1_l*lanm1_l, double);
89  /*-------------------- trlen = dim. of thick restart set */
90  int trlen = 0;
91  /*-------------------- Ritz values and vectors of p(A) */
92  double *Rval, *Rvec, *BRvec=NULL;
93  Malloc(Rval, lanm, double);
94  /*-------------------- Only compute 2 Ritz vectors */
95  Malloc(Rvec, n_l*2, double);
96  if (ifGenEv) {
97  Malloc(BRvec, n_l*2, double);
98  }
99  /*-------------------- Eigen vectors of T */
100  double *EvecT;
101  Malloc(EvecT, lanm1_l*lanm1_l, double);
102  /*-------------------- s used by TR (the ``spike'' of 1st block in Tm)*/
103  double s[3];
104  /*-------------------- alloc some work space */
105  double *work;
106  size_t work_size = 2*n_l;
107  Malloc(work, work_size, double);
108  /*-------------------- copy initial vector to V(:,1) */
109  DCOPY(&n, vinit, &one, V, &one);
110  /*-------------------- normalize it */
111  if (ifGenEv) {
112  /* B norm */
113  matvec_B(V, Z);
114  t = 1.0 / sqrt(DDOT(&n, V, &one, Z, &one));
115  /* z = B*v */
116  DSCAL(&n, &t, Z, &one);
117  } else {
118  /* 2-norm */
119  t = 1.0 / DNRM2(&n, V, &one);
120  }
121  /* unit B-norm or 2-norm */
122  DSCAL(&n, &t, V, &one);
123  /*-------------------- main (restarted Lan) outer loop */
124  while (it < maxit) {
125  /*-------------------- for ortho test */
126  double wn = 0.0;
127  int nwn = 0;
128  /* beta */
129  double beta = 0.0;
130  /* start with V(:,k) */
131  int k = trlen > 0 ? trlen + 1 : 0;
132  /* ! add a test if dimension exceeds (m+1)
133  * (trlen + 1) + min_inner_step <= lanm + 1 */
134  if (k+min_inner_step > lanm1) {
135  fprintf(stderr, "Krylov dim too small for this problem. Try a larger dim\n");
136  exit(1);
137  }
138  /*-------------------- thick restart special step */
139  if (trlen > 0) {
140  int k1 = k-1;
141  /*------------------ a quick reference to V(:,k) */
142  double *v = &V[k1*n_l];
143  double *z = &Z[k1*n_l];
144  /*------------------ next Lanczos vector */
145  double *vnew = v + n;
146  double *znew = z + n;
147  /*------------------ znew = A * v */
148  matvec_A(v, znew);
149  /*-------------------- restart with 'trlen' Ritz values/vectors
150  T = diag(Rval(1:trlen)) */
151  for (i=0; i<trlen; i++) {
152  T[i*lanm1_l+i] = Rval[i];
153  wn += fabs(Rval[i]);
154  }
155  /*--------------------- s(k) = V(:,k)'* znew */
156  s[k1] = DDOT(&n, v, &one, znew, &one);
157  /*--------------------- znew = znew - Z(:,1:k)*s(1:k) */
158  DGEMV(&cN, &n, &k, &dmone, Z, &n, s, &one, &done, znew, &one);
159  /*-------------------- expand T matrix to k-by-k, arrow-head shape
160  T = [T, s(1:k-1)] then T = [T; s(1:k)'] */
161  for (i=0; i<k1; i++) {
162  T[trlen*lanm1_l+i] = s[i];
163  T[i*lanm1_l+trlen] = s[i];
164  wn += 2.0 * fabs(s[i]);
165  }
166  T[trlen*lanm1_l+trlen] = s[k1];
167  wn += fabs(s[k1]);
168  if (ifGenEv) {
169  /*-------------------- vnew = B \ znew */
170  solve_B(znew, vnew);
171  /*-------------------- beta = (vnew, znew)^{1/2} */
172  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
173  } else {
174  /*-------------------- beta = norm(w) */
175  beta = DNRM2(&n, vnew, &one);
176  }
177  wn += 2.0 * beta;
178  nwn += 3*k;
179  /* beta ~ 0 */
180  if (beta*nwn < orthTol*wn) {
181  rand_double(n, vnew);
182  if (ifGenEv) {
183  /* vnew = vnew - V(:,1:k)*Z(:,1:k)'*vnew */
184  CGS_DGKS2(n, k, NGS_MAX, V, Z, vnew, work);
185  matvec_B(vnew, znew);
186  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
187  double ibeta = 1.0 / beta;
188  DSCAL(&n, &ibeta, vnew, &one);
189  DSCAL(&n, &ibeta, znew, &one);
190  beta = 0.0;
191  } else {
192  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
193  /* beta = norm(w) */
194  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
195  double ibeta = 1.0 / beta;
196  DSCAL(&n, &ibeta, vnew, &one);
197  beta = 0.0;
198  }
199  } else {
200  /*------------------- w = w / beta */
201  double ibeta = 1.0 / beta;
202  DSCAL(&n, &ibeta, vnew, &one);
203  if (ifGenEv) {
204  DSCAL(&n, &ibeta, znew, &one);
205  }
206  }
207  /*------------------- T(k+1,k) = beta; T(k,k+1) = beta; */
208  T[k1*lanm1_l+k] = beta;
209  T[k*lanm1_l+k1] = beta;
210  } /* if (trlen > 0) */
211  /*-------------------- Done with TR step. Rest of Lanczos step */
212  /*-------------------- regardless of trlen, *(k+1)* is the current
213  * number of Lanczos vectors in V */
214  /*-------------------- pointer to the previous Lanczos vector */
215  double *zold = k > 0 ? Z+(k-1)*n_l : NULL;
216  /*------------------------------------------------------*/
217  /*------------------ Lanczos inner loop ----------------*/
218  /*------------------------------------------------------*/
219  while (k < lanm && it < maxit) {
220  k++;
221  /*---------------- a quick reference to V(:,k) */
222  double *v = &V[(k-1)*n_l];
223  double *z = &Z[(k-1)*n_l];
224  /*---------------- next Lanczos vector */
225  double *vnew = v + n;
226  double *znew = z + n;
227  /*------------------ znew = A * v */
228  matvec_A(v, znew);
229  it++;
230  /*-------------------- znew = znew - beta*zold */
231  if (zold) {
232  double nbeta = -beta;
233  DAXPY(&n, &nbeta, zold, &one, znew, &one);
234  }
235  /*-------------------- alpha = znew'*v */
236  double alpha = DDOT(&n, v, &one, znew, &one);
237  /*-------------------- T(k,k) = alpha */
238  T[(k-1)*lanm1_l+(k-1)] = alpha;
239  wn += fabs(alpha);
240  /*-------------------- znew = znew - alpha*z */
241  double nalpha = -alpha;
242  DAXPY(&n, &nalpha, z, &one, znew, &one);
243  /*-------------------- FULL reortho to all previous Lan vectors */
244  if (ifGenEv) {
245  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
246  CGS_DGKS2(n, k, NGS_MAX, Z, V, znew, work);
247  /* vnew = B \ znew */
248  solve_B(znew, vnew);
249  /*-------------------- beta = (vnew, znew)^{1/2} */
250  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
251  } else {
252  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
253  /* beta = norm(w) */
254  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
255  }
256  wn += 2.0 * beta;
257  nwn += 3;
258  /*-------------------- zold = z */
259  zold = z;
260  /*-------------------- lucky breakdown test */
261  if (beta*nwn < orthTol*wn) {
262  if (do_print) {
263  fprintf(fstats, "it %4d: Lucky breakdown, beta = %.15e\n", it, beta);
264  }
265  /* generate a new init vector*/
266  rand_double(n, vnew);
267  if (ifGenEv) {
268  /* vnew = vnew - V(:,1:k)*Z(:,1:k)'*vnew */
269  CGS_DGKS2(n, k, NGS_MAX, V, Z, vnew, work);
270  matvec_B(vnew, znew);
271  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
272  double ibeta = 1.0 / beta;
273  DSCAL(&n, &ibeta, vnew, &one);
274  DSCAL(&n, &ibeta, znew, &one);
275  beta = 0.0;
276  } else {
277  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
278  /* beta = norm(w) */
279  CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
280  double ibeta = 1.0 / beta;
281  DSCAL(&n, &ibeta, vnew, &one);
282  beta = 0.0;
283  }
284  } else {
285  /*---------------------- vnew = vnew / beta */
286  double ibeta = 1.0 / beta;
287  DSCAL(&n, &ibeta, vnew, &one);
288  if (ifGenEv) {
289  /*-------------------- znew = znew / beta */
290  DSCAL(&n, &ibeta, znew, &one);
291  }
292  }
293  /*-------------------- T(k,k+1) = T(k+1,k) = beta */
294  T[k*lanm1_l+(k-1)] = beta;
295  T[(k-1)*lanm1_l+k] = beta;
296  } /* while (k<mlan) loop */
297 
298  /*-------------------- solve eigen-problem for T(1:k,1:k)
299  vals in Rval, vecs in EvecT */
300  SymEigenSolver(k, T, lanm1, EvecT, lanm1, Rval);
301 
302  /*-------------------- Rval is in ascending order */
303  /*-------------------- Rval[0] is smallest, Rval[k-1] is largest */
304  /*-------------------- special vector for TR that is the bottom row of
305  eigenvectors of Tm */
306  s[0] = beta * EvecT[k-1];
307  s[1] = beta * EvecT[(k-1)*lanm1_l+(k-1)];
308  /*---------------------- bounds */
309  if (bndtype <= 1) {
310  /*-------------------- BOUNDS type 1 (simple) */
311  t1 = fabs(s[0]);
312  t2 = fabs(s[1]);
313  } else {
314  /*-------------------- BOUNDS type 2 (Kato-Temple) */
315  t1 = 2.0*s[0]*s[0] / (Rval[1] - Rval[0]);
316  t2 = 2.0*s[1]*s[1] / (Rval[k-1] - Rval[k-2]);
317  }
318  lmin = Rval[0] - t1;
319  lmax = Rval[k-1] + t2;
320  /*---------------------- Compute two Ritz vectors:
321  * Rvec(:,1) = V(:,1:k) * EvecT(:,1)
322  * Rvec(:,end) = V(:,1:k) * EvecT(:,end) */
323  DGEMV(&cN, &n, &k, &done, V, &n, EvecT, &one, &dzero, Rvec, &one);
324  DGEMV(&cN, &n, &k, &done, V, &n, EvecT+(k-1)*lanm1_l, &one, &dzero, Rvec+n, &one);
325  if (ifGenEv) {
326  DGEMV(&cN, &n, &k, &done, Z, &n, EvecT, &one, &dzero, BRvec, &one);
327  DGEMV(&cN, &n, &k, &done, Z, &n, EvecT+(k-1)*lanm1_l, &one, &dzero, BRvec+n, &one);
328  }
329  /*---------------------- Copy two Rval and Rvec to TR set */
330  trlen = 2;
331  for (i=0; i<2; i++) {
332  double *y = Rvec + i*n_l;
333  DCOPY(&n, y, &one, V+i*n_l, &one);
334  if (ifGenEv) {
335  double *By = BRvec + i*n_l;
336  DCOPY(&n, By, &one, Z+i*n_l, &one);
337  }
338  }
339  Rval[1] = Rval[k-1];
340  /*-------------------- recompute residual norm for debug only */
341 #if COMP_RES
342  /*-------------------- compute residual */
343  double *w1 = work;
344  double *w2 = work + n;
345  double rr[2];
346  for (i=0; i<2; i++) {
347  double *y = Rvec + i*n_l;
348  double nt = -Rval[i];
349  matvec_A(y, w1);
350  if (ifGenEv) {
351  matvec_B(y, w2);
352  DAXPY(&n, &nt, w2, &one, w1, &one);
353  solve_B(w1, w2);
354  rr[i] = sqrt(DDOT(&n, w1, &one, w2, &one));
355  } else {
356  DAXPY(&n, &nt, y, &one, w1, &one);
357  rr[i] = DNRM2(&n, w1, &one);
358  }
359  }
360  if (do_print) {
361  fprintf(fstats,"it %4d, k %3d: ritz %.15e %.15e, t1,t2 %e %e, res %.15e %.15e, comp res %.15e %.15e\n",
362  it, k, Rval[0], Rval[1], t1, t2, fabs(s[0]), fabs(s[1]), rr[0], rr[1]);
363  }
364 #else
365  if (do_print) {
366  fprintf(fstats,"it %4d, k %3d: ritz %.15e %.15e, t1,t2 %e %e, res %.15e %.15e\n",
367  it, k, Rval[0], Rval[1], t1, t2, fabs(s[0]), fabs(s[1]));
368  }
369 #endif
370  /*---------------------- test convergence */
371  if (t1+t2 < tol*(fabs(lmin)+fabs(lmax))) {
372  break;
373  }
374  /*-------------------- prepare to restart. First zero out all T */
375  memset(T, 0, lanm1_l*lanm1_l*sizeof(double));
376  /*-------------------- move starting vector vector V(:,k+1); V(:,trlen+1) = V(:,k+1) */
377  DCOPY(&n, V+k*n_l, &one, V+trlen*n_l, &one);
378  if (ifGenEv) {
379  DCOPY(&n, Z+k*n_l, &one, Z+trlen*n_l, &one);
380  }
381  } /* outer loop (it) */
382 
383  /*-------------------- Done. output : */
384  *lammin = lmin;
385  *lammax = lmax;
386  /*-------------------- free arrays */
387  free(V);
388  free(T);
389  free(Rval);
390  free(EvecT);
391  free(Rvec);
392  free(work);
393  if (ifGenEv) {
394  free(Z);
395  free(BRvec);
396  }
397 
398  return 0;
399 }
400 
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
#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
int LanTrbounds(int lanm, int maxit, double tol, double *vinit, int bndtype, double *lammin, double *lammax, FILE *fstats)
Lanczos process for eigenvalue bounds [Thick restart version].
Definition: lanTrbounds.c:38
#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)
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)
double DNRM2(int *n, double *x, int *incx)
Defs for blaslapack routines.
evslData evsldata
global variable of EVSL
Definition: evsl.c:15
void DCOPY(int *n, double *dx, int *incx, double *dy, int *incy)