EVSL  1.1.0 EigenValues Slicing Library
landos.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include "blaslapack.h"
5 #include "def.h"
6 #include "evsl.h"
7 #include "internal_proto.h"
8 #include "string.h" //for memset
9 #include "struct.h"
40 int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos,
41  double *neig, const double *const intv) {
42  //--------------------
43  const int ifGenEv = evsldata.ifGenEv;
44  if (ifGenEv) {
45  fprintf(stderr,
46  "ERROR: landos.c/LanDos only works with the standard eigenvalue "
47  "problem. Please use landosG.c/LanDos for the generalized problem. "
48  "\n");
49  exit(-1);
50  }
51  double *alp, *bet, nbet, nalp, t, *V;
52  int one = 1;
53  int n, m, i, j;
54  size_t n_l;
55
56  n = evsldata.n;
57  n_l = n;
58  //-------------------- Variables that persist through iterations
59  double *v, *y; // v=Vector for current iteration; y Stores y values
60  int *ind;
61  Malloc(v, n, double);
62  Calloc(y, npts, double);
63  Malloc(ind, npts, int);
64  /*-------------------- for tridiag. eigenvalue problem + lanczos
66  double *S, *ritzVal, *gamma2;
67  // S will contain a matrix compressed into a single array.
68  size_t msteps_l = msteps;
69  Malloc(S, msteps_l * msteps_l, double);
70  Malloc(gamma2, msteps, double);
71  Malloc(ritzVal, msteps, double);
72  const double lm = intv[2];
73  const double lM = intv[3];
74  const double tolBdwn = 1.e-13 * (abs(lM) + abs(lm));
75  const double aa = max(intv[0], intv[2]);
76  const double bb = min(intv[1], intv[3]);
77  const double kappa = 1.25;
78  const int M = min(msteps, 30);
79  const double H = (lM - lm) / (M - 1);
80  const double sigma = H / sqrt(8 * log(kappa));
81  double sigma2 = 2 * sigma * sigma;
82  //-------------------- If gaussian small than tol ignore point.
83  const double tol = 1e-08;
84  double width = sigma * sqrt(-2.0 * log(tol));
85  linspace(aa, bb, npts, xdos); // xdos = linspace(lm,lM, npts);
86
87  Malloc(alp, msteps, double);
88  Malloc(bet, msteps, double);
89  Malloc(V, (msteps + 1) * n_l, double);
90  //-------------------- Lanczos loop for this vector
91  for (m = 0; m < nvec; m++) {
92  randn_double(n, v); // w = randn(size(A,1),1);
93  //--------------------Start of bulk of lanbound.c code
94  t = DDOT(&n, v, &one, v, &one);
95  //-------------------- normalize vector
96  // v = can also use DNRM2 instead.
97  t = 1.0 / sqrt(t);
98  DSCAL(&n, &t, v, &one);
99  DCOPY(&n, v, &one, V, &one);
100  double wn = 0.0;
101  /*-------------------- main Lanczos loop */
102  for (j = 0; j < msteps; j++) {
103  // w = A*v
104  matvec_A(&V[j * n_l], &V[(j + 1) * n_l]);
105  // w = w - bet * vold
106  if (j) {
107  nbet = -bet[j - 1];
108  DAXPY(&n, &nbet, &V[(j - 1) * n_l], &one, &V[(j + 1) * n_l], &one);
109  }
110  /*-------------------- alp = w' * v */
111  alp[j] = DDOT(&n, &V[(j + 1) * n_l], &one, &V[j * n_l], &one);
112  wn += alp[j] * alp[j];
113  //-------------------- w = w - alp * v
114  nalp = -alp[j];
115  DAXPY(&n, &nalp, &V[j * n_l], &one, &V[(j + 1) * n_l], &one);
116  //-------------------- full reortho
117  for (i = 0; i <= j; i++) {
118  t = DDOT(&n, &V[(j + 1) * n_l], &one, &V[i * n_l], &one);
119  double mt = -t;
120  DAXPY(&n, &mt, &V[i * n_l], &one, &V[(j + 1) * n_l], &one);
121  }
122  bet[j] = DDOT(&n, &V[(j + 1) * n_l], &one, &V[(j + 1) * n_l], &one);
123  if (bet[j] * (j + 1) < orthTol * wn) {
124  fprintf(stdout, "lanbounds: lucky break, j=%d, beta=%e, break\n", j,
125  bet[j]);
126  msteps = j + 1;
127  break;
128  }
129  if (bet[j] > tolBdwn) { // If it's not zero, continue as normal
130  wn += 2.0 * bet[j];
131  bet[j] = sqrt(bet[j]);
132  t = 1.0 / bet[j];
133  DSCAL(&n, &t, &V[(j + 1) * n_l], &one);
134  } else { // Otherwise generate a new vector and redo the previous
135  // calculations on it
136  randn_double(n, v); // w = randn(size(A,1),1);
137  for (i = 0; i <= j; i++) {
138  t = DDOT(&n, &V[(j + 1) * n_l], &one, &V[i * n_l], &one);
139  double mt = -t;
140  DAXPY(&n, &mt, &V[i * n_l], &one, &V[(j + 1) * n_l], &one);
141  }
142  bet[j] = DDOT(&n, &V[(j + 1) * n_l], &one, &V[(j + 1) * n_l], &one);
143  wn += 2.0 * bet[j];
144  bet[j] = sqrt(bet[j]);
145  t = 1.0 / bet[j];
146  DSCAL(&n, &t, &V[(j + 1) * n_l], &one);
147  bet[j] = 0;
148  }
149  }
150  /*-------------------- end Lanczos loop for this vector
151  -------------------- diagonalize tridiagonal matrix */
152  SymmTridEig(ritzVal, S, msteps, alp, bet);
153  /* S = -eigvec
154  ritzVal = diags of D */
155
156  /* theta = ritzVal = sorted eigenvalues IN ASCENDING ORDER */
157  for (i = 0; i < msteps; i++) {
158  /*-------------------- weights for Lanczos quadrature
159  Gamma2(i) = elementwise square of top entry of i-th eginvector*/
160  gamma2[i] = S[i * msteps_l] * S[i * msteps_l];
161  }
162  /*-------------------- dos curve parameters
163  Generate DOS from small gaussians centered at the ritz values */
164  for (i = 0; i < msteps; i++) {
165  /* As msteps is width of ritzVal -> we get msteps eigenvectors */
166  const double t = ritzVal[i];
167  int numPlaced = 0;
168  /*-------------------- Place elements close to t in ind */
169  for (j = 0; j < npts; j++) {
170  if (fabs(xdos[j] - t) < width) ind[numPlaced++] = j;
171  }
172
173  for (j = 0; j < numPlaced; j++)
174  y[ind[j]] += gamma2[i] *
175  exp(-((xdos[ind[j]] - t) * (xdos[ind[j]] - t)) / sigma2);
176  }
177  /*-------------------- end vector loop */
178  }
179
180  double scaling = 1.0 / (nvec * sqrt(sigma2 * PI));
181  /* y = ydos * scaling */
182  DSCAL(&npts, &scaling, y, &one);
183  DCOPY(&npts, y, &one, ydos, &one);
184  simpson(xdos, y, npts);
185
186  *neig = y[npts - 1] * n;
187  free(gamma2);
188  free(S);
189  free(ritzVal);
190
191  free(alp);
192  free(bet);
193  free(V);
194
195  free(v);
196  free(y);
197  free(ind);
198
199  return 0;
200 }
int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos, double *neig, const double *const intv)
Definition: landos.c:40
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)
#define Calloc(base, nmem, type)
Definition: def.h:32
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 randn_double(int n, double *v)
Definition: vect.c:26
structs used in evsl
void DAXPY(int *n, double *alpha, double *x, int *incx, double *y, int *incy)
This file contains function prototypes and constant definitions for EVSL.
#define max(a, b)
Definition: def.h:56
Defs for blaslapack routines.
#define PI
Definition: def.h:15
evslData evsldata
global variable of EVSL
Definition: evsl.c:15
void linspace(double a, double b, int num, double *arr)
Definition: vect.c:54
void simpson(double *xi, double *yi, int npts)
Definition: simpson.c:27
void DCOPY(int *n, double *dx, int *incx, double *dy, int *incy)