EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
spslice.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdlib.h>
3 #include <stdio.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 
32 int kpmdos(int Mdeg, int damping, int nvec, double *intv,
33  double *mu, double *ecnt) {
34  /*-------------------- initialize variables */
35  const int ifGenEv = evsldata.ifGenEv;
36  int n = evsldata.n;
37  double *vkp1, *v, *vkm1, *vk, *jac;
38  double *w = NULL;
39  /*-------------------- workspace for generalized eigenvalue prob */
40  if (ifGenEv) {
41  Malloc(w, n, double);
42  }
43  Malloc(vkp1, n, double);
44  Malloc(v, n, double);
45  Malloc(vkm1, n, double);
46  Malloc(vk, n, double);
47  Malloc(jac, Mdeg+1, double);
48  double *tmp, ctr, wid;
49  double scal, t, tcnt, beta1, beta2, aa, bb;
50  int k, k1, i, m, mdegp1, one=1;
51  //-------------------- check if the interval is valid
52  if (check_intv(intv, stdout) < 0) {
53  return -1;
54  }
55 
56  aa = max(intv[0], intv[2]); bb = min(intv[1], intv[3]);
57  if (intv[0] < intv[2] || intv[1] > intv[3]) {
58  fprintf(stdout, " warning [%s (%d)]: interval (%e, %e) is adjusted to (%e, %e)\n",
59  __FILE__, __LINE__, intv[0], intv[1], aa, bb);
60  }
61  /*-------------------- some needed constants */
62  ctr = (intv[3]+intv[2])/2.0;
63  wid = (intv[3]-intv[2])/2.0;
64  t = max(-1.0+DBL_EPSILON, (aa-ctr)/wid);
65  beta1 = acos(t);
66  t = min(1.0-DBL_EPSILON, (bb-ctr)/wid);
67  beta2 = acos(t);
68  /*-------------------- compute damping coefs. */
69  dampcf(Mdeg, damping, jac);
70  /*-------------------- readjust jac[0] it was divided by 2 */
71  jac[0] = 1.0;
72  memset(mu,0,(Mdeg+1)*sizeof(double));
73  /*-------------------- seed the random generator */
74  i = time_seeder();
75  srand(i);
76  tcnt = 0.0;
77  /*-------------------- for random loop */
78  for (m=0; m<nvec; m++) {
79  if (ifGenEv) {
80  /* unit 2-norm v */
81  rand_double(n, v);
82  t = 1.0 / DNRM2(&n, v, &one);
83  DSCAL(&n, &t, v, &one);
84  /* w = L^{-T}*v */
85  solve_LT(v, w);
86  /* v = B*w */
87  matvec_B(w, v);
88  t = DDOT(&n, v, &one, w, &one);
89  memcpy(vk, w, n*sizeof(double));
90  } else {
91  /* unit 2-norm */
92  rand_double(n, v);
93  t = 1.0 / DNRM2(&n, v, &one);
94  DSCAL(&n, &t, v, &one);
95  memcpy(vk, v, n*sizeof(double));
96  }
97 
98  memset(vkm1, 0, n*sizeof(double));
99  mu[0] += jac[0];
100  //-------------------- for eigCount
101  tcnt -= jac[0]*(beta2-beta1);
102  /*-------------------- Chebyshev (degree) loop */
103  for (k=0; k<Mdeg; k++){
104  /*-------------------- Cheb. recurrence */
105  if (ifGenEv) {
106  /* v_{k+1} := B \ A * v_k (partial result) */
107  matvec_A(vk, w);
108  solve_B(w, vkp1);
109  } else {
110  matvec_A(vk, vkp1);
111  }
112  scal = k==0 ? 1.0 : 2.0;
113  scal /= wid;
114  for (i=0; i<n; i++) {
115  vkp1[i] = scal*(vkp1[i]-ctr*vk[i]) - vkm1[i];
116  }
117  //-------------------- rotate pointers to exchange vectors
118  tmp = vkm1;
119  vkm1 = vk;
120  vk = vkp1;
121  vkp1 = tmp;
122  /*-------------------- accumulate dot products for DOS expansion */
123  k1 = k+1;
124  t = 2*jac[k1] * DDOT(&n, vk, &one, v, &one);
125  mu[k1] += t;
126  /*-------------------- for eig. counts */
127  tcnt -= t*(sin(k1*beta2)-sin(k1*beta1))/k1;
128  }
129  }
130  //--------------------change of interval + scaling in formula
131  t = 1.0 /(((double)nvec)*PI) ;
132  mdegp1 = Mdeg+1;
133  DSCAL(&mdegp1, &t, mu, &one) ;
134  tcnt *= t*((double) n);
135  *ecnt = tcnt;
136  /*-------------------- free memory */
137  free(vkp1);
138  free(v);
139  free(vkm1);
140  free(vk);
141  free(jac);
142  if (ifGenEv) {
143  free(w);
144  }
145 
146  return 0;
147 }
148 
154 void intChx(const int Mdeg, double *mu, const int npts, double *xi, double *yi) {
155  //
156  int ndp1, j, k;
157  if(npts <= 0) {
158  fprintf(stderr, "Must have more than 0 points");
159  exit (1);
160  }
161  double val0, theta0, *thetas;
162  Malloc(thetas, npts, double);
163  ndp1 = Mdeg+1;
164  // if (xi[0]<-1.0) xi[0] = -1;
165  //if (xi[npts-1]> 1.0) xi[npts-1] = 1;
166 
167  for (j=0; j<npts; j++)
168  thetas[j] = acos(xi[j]);
169  theta0 = thetas[0];
170  for (j=0; j<npts; j++)
171  yi[j] = mu[0]*(theta0 - thetas[j]);
172  //-------------------- degree loop
173  for (k=1; k<ndp1; k++){
174  val0 = sin(k*theta0)/k;
175  //-------------------- points loop
176  for (j=0; j<npts; j++)
177  yi[j] += mu[k]*(val0 - sin(k*thetas[j])/k);
178  }
179  free (thetas);
180 }
181 
204 int spslicer(double *sli, double *mu, int Mdeg, double *intv, int n_int, int npts) {
205  int ls, ii, err=0;
206  double ctr, wid, aL, bL, target, aa, bb;
207 
208  if (check_intv(intv, stdout) < 0) {
209  return -1;
210  }
211 
212  // adjust a, b: intv[0], intv[1]
213  aa = max(intv[0], intv[2]); bb = min(intv[1], intv[3]);
214  if (intv[0] < intv[2] || intv[1] > intv[3]) {
215  fprintf(stdout, " warning [%s (%d)]: interval (%e, %e) is adjusted to (%e, %e)\n",
216  __FILE__, __LINE__, intv[0], intv[1], aa, bb);
217  }
218 
219  //--------------------
220  memset(sli,0,(n_int+1)*sizeof(double));
221  //-------------------- transform to ref interval [-1 1]
222  //-------------------- take care of trivial case n_int==1
223  if (n_int == 1){
224  sli[0] = intv[0];
225  sli[1] = intv[1];
226  return 0;
227  }
228  //-------------------- general case
229  ctr = (intv[3] + intv[2])/2;
230  wid = (intv[3] - intv[2])/2;
231  aL = (aa - ctr)/wid; // (a - ctr)/wid
232  bL = (bb - ctr)/wid; // (b - ctr)/wid
233  aL = max(aL,-1.0);
234  bL = min(bL,1.0);
235  npts = max(npts,2*n_int+1);
236  double *xi, *yi;
237  Malloc(xi, npts, double);
238  Malloc(yi, npts, double);
239  linspace(aL, bL, npts, xi);
240  //printf(" aL %15.3e bL %15.3e \n",aL,bL);
241  //-------------------- get all integrals at the xi's
242  //-------------------- exact integrals used.
243  intChx(Mdeg, mu, npts, xi, yi) ;
244  //-------------------- goal: equal share of integral per slice
245  target = yi[npts-1] / (double)n_int;
246  ls = 0;
247  ii = 0;
248  // use the unadjust left boundary
249  sli[ls] = intv[0];
250  //-------------------- main loop
251  while (++ls < n_int) {
252  while (ii < npts && yi[ii] < target) {
253  ii++;
254  }
255  if (ii == npts) {
256  break;
257  }
258  //-------------------- take best of 2 points in interval
259  if ( (target-yi[ii-1]) < (yi[ii]-target)) {
260  ii--;
261  }
262  sli[ls] = ctr + wid*xi[ii];
263  //-------------------- update target.. Slice size adjusted
264  target = yi[ii] + (yi[npts-1] - yi[ii])/(n_int-ls);
265  //printf("ls %d, n_int %d, target %e\n", ls, n_int, target);
266  }
267 
268  // use the unadjust left boundary
269  sli[n_int] = intv[1];
270 
271  //-------------------- check errors
272  if (ls != n_int) {
273  err = 1;
274  }
275  for (ii=1; ii<=n_int; ii++) {
276  if (sli[ii] < sli[ii-1]) {
277  err += 2;
278  break;
279  }
280  }
281 
282  if (err) {
283  printf("sli:\n");
284  for (ii=0; ii<=n_int; ii++) {
285  printf("%.15f\n", sli[ii]);
286  }
287  printf("\n");
288  save_vec(npts, xi, "OUT/xi.out");
289  save_vec(npts, yi, "OUT/yi.out");
290  }
291 
292  /*-------------------- free arrays */
293  free(xi);
294  free(yi);
295  return err;
296 }
297 
298 
299 
int kpmdos(int Mdeg, int damping, int nvec, double *intv, double *mu, double *ecnt)
This function computes the coefficients of the density of states in the chebyshev basis...
Definition: spslice.c:32
defs in EVSL
int ifGenEv
Definition: struct.h:166
int time_seeder()
Uses the timer to generate a seed to be used for srand.
Definition: mactime.c:30
double DDOT(int *n, double *x, int *incx, double *y, int *incy)
void rand_double(int n, double *v)
Definition: vect.c:11
void save_vec(int n, const double *x, const char fn[])
Definition: dumps.c:38
void intChx(const int Mdeg, double *mu, const int npts, double *xi, double *yi)
Computes the integrals where p(t) is the approximate DOS as given in the KPM method in the expanded ...
Definition: spslice.c:154
int n
Definition: struct.h:165
void DSCAL(int *n, double *a, double *x, int *incx)
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
#define max(a, b)
Definition: def.h:56
int spslicer(double *sli, double *mu, int Mdeg, double *intv, int n_int, int npts)
given the dos function defined by mu find a partitioning of sub-interval [a,b] of the spectrum so eac...
Definition: spslice.c:204
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
int dampcf(int m, int damping, double *jac)
Computes damping coefficient for cheb. expansions.
Definition: chebpoly.c:40