EVSL  1.1.0 EigenValues Slicing Library
LapPLanN.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <sys/stat.h>
5 #include <math.h>
6 #include "evsl.h"
7 #include "io.h"
8
9 #define max(a, b) ((a) > (b) ? (a) : (b))
10 #define min(a, b) ((a) < (b) ? (a) : (b))
11 int findarg(const char *argname, ARG_TYPE type, void *val, int argc, char **argv);
12 int lapgen(int nx, int ny, int nz, cooMat *Acoo);
13 int exeiglap3(int nx, int ny, int nz, double a, double b, int *m, double **vo);
14
15 int main(int argc, char *argv[]) {
16  /*------------------------------------------------------------
17  generates a laplacean matrix on an nx x ny x nz mesh
18  and computes all eigenvalues in a given interval [a b]
19  The default set values are
20  nx = 41; ny = 53; nz = 1;
21  a = 0.4; b = 0.8;
22  nslices = 1 [one slice only]
23  other parameters
24  tol [tolerance for stopping - based on residual]
25  Mdeg = pol. degree used for DOS
26  nvec = number of sample vectors used for DOS
27  This uses:
28  Non-restart Lanczos with polynomial filtering
29  ------------------------------------------------------------*/
30  int n, nx, ny, nz, i, j, npts, nslices, nvec, Mdeg, nev,
31  mlan, ev_int, sl, flg, ierr;
32  /* find the eigenvalues of A in the interval [a,b] */
33  double a, b, lmax, lmin, ecount, tol, *sli, *mu;
34  double xintv[4];
35  double *vinit;
36  polparams pol;
37
38  struct stat st = {0}; /* Make sure OUT directory exists */
39  if (stat("OUT", &st) == -1) {
40  mkdir("OUT", 0750);
41  }
42
43  FILE *fstats = NULL;
44  if (!(fstats = fopen("OUT/LapPLanN.out","w"))) {
45  printf(" failed in opening output file in OUT/\n");
46  fstats = stdout;
47  }
48  /*-------------------- matrix A: coo format and csr format */
49  cooMat Acoo;
50  csrMat Acsr;
51  /*-------------------- default values */
52  nx = 16;
53  ny = 16;
54  nz = 20;
55  a = 0.4;
56  b = 0.8;
57  nslices = 4;
58  //-----------------------------------------------------------------------
59  //-------------------- reset some default values from command line
60  /* user input from command line */
61  flg = findarg("help", NA, NULL, argc, argv);
62  if (flg) {
63  printf("Usage: ./testL.ex -nx [int] -ny [int] -nz [int] -a [double] -b [double] -nslices [int]\n");
64  return 0;
65  }
66  /* print output on screen */
67  flg = findarg("stdout", NA, NULL, argc, argv);
68  if (flg) {
69  fstats = stdout;
70  }
71  findarg("nx", INT, &nx, argc, argv);
72  findarg("ny", INT, &ny, argc, argv);
73  findarg("nz", INT, &nz, argc, argv);
74  findarg("a", DOUBLE, &a, argc, argv);
75  findarg("b", DOUBLE, &b, argc, argv);
76  findarg("nslices", INT, &nslices, argc, argv);
77  fprintf(fstats,"used nx = %3d ny = %3d nz = %3d",nx,ny,nz);
78  fprintf(fstats," [a = %f b = %f], nslices=%2d \n",a,b,nslices);
79  //-------------------- eigenvalue bounds set by hand.
80  lmin = 0.0;
81  lmax = nz == 1 ? 8.0 : 12.0;
82  xintv[0] = a;
83  xintv[1] = b;
84  xintv[2] = lmin;
85  xintv[3] = lmax;
86  tol = 1e-6;
87  n = nx * ny * nz;
88  /*-------------------- output the problem settings */
89  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
90  fprintf(fstats, "Laplacian: %d x %d x %d, n = %d\n", nx, ny, nz, n);
91  fprintf(fstats, "Interval: [%20.15f, %20.15f] -- %d slices \n", a, b, nslices);
92  /*-------------------- generate 2D/3D Laplacian matrix
93  * saved in coo format */
94  ierr = lapgen(nx, ny, nz, &Acoo);
95  /*-------------------- convert coo to csr */
96  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
97  /*-------------------- step 0: get eigenvalue bounds */
98  fprintf(fstats, "Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n", lmin, lmax);
99  /*-------------------- call landos to get the DOS for dividing the spectrum*/
100  /*-------------------- define landos parameters */
101  Mdeg = 300;
102  nvec = 60;
103  /*-------------------- start EVSL */
104  EVSLStart();
105  /*-------------------- set matrix A */
106  SetAMatrix(&Acsr);
107  /*-------------------- call LanDos */
108  mu = (double *) malloc((Mdeg+1)*sizeof(double));
109  double t = evsl_timer();
110 #if 0
111  //-------------------- number of points - determines fine-ness of slices
112  npts = 100*nslices;
113  double *xdos, *ydos;
114  xdos = (double *) malloc((npts)*sizeof(double));
115  ydos = (double *) malloc((npts)*sizeof(double));
116
117  //fprintf(stdout," %d %d \n",npts,nslices);
118  LanDos(nvec, Mdeg, npts, xdos, ydos, &ecount, xintv);
119  //fprintf(stdout," %f \n",ecount);
120
121  //for (j=0; j<npts;j++) {
122  // printf(" %10.4f %10.4f \n",xdos[j],ydos[j]);}
123
124  t = evsl_timer() - t;
125
126  fprintf(fstats, " Time to build DOS (Landos) was : %10.2f \n",t);
127  fprintf(fstats, " estimated eig count in interval: %.15e \n",ecount);
128  //-------------------- call splicer to slice the spectrum
129
130  sli = malloc((nslices+1)*sizeof(double));
131
132  fprintf(fstats,"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n",
133  Mdeg, nvec, npts);
134  spslicer2(xdos, ydos, nslices, npts, sli);
135  free (xdos);
136  free (ydos);
137 #else
138  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
139  t = evsl_timer() - t;
140  if (ierr) {
141  printf("kpmdos error %d\n", ierr);
142  return 1;
143  }
144  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n",t);
145  fprintf(fstats, " estimated eig count in interval: %.15e \n",ecount);
146  fprintf(stdout, " estimated eig count in interval: %.15e \n",ecount);
147  //-------------------- call splicer to slice the spectrum
148  npts = 10 * ecount;
149  sli = malloc((nslices+1)*sizeof(double));
150
151  fprintf(fstats,"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n",Mdeg, nvec, npts);
152  ierr = spslicer(sli, mu, Mdeg, xintv, nslices, npts);
153 #endif
154  //-------------------- slicing done
155  if (ierr) {
156  printf("spslicer error %d\n", ierr);
157  return 1;
158  }
159  printf("==================== SLICES FOUND ====================\n");
160  for (j=0; j<nslices;j++) {
161  printf(" %2d: [% .15e , % .15e]\n", j+1, sli[j],sli[j+1]);
162  }
163  //-------------------- # eigs per slice
164  ev_int = (int) (1 + ecount / ((double) nslices));
165  //-------------------- initial vector
166  vinit = (double *) malloc(n*sizeof(double));
167  rand_double(n, vinit);
168  //-------------------- debug only :
169  // save_vec(n, vinit, "OUT/vinit.mtx");
170  //-------------------- For each slice call ChebLanr
171  for (sl=0; sl<nslices; sl++){
172  printf("======================================================\n");
173  int nev2;
174  double *lam, *Y, *res;
175  int *ind;
176  int nev_ex;
177  double *lam_ex;
178  //--------------------
179  StatsReset();
180  a = sli[sl];
181  b = sli[sl+1];
182  printf(" subinterval: [%.4e , %.4e]\n", a, b);
183  //-------------------- approximate number of eigenvalues wanted
184  nev = ev_int+2;
185  //-------------------- Dimension of Krylov subspace
186  mlan = max(5*nev, 300);
187  mlan = min(mlan, n);
188  //-------------------- ChebLanTr
189  xintv[0] = a; xintv[1] = b;
190  xintv[2] = lmin; xintv[3] = lmax;
191  //-------------------- set up default parameters for pol.
192  set_pol_def(&pol);
193  //-------------------- this is to show how you can reset some of the
194  // parameters to determine the filter polynomial
195  pol.damping = 2;
196  //-------------------- use a stricter requirement for polynomial
197  pol.thresh_int = 0.8;
198  pol.thresh_ext = 0.2;
199  pol.max_deg = 3000;
200  // pol.deg = 20 //<< this will force this exact degree . not recommended
201  // it is better to change the values of the thresholds
202  // pol.thresh_ext and plot.thresh_int
203  //-------------------- Now determine polymomial to use
204  find_pol(xintv, &pol);
205
206  fprintf(fstats, " polynomial [type = %d], deg %d, bar %e gam %e\n",
207  pol.type, pol.deg, pol.bar, pol.gam);
208  //-------------------- then call ChenLanNr
209  ierr = ChebLanNr(xintv, mlan, tol, vinit, &pol, &nev2, &lam, &Y, &res, fstats);
210  if (ierr) {
211  printf("ChebLanNr error %d\n", ierr);
212  return 1;
213  }
214
215  /* [compute residual] already computed in res */
216  /* sort the eigenvals: ascending order
217  * ind: keep the orginal indices */
218  ind = (int *) malloc(nev2*sizeof(int));
219  sort_double(nev2, lam, ind);
220  /* compute exact eigenvalues */
221  exeiglap3(nx, ny, nz, a, b, &nev_ex, &lam_ex);
222  printf(" number of eigenvalues: %d, found: %d\n", nev_ex, nev2);
223
224  /* print eigenvalues */
225  fprintf(fstats, " Eigenvalues in [a, b]\n");
226  fprintf(fstats, " Computed [%d] ||Res|| Exact [%d]", nev2, nev_ex);
227  if (nev2 == nev_ex) {
228  fprintf(fstats, " Err");
229  }
230  fprintf(fstats, "\n");
231  for (i=0; i<max(nev2, nev_ex); i++) {
232  if (i < nev2) {
233  fprintf(fstats, "% .15e %.1e", lam[i], res[ind[i]]);
234  } else {
235  fprintf(fstats, " ");
236  }
237  if (i < nev_ex) {
238  fprintf(fstats, " % .15e", lam_ex[i]);
239  }
240  if (nev2 == nev_ex) {
241  fprintf(fstats, " % .1e", lam[i]-lam_ex[i]);
242  }
243  fprintf(fstats,"\n");
244  if (i>50) {
245  fprintf(fstats," -- More not shown --\n");
246  break;
247  }
248  }
249  //-------------------- free allocated space withing this scope
250  if (lam) free(lam);
251  if (Y) free(Y);
252  if (res) free(res);
253  free_pol(&pol);
254  free(ind);
255  free(lam_ex);
256  StatsPrint(fstats);
257  }
258  //-------------------- free other allocated space
259  free(vinit);
260  free(sli);
261  free_coo(&Acoo);
262  free_csr(&Acsr);
263  free(mu);
264  fclose(fstats);
265  /*-------------------- finalize EVSL */
266  EVSLFinish();
267  return 0;
268 }
269
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
int deg
Definition: struct.h:70
void free_pol(polparams *pol)
Definition: chebpoly.c:488
void free_coo(cooMat *coo)
memory deallocation for coo matrix
Definition: spmat.c:126
int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos, double *neig, const double *const intv)
Definition: landos.c:40
ARG_TYPE
Definition: io.h:33
Definition: io.h:35
int SetAMatrix(csrMat *A)
Set the matrix A.
Definition: evsl.c:68
void StatsPrint(FILE *fstats)
Definition: stats.c:14
void rand_double(int n, double *v)
Definition: vect.c:11
int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr)
convert coo to csr
Definition: spmat.c:135
Definition: io.h:34
int find_pol(double *intv, polparams *pol)
Sets the values in pol.
Definition: chebpoly.c:361
double gam
Definition: struct.h:64
void sort_double(int n, double *v, int *ind)
Definition: vect.c:92
void spslicer2(double *xi, double *yi, int n_int, int npts, double *sli)
Definition: spslice2.c:29
int EVSLStart()
Initialize evslData.
Definition: evsl.c:27
int findarg(const char *argname, ARG_TYPE type, void *val, int argc, char **argv)
Definition: io.c:171
void free_csr(csrMat *csr)
memory deallocation for csr matrix
Definition: spmat.c:91
#define max(a, b)
Definition: LapPLanN.c:9
double thresh_ext
Definition: struct.h:51
#define min(a, b)
Definition: LapPLanN.c:10
void set_pol_def(polparams *pol)
set default values for polparams struct.
Definition: chebpoly.c:18
double thresh_int
Definition: struct.h:52
int main(int argc, char *argv[])
Definition: LapPLanN.c:15
int type
Definition: struct.h:59
sparse matrix format: the compressed sparse row (CSR) format, 0-based
Definition: struct.h:31
This file contains function prototypes and constant definitions for EVSL.
int max_deg
Definition: struct.h:48
int EVSLFinish()
Finish EVSL.
Definition: evsl.c:48
int ChebLanNr(double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
Definition: cheblanNr.c:60
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
int damping
Definition: struct.h:50
void StatsReset()
Definition: stats.c:79
parameters for polynomial filter
Definition: struct.h:45
int exeiglap3(int nx, int ny, int nz, double a, double b, int *m, double **vo)
Exact eigenvalues of Laplacean in interval [a b].
Definition: lapl.c:83
int lapgen(int nx, int ny, int nz, cooMat *Acoo)
Laplacean Matrix generator.
Definition: lapl.c:17
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
double bar
Definition: struct.h:65
Definition: io.h:37