EVSL  1.1.0 EigenValues Slicing Library
lanbounds.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include "def.h"
5 #include "blaslapack.h"
6 #include "struct.h"
7 #include "internal_proto.h"
8
21 int LanBounds(int msteps, double *v, double *lmin, double *lmax) {
22  const int ifGenEv = evsldata.ifGenEv;
23  double *alp, *bet, nbet, nalp, t, *V, *Z;
24  int j, one=1, n;
25  size_t n_l;
26
27  n = evsldata.n;
28  n_l = n;
29
30  msteps = min(n, msteps);
31  Malloc(alp, msteps, double);
32  Malloc(bet, msteps, double);
33  Malloc(V, (msteps+1)*n_l, double);
34  if (ifGenEv) {
35  /* storage for Z = B * V */
36  Malloc(Z, (msteps+1)*n_l, double);
37  } else {
38  /* Z and V are the same */
39  Z = V;
40  }
41  /* init vector */
42  if (ifGenEv) {
43  /* B norm */
44  matvec_B(v, Z);
45  t = 1.0 / sqrt(DDOT(&n, v, &one, Z, &one));
46  DSCAL(&n, &t, Z, &one);
47  } else {
48  /* 2-norm */
49  t = 1.0 / DNRM2(&n, v, &one);
50  }
51  /* starting vector */
52  DCOPY(&n, v, &one, V, &one);
53  /* unit B-norm or 2-norm */
54  DSCAL(&n, &t, V, &one);
55  double wn = 0.0;
56  /*-------------------- main Lanczos loop */
57  for (j=0; j<msteps; j++) {
58  int i;
59  /* znew = A*v */
60  /* vnew = A*v */
61  matvec_A(&V[j*n_l], &Z[(j+1)*n_l]);
62  /* znew = znew - bet * zold */
63  /* vnew = vnew - bet * vold */
64  if (j) {
65  nbet = -bet[j-1];
66  DAXPY(&n, &nbet, &Z[(j-1)*n_l], &one, &Z[(j+1)*n_l], &one);
67  }
68  /* alpha */
69  /* alp = znew' * v */
70  /* alp = vnew' * v */
71  alp[j] = DDOT(&n, &Z[(j+1)*n_l], &one, &V[j*n_l], &one);
72  wn += alp[j] * alp[j];
73  /* znew = znew - alp * z */
74  /* vnew = vnew - alp * v */
75  nalp = -alp[j];
76  DAXPY(&n, &nalp, &Z[j*n_l], &one, &Z[(j+1)*n_l], &one);
77  /* full reortho for znew */
78  for (i=0; i<=j; i++) {
79  /* (znew, v) */
80  /* (vnew, v) */
81  t = DDOT(&n, &Z[(j+1)*n_l], &one, &V[i*n_l], &one);
82  double mt = -t;
83  DAXPY(&n, &mt, &Z[i*n_l], &one, &Z[(j+1)*n_l], &one);
84  }
85  if (ifGenEv) {
86  /* vnew = B \ znew */
87  solve_B(&Z[(j+1)*n_l], &V[(j+1)*n_l]);
88  }
89  /* beta = (vnew, znew) */
90  /* beta = (vnew, vnew) */
91  bet[j] = DDOT(&n, &V[(j+1)*n_l], &one, &Z[(j+1)*n_l], &one);
92  if (bet[j]*(j+1) < orthTol*wn) {
93  fprintf(stdout, "lanbounds: lucky break, j=%d, beta=%e, break\n", j, bet[j]);
94  msteps = j + 1;
95  break;
96  }
97  wn += 2.0 * bet[j];
98  bet[j] = sqrt(bet[j]);
99  t = 1.0 / bet[j];
100  /* vnew = vnew / bet */
101  DSCAL(&n, &t, &V[(j+1)*n_l], &one);
102  if (ifGenEv) {
103  /* znew = znew / bet */
104  DSCAL(&n, &t, &Z[(j+1)*n_l], &one);
105  }
106  } /* main loop */
107
108  double bottomBeta = bet[msteps-1];
109  double *S, *ritzVal;
110  size_t msteps_l = msteps;
111  Malloc(S, msteps_l*msteps_l, double);
112  Malloc(ritzVal, msteps, double);
113  /*-------------------- diagonalize tridiagonal matrix */
114  SymmTridEig(ritzVal, S, msteps, alp, bet);
115 #if 1
116  *lmin = ritzVal[0] - fabs(bottomBeta * S[msteps-1]);
117  *lmax = ritzVal[msteps-1] + fabs(bottomBeta * S[msteps_l*msteps_l-1]);
118 #else
119  /*-------------------- 'safe' bounds */
120  double amin, amax, x;
121  amax = -INFINITY;
122  amin = INFINITY;
123  for (j=0; j<msteps; j++) {
124  t = fabs(bottomBeta * S[(j+1)*msteps-1]);
125  x = ritzVal[j]-t;
126  if (x<amin) amin = x;
127  x = ritzVal[j]+t;
128  if (x>amax) amax = x;
129  }
130  *lmin = amin;
131  *lmax = amax;
132 #endif
133  /*-------------------- done */
134  free(alp);
135  free(bet);
136  free(V);
137  if (ifGenEv) {
138  free(Z);
139  }
140  free(S);
141  free(ritzVal);
142
143  return 0;
144 }
145
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)
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
structs used in evsl
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
int LanBounds(int msteps, double *v, double *lmin, double *lmax)
Definition: lanbounds.c:21
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)