EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
LapRLanR.c
Go to the documentation of this file.
1 #include <complex.h>
2 #include <math.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <sys/stat.h>
7 #include "evsl.h"
8 #include "evsl_direct.h"
9 #include "io.h"
10 
11 #define max(a, b) ((a) > (b) ? (a) : (b))
12 #define min(a, b) ((a) < (b) ? (a) : (b))
13 
14 int findarg(const char *argname, ARG_TYPE type, void *val, int argc,
15  char **argv);
16 int lapgen(int nx, int ny, int nz, cooMat *Acoo);
17 int exeiglap3(int nx, int ny, int nz, double a, double b, int *m, double **vo);
18 
19 int main(int argc, char *argv[]) {
20  /*------------------------------------------------------------
21  generates a laplacean matrix on an nx x ny x nz mesh
22  and computes all eigenvalues in a given interval [a b]
23  The default set values are
24  nx = 41; ny = 53; nz = 1;
25  a = 0.4; b = 0.8;
26  nslices = 1 [one slice only]
27  other parameters
28  tol [tolerance for stopping - based on residual]
29  Mdeg = pol. degree used for DOS
30  nvec = number of sample vectors used for DOS
31  This uses:
32  Thick-restart Lanczos with rational filtering
33  ------------------------------------------------------------*/
34  int n, nx, ny, nz, i, j, npts, nslices, nvec, Mdeg, nev, mlan, max_its,
35  ev_int, sl, flg, ierr;
36  /* find the eigenvalues of A in the interval [a,b] */
37  double a, b, lmax, lmin, ecount, tol, *sli, *mu;
38  double xintv[4];
39  /* initial vector: random */
40  double *vinit;
41  /* parameters for rational filter */
42  int num = 1; // number of poles used for each slice
43  int pow = 2; // multiplicity of each pole
44  double beta = 0.01; // beta in the LS approximation
45 
46  struct stat st = {0}; /* Make sure OUT directory exists */
47  if (stat("OUT", &st) == -1) {
48  mkdir("OUT", 0750);
49  }
50 
51  FILE *fstats = NULL;
52  if (!(fstats = fopen("OUT/LapRLanR.out", "w"))) {
53  printf(" failed in opening output file in OUT/\n");
54  fstats = stdout;
55  }
56  /*-------------------- matrix A: coo format and csr format */
57  cooMat Acoo;
58  csrMat Acsr;
59 #if CXSPARSE == 1
60  printf("-----------------------------------------\n");
61  printf(
62  "Note: You are using CXSparse for the direct solver. \n We recommend a "
63  "more performance based direct solver for anything more than basic "
64  "tests. \n SuiteSparse is supported with a makefile change. \n Using "
65  "SuiteSparse can result in magnitudes faster times. \n\n ");
66  printf("-----------------------------------------\n");
67 #endif
68  /*-------------------- default values */
69  nx = 41;
70  ny = 53;
71  nz = 8;
72  a = 0.4;
73  b = 0.8;
74  nslices = 4;
75  //-----------------------------------------------------------------------
76  //-------------------- reset some default values from command line [Yuanzhe/]
77  /* user input from command line */
78  flg = findarg("help", NA, NULL, argc, argv);
79  if (flg) {
80  printf(
81  "Usage: ./testL.ex -nx [int] -ny [int] -nz [int] -a [double] -b "
82  "[double] -nslices [int]\n");
83  return 0;
84  }
85  findarg("nx", INT, &nx, argc, argv);
86  findarg("ny", INT, &ny, argc, argv);
87  findarg("nz", INT, &nz, argc, argv);
88  findarg("a", DOUBLE, &a, argc, argv);
89  findarg("b", DOUBLE, &b, argc, argv);
90  findarg("nslices", INT, &nslices, argc, argv);
91  fprintf(fstats, "used nx = %3d ny = %3d nz = %3d", nx, ny, nz);
92  fprintf(fstats, " [a = %4.2f b= %4.2f], nslices=%2d \n", a, b, nslices);
93  //-------------------- eigenvalue bounds set by hand.
94  lmin = 0.0;
95  lmax = nz == 1 ? 8.0 : 12.0;
96  xintv[0] = a;
97  xintv[1] = b;
98  xintv[2] = lmin;
99  xintv[3] = lmax;
100  tol = 1e-8;
101  n = nx * ny * nz;
102  /*-------------------- generate 2D/3D Laplacian matrix
103  * saved in coo format */
104  ierr = lapgen(nx, ny, nz, &Acoo);
105  /*-------------------- convert coo to csr */
106  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
107  /* output the problem settings */
108  fprintf(fstats,
109  "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
110  fprintf(fstats, "Laplacian: %d x %d x %d, n = %d, nnz = %d\n", nx, ny, nz, n,
111  Acoo.nnz);
112  fprintf(fstats, "Interval: [%20.15f, %20.15f] -- %d slices \n", a, b,
113  nslices);
114  fprintf(fstats, "Step 0: Eigenvalue bound s for A: [%.15e, %.15e]\n", lmin,
115  lmax);
116  /*-------------------- call kpmdos to get the DOS for dividing the spectrum*/
117  /*-------------------- define kpmdos parameters */
118  Mdeg = 300;
119  nvec = 60;
120  /*-------------------- start EVSL */
121  EVSLStart();
122  /*-------------------- set the left-hand side matrix A */
123  SetAMatrix(&Acsr);
124  //-------------------- call kpmdos
125  mu = malloc((Mdeg + 1) * sizeof(double));
126  double t = evsl_timer();
127  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
128  t = evsl_timer() - t;
129  if (ierr) {
130  printf("kpmdos error %d\n", ierr);
131  return 1;
132  }
133  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n", t);
134  fprintf(fstats, " estimated eig count in interval: %10.2e \n", ecount);
135  //-------------------- call splicer to slice the spectrum
136  npts = 10 * ecount;
137  sli = malloc((nslices + 1) * sizeof(double));
138  ierr = spslicer(sli, mu, Mdeg, xintv, nslices, npts);
139  if (ierr) {
140  printf("spslicer error %d\n", ierr);
141  return 1;
142  }
143  printf("==================== SLICES FOUND ====================\n");
144  for (j = 0; j < nslices; j++) {
145  printf(" %2d: [% .15e , % .15e]\n", j + 1, sli[j], sli[j + 1]);
146  }
147  //-------------------- # eigs per slice
148  ev_int = (int)(1 + ecount / ((double)nslices));
149  //-------------------- initial vector
150  vinit = (double *)malloc(n * sizeof(double));
151  rand_double(n, vinit);
152  //-------------------- For each slice call RatLanrTr
153  for (sl = 0; sl < nslices; sl++) {
154  printf("======================================================\n");
155  int nev2;
156  double *lam, *Y, *res;
157  int *ind;
158  int nev_ex;
159  double *lam_ex;
160  //--------------------
161  a = sli[sl];
162  b = sli[sl + 1];
163  printf(" subinterval: [%.4e , %.4e]\n", a, b);
164  double intv[4];
165  intv[0] = a;
166  intv[1] = b;
167  intv[2] = lmin;
168  intv[3] = lmax;
169  // find the rational filter on this slice
170  ratparams rat;
171  // setup default parameters for rat
172  set_ratf_def(&rat);
173  // change some default parameters here:
174  rat.pw = pow;
175  rat.num = num;
176  rat.beta = beta;
177  // now determine rational filter
178  find_ratf(intv, &rat);
179  // use direct solver function
180  void **solshiftdata = (void **)malloc(num * sizeof(void *));
181  /*------------ factoring the shifted matrices and store the factors */
182  SetupASIGMABSolDirect(&Acsr, NULL, num, rat.zk, solshiftdata);
183  /*------------ give the data to rat */
184  SetASigmaBSol(&rat, NULL, ASIGMABSolDirect, solshiftdata);
185  //-------------------- approximate number of eigenvalues wanted
186  nev = ev_int + 2;
187  //-------------------- Dimension of Krylov subspace and maximal iterations
188  mlan = max(4 * nev, 300);
189  mlan = min(mlan, n);
190  max_its = 3 * mlan;
191  //-------------------- RationalLanTr
192  ierr = RatLanTr(mlan, nev, intv, max_its, tol, vinit, &rat, &nev2, &lam, &Y,
193  &res, fstats);
194  if (ierr) {
195  printf("RatLanTr error %d\n", ierr);
196  return 1;
197  }
198  /* [compute residual] already computed in res */
199  /* sort the eigenvals: ascending order
200  * ind: keep the orginal indices */
201  ind = (int *)malloc(nev2 * sizeof(int));
202  sort_double(nev2, lam, ind);
203  /* compute exact eigenvalues */
204  exeiglap3(nx, ny, nz, a, b, &nev_ex, &lam_ex);
205  printf(" number of eigenvalues: %d, found: %d\n", nev_ex, nev2);
206  /* print eigenvalues */
207  fprintf(fstats,
208  " Eigenvalues in [a, b]\n");
209  fprintf(fstats, " Computed [%d] ||Res|| Exact [%d]",
210  nev2, nev_ex);
211  if (nev2 == nev_ex) {
212  fprintf(fstats, " Err");
213  }
214  fprintf(fstats, "\n");
215  for (i = 0; i < max(nev2, nev_ex); i++) {
216  if (i < nev2) {
217  fprintf(fstats, "% .15e %.1e", lam[i], res[ind[i]]);
218  } else {
219  fprintf(fstats, " ");
220  }
221  if (i < nev_ex) {
222  fprintf(fstats, " % .15e", lam_ex[i]);
223  }
224  if (nev2 == nev_ex) {
225  fprintf(fstats, " % .1e", lam[i] - lam_ex[i]);
226  }
227  fprintf(fstats, "\n");
228  if (i > 50) {
229  fprintf(fstats, " -- More not shown --\n");
230  break;
231  }
232  }
233  //-------------------- free allocated space withing this scope
234  if (lam) free(lam);
235  if (Y) free(Y);
236  if (res) free(res);
237  FreeASIGMABSolDirect(rat.num, solshiftdata);
238  free(solshiftdata);
239  free(ind);
240  free(lam_ex);
241  free_rat(&rat);
242  } // for (sl=0; sl<nslices; sl++)
243  //-------------------- free other allocated space
244  free(vinit);
245  free(sli);
246  free_coo(&Acoo);
247  free_csr(&Acsr);
248  free(mu);
249  fclose(fstats);
250  /*-------------------- finalize EVSL */
251  EVSLFinish();
252  return 0;
253 }
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
void free_coo(cooMat *coo)
memory deallocation for coo matrix
Definition: spmat.c:126
void ASIGMABSolDirect(int n, double *br, double *bi, double *xr, double *xz, void *data)
complex linear solver routine passed to evsl
int main(int argc, char *argv[])
Definition: LapRLanR.c:19
complex double * zk
Definition: struct.h:126
int num
Definition: struct.h:114
void free_rat(ratparams *rat)
Definition: ratfilter.c:424
ARG_TYPE
Definition: io.h:33
int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num, complex double *zk, void **data)
setup CXsparse solver for A - SIGMA B
int nnz
Definition: struct.h:17
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 cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr)
convert coo to csr
Definition: spmat.c:135
Definition: io.h:34
void FreeASIGMABSolDirect(int num, void **data)
free the data needed by CXSparse
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
int pw
Definition: struct.h:115
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
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 findarg(const char *argname, ARG_TYPE type, void *val, int argc, char **argv)
Definition: io.c:171
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 RatLanTr(int lanm, int nev, double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nev2, double **vals, double **W, double **resW, FILE *fstats)
RatLanTR filtering Lanczos process [Thick restart version].
Definition: ratlanTr.c:61
#define max(a, b)
Definition: LapRLanR.c:11
#define min(a, b)
Definition: LapRLanR.c:12
double beta
Definition: struct.h:117
void set_ratf_def(ratparams *rat)
Sets default values for ratparams struct.
Definition: ratfilter.c:345
int SetASigmaBSol(ratparams *rat, SolFuncC *func, SolFuncC allf, void **data)
Set the solve routine and the associated data for A-SIGMA*B if func == NULL, set all functions to be ...
Definition: evsl.c:168
Definitions used for direct solver interface.
int find_ratf(double *intv, ratparams *rat)
Definition: ratfilter.c:375
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
parameters for rational filter
Definition: struct.h:112
Definition: io.h:37