EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
dos_utils.c
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include "blaslapack.h"
6 #include "def.h"
7 #include "evsl.h"
8 #include "internal_proto.h"
9 #include "struct.h"
10 
15 double rec(const double a) { return 1.0 / a; } // Reciprocal
16 
17 double isqrt(const double a) { return 1.0 / sqrt(a); } // Inverse square root
18 
19 /*
20  * Initalize BSolDataPol: use L-S polynomial to approximate a function 'ffun'
21  */
22 void SetupBPol(int n, int max_deg, double tol, double lmin, double lmax,
23  double (*ffun)(double), BSolDataPol *data) {
24  data->max_deg = max_deg;
25  data->intv[0] = lmin;
26  data->intv[1] = lmax;
27  data->tol = tol;
28  Calloc(data->mu, max_deg, double);
29  Malloc(data->wk, 3*n, double);
30  /* determine degree and compute poly coeff */
31  lsPol(ffun, data);
32 }
33 
34 /*
35  * Initialize the member of BSolDataPol struct for solving B^{-1}
36  */
37 void SetupPolRec(int n, int max_deg, double tol, double lmin, double lmax,
38  BSolDataPol *data) {
39  SetupBPol(n, max_deg, tol, lmin, lmax, rec, data);
40 }
41 
42 /*
43  * Initialize the member of BSolDataPol struct for solving B^{1/2}
44  */
45 void SetupPolSqrt(int n, int max_deg, double tol, double lmin, double lmax,
46  BSolDataPol *data) {
47  SetupBPol(n, max_deg, tol, lmin, lmax, isqrt, data);
48 }
49 
50 /*
51  * Free the BSolDataPol struct
52  */
54  free(data->wk);
55  free(data->mu);
56 }
57 
58 /*
59  * Setup the function pointer for evsl struct to call B_sol function
60  */
61 void BSolPol(double *b, double *x, void *data) {
62  BSolDataPol *pol = (BSolDataPol *)data;
63  pnav(pol->mu, pol->deg, pol->cc, pol->dd, b, x, pol->wk);
64 }
65 
83 int apfun(const double c, const double h, const double *const xi,
84  double (*ffun)(double), const int npts, double *yi) {
85  int i = 0;
86  for (i = 0; i < npts; i++) {
87  yi[i] = ffun(c + h * xi[i]);
88  }
89  return 0;
90 }
91 
109 int pnav(double *mu, const int m, const double cc, const double dd, double *v,
110  double *y, double *w) { // Really just ChebAv
111  const int n = evsldata.n;
112  /*-------------------- pointers to v_[k-1],v_[k], v_[k+1] from w */
113  double *vk = w;
114  double *vkp1 = w + n;
115  double *vkm1 = vkp1 + n;
116  /*-------------------- */
117  int k, i;
118  double t, s, *tmp;
119  const double t1 = 1.0 / dd;
120  const double t2 = 2.0 / dd;
121  /*-------------------- vk <- v; vkm1 <- zeros(n,1) */
122  memcpy(vk, v, n * sizeof(double));
123  memset(vkm1, 0, n * sizeof(double));
124  /*-------------------- special case: k == 0 */
125  s = mu[0];
126  for (i = 0; i < n; i++) {
127  y[i] = s * vk[i];
128  }
129  /*-------------------- degree loop. k IS the degree */
130  for (k = 1; k <= m; k++) {
131  /*-------------------- y = mu[k]*Vk + y */
132  t = k == 1 ? t1 : t2;
133  /*-------------------- */
134  s = mu[k];
135  matvec_B(vk, vkp1);
136 
137  for (i = 0; i < n; i++) {
138  vkp1[i] = t * (vkp1[i] - cc * vk[i]) - vkm1[i];
139  /*-------------------- for degree 2 and up: */
140  y[i] += s * vkp1[i];
141  }
142  /*-------------------- next: rotate vectors via pointer exchange */
143  tmp = vkm1;
144  vkm1 = vk;
145  vk = vkp1;
146  vkp1 = tmp;
147  }
148 
149  return 0;
150 }
151 
168 int lsPol(double (*ffun)(double), BSolDataPol *pol) {
169  const double a = pol->intv[0];
170  const double b = pol->intv[1];
171  pol->cc = (a + b) / 2;
172  pol->dd = (b - a) / 2;
173  /*------------------------- Number of points for Gauss-Chebyshev integration */
174  int maxDeg = pol->max_deg;
175  const int npts = maxDeg * 4;
176  double *theti;
177  Malloc(theti, npts, double);
178 
179  int i = 0;
180  for (i = 0; i < npts * 2; i += 2) {
181  theti[i / 2] = (i + 1) * (PI / (2 * npts));
182  }
183 
184  double *xi;
185  Malloc(xi, npts, double);
186 
187  for (i = 0; i < npts; i++) {
188  xi[i] = cos(theti[i]);
189  }
190 
191  double *gi;
192  Malloc(gi, npts, double);
193  apfun(pol->cc, pol->dd, xi, ffun, npts, gi);
194  free(xi);
195 
196  /* ----------------------- Degree loop
197  * ----------------------- Each coefficient generated by gauss-Chebyshev quadature */
198  const int ptsNrm = 50 * maxDeg; // Number of points for norm infinity
199  Malloc(xi, ptsNrm, double);
200  linspace(-1, 1, ptsNrm, xi);
201 
202  /* Compute f(x) once */
203  double *yx;
204  Malloc(yx, ptsNrm, double);
205  apfun(pol->cc, pol->dd, xi, ffun, ptsNrm, yx);
206 
207  double *yi;
208  Malloc(yi, npts, double);
209 
210  double *ya;
211  Malloc(ya, ptsNrm, double);
212 
213  double na;
214 
215  int k = 0;
216  const double t1 = 1.0 / npts;
217  const double t2 = 2.0 / npts;
218 
219  /* degree loop */
220  for (k = 0; k < maxDeg; k++) {
221  for (i = 0; i < npts; i++) {
222  yi[i] = cos(k * theti[i]);
223  }
224  const double t = k == 0 ? t1 : t2;
225  double sum = 0;
226  for (i = 0; i < npts; i++) {
227  sum += yi[i] * gi[i];
228  }
229  pol->mu[k] = t * sum;
230  chebxPltd(k, pol->mu, ptsNrm, xi, ya);
231  double nrm = 0;
232  for (i = 0; i < ptsNrm; i++) {
233  na = (ya[i] - yx[i]) / yx[i];
234  nrm += fabs(na);
235  }
236  if (nrm < pol->tol) {
237  k++;
238  break;
239  }
240  }
241  /* mu is of size k, so deg = k - 1 */
242  pol->deg = k - 1;
243 
244  free(xi);
245  free(yi);
246  free(ya);
247  free(yx);
248  free(theti);
249  free(gi);
250 
251  return 0;
252 }
253 
254 #if 0
255 /*
256  * Recover the eigenvectors This is needed for GEN_MM folder only not DOS folder
257  */
258 int scalEigVec(int n, int nev, double *Y, double* sqrtdiag) {
259  // Formula (2.19) in the paper. Y(:,i) is x in the paper and D^{1/2} is sqrtdiag here
260  // n: matrix size
261  // nev: number of eigenvectors V: n*nev
262  // V: computed eigenvectors
263  // sqrtdiag: square root of diagonal entries of B
264  int i, j;
265  double *v;
266  for (i=0; i<nev; i++){
267  v = &Y[i*n];
268  for (j=0; j<n; j++) {
269  v[j] = v[j]*sqrtdiag[j];
270  }
271  }
272  return 0;
273 }
274 #endif
275 
double tol
Definition: struct.h:183
int max_deg
Definition: struct.h:179
double * mu
Definition: struct.h:184
double cc
Definition: struct.h:182
defs in EVSL
int lsPol(double(*ffun)(double), BSolDataPol *pol)
Definition: dos_utils.c:168
void SetupBPol(int n, int max_deg, double tol, double lmin, double lmax, double(*ffun)(double), BSolDataPol *data)
Definition: dos_utils.c:22
double dd
Definition: struct.h:182
int chebxPltd(int m, double *mu, int npts, double *xi, double *yi)
function yi = chebxPltd computes yi = p_mu (xi),
Definition: chebpoly.c:106
int pnav(double *mu, const int m, const double cc, const double dd, double *v, double *y, double *w)
Definition: dos_utils.c:109
int apfun(const double c, const double h, const double *const xi, double(*ffun)(double), const int npts, double *yi)
Definition: dos_utils.c:83
int n
Definition: struct.h:165
#define Calloc(base, nmem, type)
Definition: def.h:32
This file contains function prototypes and constant definitions internally used in EVSL...
real *8 function ffun(x, y, z)
Definition: functns.f90:40
double rec(const double a)
Definition: dos_utils.c:15
double * wk
Definition: struct.h:185
double intv[2]
Definition: struct.h:181
void SetupPolSqrt(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
Definition: dos_utils.c:45
#define Malloc(base, nmem, type)
Definition: def.h:22
structs used in evsl
This file contains function prototypes and constant definitions for EVSL.
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
double isqrt(const double a)
Definition: dos_utils.c:17
void BSolPol(double *b, double *x, void *data)
Definition: dos_utils.c:61
void SetupPolRec(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
Definition: dos_utils.c:37
int scalEigVec(int n, int nev, double *Y, double *sqrtdiag)
void FreeBSolPolData(BSolDataPol *data)
Definition: dos_utils.c:53