EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
LapPSI.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 
16 int main(int argc, char *argv[]) {
17  /*------------------------------------------------------------
18  generates a laplacean matrix on an nx x ny x nz mesh
19  and computes all eigenvalues in a given interval [a b]
20  The default set values are
21  nx = 41; ny = 53; nz = 1;
22  a = 0.4; b = 0.8;
23  nslices = 1 [one slice only]
24  other parameters
25  tol [tolerance for stopping - based on residual]
26  Mdeg = pol. degree used for DOS
27  nvec = number of sample vectors used for DOS
28  ...
29  This uses:
30  Subspace iteration with polynomial filtering
31  ------------------------------------------------------------*/
32  FILE *fstats = NULL;
33 
34  struct stat st = {0}; /* Make sure OUT directory exists */
35  if (stat("OUT", &st) == -1) {
36  mkdir("OUT", 0750);
37  }
38 
39  if (!(fstats = fopen("OUT/LapPSI.out","w"))) {
40  printf(" failed in opening output file in OUT/\n");
41  fstats = stdout;
42  }
43  int n, nx, ny, nz, i, j, npts, nslices, nvec, Mdeg, nev,
44  max_its, sl, flg, ierr;
45  //int ev_int;
46  /* find the eigenvalues of A in the interval [a,b] */
47  double a, b, lmax, lmin, ecount, tol, *sli, *mu;
48  double xintv[4];
49  double *vinit;
50  polparams pol;
51  /*-------------------- matrix A: coo format and csr format */
52  cooMat Acoo;
53  csrMat Acsr;
54  /*-------------------- default values */
55  nx = 41;
56  ny = 53;
57  nz = 1;
58  a = 0.4;
59  b = 0.8;
60  nslices = 4;
61  //-----------------------------------------------------------------------
62  //-------------------- reset some default values from command line
63  /* user input from command line */
64  flg = findarg("help", NA, NULL, argc, argv);
65  if (flg) {
66  printf("Usage: ./testL.ex -nx [int] -ny [int] -nz [int] -a [double] -b [double] -nslices [int]\n");
67  return 0;
68  }
69  findarg("nx", INT, &nx, argc, argv);
70  findarg("ny", INT, &ny, argc, argv);
71  findarg("nz", INT, &nz, argc, argv);
72  findarg("a", DOUBLE, &a, argc, argv);
73  findarg("b", DOUBLE, &b, argc, argv);
74  findarg("nslices", INT, &nslices, argc, argv);
75  fprintf(fstats,"used nx = %3d ny = %3d nz = %3d",nx,ny,nz);
76  fprintf(fstats," [ a = %4.2f b= %4.2f], nslices=%2d \n",a,b,nslices);
77  //-------------------- eigenvalue bounds set by hand.
78  lmin = 0.0;
79  lmax = nz == 1 ? 8.0 : 12.0;
80  xintv[0] = a;
81  xintv[1] = b;
82  xintv[2] = lmin;
83  xintv[3] = lmax;
84  tol = 1e-8;
85  n = nx * ny * nz;
86 
87  /* output the problem settings */
88  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
89  fprintf(fstats, "Laplacian: %d x %d x %d, n = %d\n", nx, ny, nz, n);
90  fprintf(fstats, "Interval: [%20.15f, %20.15f] -- %d slices \n", a, b, nslices);
91 
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: short Lanczos run to get bounds */
98  fprintf(fstats, "Step 0: Eigenvalue bound s for A: [%.15e, %.15e]\n", lmin, lmax);
99  /*-------------------- call kpmdos to get the DOS for dividing the spectrum*/
100  /*-------------------- define kpmdos parameters */
101  Mdeg = 300;
102  nvec = 60;
103  /*-------------------- start EVSL */
104  EVSLStart();
105  /*-------------------- set the left-hand side matrix A */
106  SetAMatrix(&Acsr);
107  //-------------------- call kpmdos
108  mu = malloc((Mdeg+1)*sizeof(double));
109  double t = evsl_timer();
110  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
111  t = evsl_timer() - t;
112  if (ierr) {
113  printf("kpmdos error %d\n", ierr);
114  return 1;
115  }
116  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n",t);
117  fprintf(fstats, " estimated eig count in interval: %10.2e \n",ecount);
118  //-------------------- call splicer to slice the spectrum
119  npts = 10 * ecount; //200;
120  sli = malloc((nslices+1)*sizeof(double));
121 
122  ierr = spslicer(sli, mu, Mdeg, xintv, nslices, npts);
123  if (ierr) {
124  printf("spslicer error %d\n", ierr);
125  return 1;
126  }
127  printf("==================== SLICES FOUND ====================\n");
128  for (j=0; j<nslices;j++)
129  printf(" %2d: [% .15e , % .15e]\n", j+1, sli[j],sli[j+1]);
130  //-------------------- # eigs per slice
131  // ev_int = (int) (1 + ecount / ((double) nslices));
132  //-------------------- initial vector
133  vinit = (double *) malloc(n*sizeof(double));
134  rand_double(n, vinit);
135  //-------------------- debug only :
136  // save_vec(n, vinit, "OUT/vinit.mtx");
137  //-------------------- For each slice call ChebLanr
138  // set polynomial defaults
139  // forced degree
140  for (sl=0; sl<nslices; sl++){
141  printf("======================================================\n");
142  int nevOut;
143  double *lam, *Y, *res;
144  int *ind;
145  int nev_ex;
146  double *lam_ex;
147  //--------------------
148  a = sli[sl];
149  b = sli[sl+1];
150  printf(" subinterval: [%.4e , %.4e]\n", a, b);
151  //-------------------- approximate number of eigenvalues wanted
152  // this factor insures that #eigs per slice is slightly overestimated
153  double fac = 1.2;
154  nev = (int) (1 + ecount / ((double) nslices)); // # eigs per slice
155  nev = (int)(fac*nev); // want an overestimate of ev_int
156  max_its = 1000;
157  //-------------------- ChebLanTr
158  xintv[0] = a; xintv[1] = b;
159  xintv[2] = lmin; xintv[3] = lmax;
160  set_pol_def(&pol);
161  //-------------------- this is to show how you can reset some of the
162  // parameters to determine the filter polynomial
163  pol.damping = 0;
164  pol.thresh_int = 0.5;
165  // pol.deg = 20 //<< this will force this exact degree . not recommended
166  // it is better to change the values of the thresholds
167  // pol.thresh_ext and plot.thresh_int
168  //-------------------- get polynomial to use :
169  find_pol(xintv, &pol);
170 
171  fprintf(fstats, " polynomial [type = %d], deg %d, bar %e gam %e\n",
172  pol.type, pol.deg, pol.bar, pol.gam);
173  //-------------------- then call ChenLanNr
174 
175  double *V0;
176  V0 = (double *) malloc(n*nev*sizeof(double));
177  rand_double(n*nev, V0);
178 
179  ierr = ChebSI(nev, xintv, max_its, tol, V0, &pol, &nevOut,
180  &lam, &Y, &res, fstats);
181  if (ierr) {
182  printf("ChebSI error %d\n", ierr);
183  return 1;
184  }
185 
186  free(V0);
187  /* sort the eigenvals: ascending order
188  * ind: keep the orginal indices */
189  ind = (int *) malloc(nevOut*sizeof(int));
190  sort_double(nevOut, lam, ind);
191  /* compute exact eigenvalues */
192  exeiglap3(nx, ny, nz, a, b, &nev_ex, &lam_ex);
193  printf(" number of eigenvalues: %d, found: %d\n", nev_ex, nevOut);
194 
195  /* print eigenvalues */
196  fprintf(fstats, " Eigenvalues in [a, b]\n");
197  fprintf(fstats, " Computed [%d] ||Res|| Exact [%d]",
198  nevOut, nev_ex);
199  if (nevOut == nev_ex) {
200  fprintf(fstats, " Err");
201  }
202  fprintf(fstats, "\n");
203  for (i=0; i<max(nevOut, nev_ex); i++) {
204  if (i < nevOut) {
205  fprintf(fstats, "% .15e %.1e", lam[i], res[ind[i]]);
206  } else {
207  fprintf(fstats, " ");
208  }
209  if (i < nev_ex) {
210  fprintf(fstats, " % .15e", lam_ex[i]);
211  }
212  if (nevOut == nev_ex) {
213  fprintf(fstats, " % .1e", lam[i]-lam_ex[i]);
214  }
215  fprintf(fstats,"\n");
216  if (i>50) {
217  fprintf(fstats," -- More not shown --\n");
218  break;
219  }
220  }
221  //-------------------- free allocated space withing this scope
222  if (lam) free(lam);
223  if (Y) free(Y);
224  if (res) free(res);
225  free_pol(&pol);
226  free(ind);
227  free(lam_ex);
228  }
229  //-------------------- free other allocated space
230  free(vinit);
231  free(sli);
232  free_coo(&Acoo);
233  free_csr(&Acsr);
234  free(mu);
235  fclose(fstats);
236  /*-------------------- finalize EVSL */
237  EVSLFinish();
238  return 0;
239 }
240 
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
#define max(a, b)
Definition: LapPSI.c:9
ARG_TYPE
Definition: io.h:33
Definition: io.h:35
int SetAMatrix(csrMat *A)
Set the matrix A.
Definition: evsl.c:68
void rand_double(int n, double *v)
Definition: vect.c:11
int findarg(const char *argname, ARG_TYPE type, void *val, int argc, char **argv)
Definition: io.c:171
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
int EVSLStart()
Initialize evslData.
Definition: evsl.c:27
int ChebSI(int nev, double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevo, double **lamo, double **Yo, double **reso, FILE *fstats)
Chebyshev polynomial filtering Subspace Iteration.
Definition: chebsi.c:42
void free_csr(csrMat *csr)
memory deallocation for csr matrix
Definition: spmat.c:91
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 main(int argc, char *argv[])
Definition: LapPSI.c:16
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
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 EVSLFinish()
Finish EVSL.
Definition: evsl.c:48
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
parameters for polynomial filter
Definition: struct.h:45
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
int lapgen(int nx, int ny, int nz, cooMat *Acoo)
Laplacean Matrix generator.
Definition: lapl.c:17
Definition: io.h:37