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