EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
landosG.c
Go to the documentation of this file.
1 #include <float.h>
2 #include <math.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include "blaslapack.h"
6 #include "def.h"
7 #include "evsl.h"
8 #include "internal_proto.h"
9 #include "string.h" //for memset
10 #include "struct.h"
11 
43 int LanDosG(const int nvec, const int msteps, int npts, double *xdos, double *ydos,
44  double *neig, const double *const intv) {
45 
46  int i, j, k;
47  int maxit = msteps, m; /* Max number of iterations */
48  /* size of the matrix */
49  int n = evsldata.n;
50  size_t n_l = n;
51  const int ifGenEv = evsldata.ifGenEv;
52 
53  double *vinit;
54  Malloc(vinit, n, double);
55 
56  int *ind;
57  Malloc(ind, npts, int);
58  double *y;
59  Calloc(y, npts, double);
60 
61  /*-------------------- frequently used constants */
62  int one = 1;
63  maxit = min(n, maxit);
64  size_t maxit_l = maxit;
65  double *gamma2;
66  Calloc(gamma2, maxit, double);
67  /*-----------------------------------------------------------------------*
68  * *Non-restarted* Lanczos iteration
69  *-----------------------------------------------------------------------
70  -------------------- Lanczos vectors V_m and tridiagonal matrix T_m */
71  double *V, *dT, *eT, *Z;
72  Calloc(V, n_l * (maxit + 1), double);
73  if (ifGenEv) {
74  /* storage for Z = B * V */
75  Calloc(Z, n_l * (maxit + 1), double);
76  } else {
77  /* Z and V are the same */
78  Z = V;
79  }
80  /*-------------------- diag. subdiag of Tridiagional matrix */
81  Malloc(dT, maxit, double);
82  Malloc(eT, maxit, double);
83  double *EvalT, *EvecT;
84  Malloc(EvalT, maxit, double); /* eigenvalues of tridia. matrix T */
85  Malloc(EvecT, maxit_l * maxit_l, double); /* Eigen vectors of T */
86  const double lm = intv[2];
87  const double lM = intv[3];
88  const double aa = max(intv[0], intv[2]);
89  const double bb = min(intv[1], intv[3]);
90  const double kappa = 1.25;
91  const int M = min(msteps, 30);
92  const double H = (lM - lm) / (M - 1);
93  const double sigma = H / sqrt(8 * log(kappa));
94  const double sigma2 = 2 * sigma * sigma;
95  /*-------------------- If gaussian small than tol ignore point. */
96  const double tol = 1e-08;
97  double width = sigma * sqrt(-2.0 * log(tol));
98  linspace(aa, bb, npts, xdos); // xdos = linspace(lm,lM, npts);
99  /*-------------------- u is just a pointer. wk == work space */
100  double *wk;
101  const size_t wk_size = ifGenEv ? 6 * n_l : 4 * n_l;
102  Malloc(wk, wk_size, double);
103  for (m = 0; m < nvec; m++) {
104  randn_double(n, vinit);
105  /*-------------------- copy initial vector to Z(:,1) */
106  /* Filter the initial vector */
107  if (ifGenEv) {
108  solve_LT(vinit, V);
109  }
110  /*-------------------- normalize it */
111  double t;
112  if (ifGenEv) {
113  /* B norm */
114  matvec_B(V, Z);
115  t = 1.0 / sqrt(DDOT(&n, V, &one, Z, &one));
116  DSCAL(&n, &t, Z, &one);
117  } else {
118  /* 2-norm */
119  t = 1.0 / DNRM2(&n, vinit, &one); // add a test here.
120  DCOPY(&n, vinit, &one, V, &one);
121  }
122  /* unit B^{-1}-norm or 2-norm */
123  DSCAL(&n, &t, V, &one);
124  /*-------------------- for ortho test */
125  double wn = 0.0;
126  int nwn = 0;
127  /*-------------------- lanczos vectors updated by rotating pointer*/
128  /*-------------------- pointers to Lanczos vectors */
129  double *zold, *z, *znew;
130  double *v, *vnew;
131  /*-------------------- Lanczos recurrence coefficients */
132  double alpha, nalpha, beta = 0.0, nbeta;
133  /* ---------------- main Lanczos loop */
134  for (k = 0; k < maxit; k++) {
135  /*-------------------- quick reference to Z(:,k-1) when k>0*/
136  zold = k > 0 ? Z + (k - 1) * n_l : NULL;
137  /*-------------------- a quick reference to V(:,k) */
138  v = &V[k * n_l];
139  /*-------------------- a quick reference to Z(:,k) */
140  z = &Z[k * n_l];
141  /*-------------------- next Lanczos vector V(:,k+1)*/
142  vnew = v + n;
143  /*-------------------- next Lanczos vector Z(:,k+1)*/
144  znew = z + n;
145  matvec_A(v, znew);
146  /*------------------ znew = znew - beta*zold */
147  if (zold) {
148  nbeta = -beta;
149  DAXPY(&n, &nbeta, zold, &one, znew, &one);
150  }
151  /*-------------------- alpha = znew'*v */
152  alpha = DDOT(&n, v, &one, znew, &one);
153  /*-------------------- T(k,k) = alpha */
154  dT[k] = alpha;
155  wn += fabs(alpha);
156  /*-------------------- znew = znew - alpha*z */
157  nalpha = -alpha;
158  DAXPY(&n, &nalpha, z, &one, znew, &one);
159  /*-------------------- FULL reortho to all previous Lan vectors */
160  if (ifGenEv) {
161  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
162  CGS_DGKS2(n, k, NGS_MAX, Z, V, znew, wk);
163  /* vnew = B \ znew */
164  solve_B(znew, vnew);
165  /*-------------------- beta = (vnew, znew)^{1/2} */
166  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
167  } else {
168  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
169  /* beta = norm(vnew) */
170  CGS_DGKS(n, k + 1, NGS_MAX, V, vnew, &beta, wk);
171  }
172  wn += 2.0 * beta;
173  nwn += 3;
174  /*-------------------- lucky breakdown test */
175  if (beta * nwn < orthTol * wn) {
176  rand_double(n, vnew);
177  if (ifGenEv) {
178  /* znew = znew - Z(:,1:k)*V(:,1:k)'*znew */
179  CGS_DGKS2(n, k + 1, NGS_MAX, V, Z, vnew, wk);
180  /* -------------- NOTE: B-matvec */
181  matvec_B(vnew, znew);
182  beta = sqrt(DDOT(&n, vnew, &one, znew, &one));
183  /*-------------------- vnew = vnew / beta */
184  t = 1.0 / beta;
185  DSCAL(&n, &t, vnew, &one);
186  /*-------------------- znew = znew / beta */
187  DSCAL(&n, &t, znew, &one);
188  beta = 0.0;
189  } else {
190  /* vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
191  /* beta = norm(vnew) */
192  CGS_DGKS(n, k + 1, NGS_MAX, V, vnew, &beta, wk);
193  /*-------------------- vnew = vnew / beta */
194  t = 1.0 / beta;
195  DSCAL(&n, &t, vnew, &one);
196  beta = 0.0;
197  }
198  } else {
199  /*-------------------- vnew = vnew / beta */
200  t = 1.0 / beta;
201  DSCAL(&n, &t, vnew, &one);
202  if (ifGenEv) {
203  /*-------------------- znew = znew / beta */
204  DSCAL(&n, &t, znew, &one);
205  }
206  }
207  /*-------------------- T(k+1,k) = beta */
208  eT[k] = beta;
209  }
210  SymmTridEig(EvalT, EvecT, maxit, dT, eT);
211  for (i = 0; i < maxit; i++) {
212  /*-------------------- weights for Lanczos quadrature */
213  /* Gamma2(i) = elementwise square of top entry of i-th eginvector
214  * */
215  gamma2[i] = EvecT[i * maxit_l] * EvecT[i * maxit_l];
216  }
217  /*-------------------- dos curve parameters
218  Generate DOS from small gaussians centered at the ritz values */
219  for (i = 0; i < msteps; i++) {
220  // As msteps is width of ritzVal -> we get msteps eigenvectors
221  const double t = EvalT[i];
222  int numPlaced = 0;
223  /*-------------------- Place elements close to t in ind */
224  for (j = 0; j < npts; j++) {
225  if (fabs(xdos[j] - t) < width) ind[numPlaced++] = j;
226  }
227  for (j = 0; j < numPlaced; j++) {
228  y[ind[j]] += gamma2[i] *
229  exp(-((xdos[ind[j]] - t) * (xdos[ind[j]] - t)) / sigma2);
230  }
231  }
232  }
233  double scaling = 1.0 / (nvec * sqrt(sigma2 * PI));
234  /* y = ydos * scaling */
235  DSCAL(&npts, &scaling, y, &one);
236  DCOPY(&npts, y, &one, ydos, &one);
237  simpson(xdos, y, npts);
238  *neig = y[npts - 1] * n;
239  free(gamma2);
240  /*-------------------- free arrays */
241  free(vinit);
242  free(V);
243  free(dT);
244  free(eT);
245  free(EvalT);
246  free(EvecT);
247  free(wk);
248  free(y);
249  free(ind);
250  if (ifGenEv) {
251  free(Z);
252  }
253 
254  return 0;
255 }
256 
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
#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
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
int LanDosG(const int nvec, const int msteps, int npts, double *xdos, double *ydos, double *neig, const double *const intv)
Definition: landosG.c:43
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)
This file contains function prototypes and constant definitions for EVSL.
#define max(a, b)
Definition: def.h:56
double DNRM2(int *n, double *x, int *incx)
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)