EVSL  1.1.0 EigenValues Slicing Library
exDOS.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 "string.h" //for memset
7
8 /*-------------------- Protos */
9 void linspace(double a, double b, int num, double *arr);
10
22 int exDOS(double *vals, int n, int npts,
23  double *x, double *y, double *intv) {
24 /* sigma = coefficient for gaussian
25  we select sigma so that at end of subinterval
26  we have exp[-0.5 [H/sigma]^2] = 1/K. */
27  int i,j=0, one=1;
28  double h, xi, t, scaling;
29  const double lm = intv[2];
30  const double lM = intv[3];
31  const double kappa = 1.25;
32  const int M = min(n, 30);
33  const double H = (lM - lm) / (M - 1);
34  const double sigma = H / sqrt(8 * log(kappa));
35  const double a = intv[0];
36  const double b = intv[1];
37  double sigma2 = 2 * sigma * sigma;
38 /*-------------------- sigma related..
39  -------------------- if gaussian smaller than tol ignore point. */
40  const double tol = 1e-08;
41  double width = sigma * sqrt(-2.0 * log(tol));
42  /* ------------------- define width of sub-intervals for DOS */
43  linspace(a, b, npts, x);
44  h = x[1]-x[0];
45  memset(y, 0.0, npts*sizeof(double));
46 /*-------------------- scan all relevant eigenvalues and add to its
47  interval [if any] */
48  for (i=0; i<n; i++){
49  t = vals[i];
50  if (t<a || t>b)
51  continue;
52  /*-------------------- first point to consider */
53  j = max((int) ((t-a-width)/h),0);
54  /*xi = a+j*h;
55  ! check! */
56  for (xi = a+j*h; xi <= min(t+width,b) ; xi+=h)
57  y[j++] += exp(-(xi-t)*(xi-t)/sigma2);
58  }
59  if (j == 0)
60  return 1;
61 /*-------------------- scale dos - [n eigenvalues in all]
62  scaling 2 -- for plots due to the missing 1/sqrt[2*pi*sigma^2] factor.
63  However, note that this does not guarantee that
64  sum[y] * h = the correct number of eigenvalues
65  y = y / sqrt[pi*sigma2] ; */
66  scaling = 1.0 / (n*sqrt(sigma2*PI));
67  DSCAL(&npts, &scaling, y, &one);
68  return 0;
69 }
defs in EVSL
void DSCAL(int *n, double *a, double *x, int *incx)
#define min(a, b)
Definition: def.h:62
int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv)
Definition: exDOS.c:22
#define max(a, b)
Definition: def.h:56
Defs for blaslapack routines.
#define PI
Definition: def.h:15
void linspace(double a, double b, int num, double *arr)
Definition: vect.c:54