EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
MMPLanR.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 <sys/stat.h>
6 #include "evsl.h"
7 #include "io.h"
8 #include "evsl_direct.h"
9 
10 #define max(a, b) ((a) > (b) ? (a) : (b))
11 #define min(a, b) ((a) < (b) ? (a) : (b))
12 
13 /*-------------------- Protos */
14 int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo);
15 int get_matrix_info(FILE *fmat, io_t *pio);
16 /*-------------------- End Protos */
17 
18 int main() {
19  /*--------------------------------------------------------------
20  * this tests the spectrum slicing idea for a generic matrix pair
21  * read in matrix format -- using
22  * Thick-Restarted Lanczos with polynomial filtering.
23  *-------------------------------------------------------------*/
24  int n = 0, i, j, npts, nslices, nvec, Mdeg, nev, mlan, max_its, ev_int, sl,
25  ierr, totcnt;
26  /* find the eigenvalues of A in the interval [a,b] */
27  double a, b, lmax, lmin, ecount, tol, *sli, *mu;
28  double xintv[4];
29  double *alleigs;
30  int *counts; // #ev computed in each slice
31  /* initial vector: random */
32  double *vinit;
33  polparams pol;
34  /*-------------------- matrices A, B: coo format and csr format */
35  cooMat Acoo, Bcoo;
36  csrMat Acsr, Bcsr;
37  /* slicer parameters */
38  Mdeg = 50;
39  nvec = 40;
40  mu = malloc((Mdeg + 1) * sizeof(double));
41  FILE *flog = stdout, *fmat = NULL;
42  FILE *fstats = NULL;
43  io_t io;
44  int numat, mat;
45  char line[MAX_LINE];
46  /*-------------------- Bsol */
47  void *Bsol;
48 #if CXSPARSE == 1
49  printf("-----------------------------------------\n");
50  printf("Note: You are using CXSparse for the direct solver. \n We recommend a more performance based direct solver for anything more than basic tests. \n SuiteSparse is supported with a makefile change. \n Using SuiteSparse can result in magnitudes faster times. \n\n");
51  printf("-----------------------------------------\n");
52 #endif
53 
54  /*-------------------- stopping tol */
55  tol = 1e-6;
56  /*-------------------- start EVSL */
57  EVSLStart();
58  /*------------------ file "matfile" contains paths to matrices */
59  if (NULL == (fmat = fopen("matfile", "r"))) {
60  fprintf(flog, "Can't open matfile...\n");
61  exit(2);
62  }
63  /*-------------------- read number of matrices ..*/
64  memset(line, 0, MAX_LINE);
65  if (NULL == fgets(line, MAX_LINE, fmat)) {
66  fprintf(flog, "error in reading matfile...\n");
67  exit(2);
68  }
69  if ((numat = atoi(line)) <= 0) {
70  fprintf(flog, "Invalid count of matrices...\n");
71  exit(3);
72  }
73  for (mat = 1; mat <= numat; mat++) {
74  if (get_matrix_info(fmat, &io) != 0) {
75  fprintf(flog, "Invalid format in matfile ...\n");
76  exit(5);
77  }
78  /*----------------input matrix and interval information -*/
79  fprintf(flog, "MATRIX A: %s...\n", io.MatNam1);
80  fprintf(flog, "MATRIX B: %s...\n", io.MatNam2);
81  a = io.a; // left endpoint of input interval
82  b = io.b; // right endpoint of input interval
83  nslices = io.n_intv;
84 
85  struct stat st = {0}; /* Make sure OUT directory exists */
86  if (stat("OUT", &st) == -1) {
87  mkdir("OUT", 0750);
88  }
89 
90  char path[1024]; // path to write the output files
91  strcpy(path, "OUT/KPM_MMPLanR_");
92  strcat(path, io.MatNam1);
93  fstats = fopen(path, "w"); // write all the output to the file io.MatNam
94  if (!fstats) {
95  printf(" failed in opening output file in OUT/\n");
96  fstats = stdout;
97  }
98  fprintf(fstats, "MATRIX A: %s...\n", io.MatNam1);
99  fprintf(fstats, "MATRIX B: %s...\n", io.MatNam2);
100  fprintf(fstats, "Partition the interval of interest [%f,%f] into %d slices\n",
101  a, b, nslices);
102  counts = malloc(nslices * sizeof(int));
103  sli = malloc((nslices + 1) * sizeof(double));
104  /*-------------------- Read matrix - case: COO/MatrixMarket formats */
105  if (io.Fmt > HB) {
106  ierr = read_coo_MM(io.Fname1, 1, 0, &Acoo);
107  if (ierr == 0) {
108  fprintf(fstats, "matrix read successfully\n");
109  // nnz = Acoo.nnz;
110  n = Acoo.nrows;
111  // printf("size of A is %d\n", n);
112  // printf("nnz of A is %d\n", nnz);
113  } else {
114  fprintf(flog, "read_coo error for A = %d\n", ierr);
115  exit(6);
116  }
117  ierr = read_coo_MM(io.Fname2, 1, 0, &Bcoo);
118  if (ierr == 0) {
119  fprintf(fstats, "matrix read successfully\n");
120  if (Bcoo.nrows != n) {
121  return 1;
122  }
123  // nnz = Bcoo.nnz;
124  // n = Bcoo.nrows;
125  // printf("size of B is %d\n", n);
126  // printf("nnz of B is %d\n", nnz);
127  } else {
128  fprintf(flog, "read_coo error for B = %d\n", ierr);
129  exit(6);
130  }
131  /*-------------------- conversion from COO to CSR format */
132  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
133  ierr = cooMat_to_csrMat(0, &Bcoo, &Bcsr);
134  }
135  if (io.Fmt == HB) {
136  fprintf(flog, "HB FORMAT not supported (yet) * \n");
137  exit(7);
138  }
139  alleigs = malloc(n * sizeof(double));
140  /*-------------------- use direct solve as the solver for B in eigenvalue
141  * computation*/
142  SetupBSolDirect(&Bcsr, &Bsol);
143  /*-------------------- set the solver for B and LT */
144  SetBSol(BSolDirect, Bsol);
145  SetLTSol(LTSolDirect, Bsol);
146  /*-------------------- set the left-hand side matrix A */
147  SetAMatrix(&Acsr);
148  /*-------------------- set the right-hand side matrix B */
149  SetBMatrix(&Bcsr);
150  /*-------------------- for generalized eigenvalue problem */
151  SetGenEig();
152  /*-------------------- step 0: get eigenvalue bounds */
153  //-------------------- initial vector
154  vinit = (double *)malloc(n * sizeof(double));
155  rand_double(n, vinit);
156  ierr = LanTrbounds(50, 200, 1e-12, vinit, 1, &lmin, &lmax, fstats);
157  fprintf(fstats, "Step 0: Eigenvalue bound s for B^{-1}*A: [%.15e, %.15e]\n",
158  lmin, lmax);
159  /*-------------------- interval and eig bounds */
160  xintv[0] = a;
161  xintv[1] = b;
162  xintv[2] = lmin;
163  xintv[3] = lmax;
164  /*-------------------- call kpmdos to get the DOS for dividing the
165  * spectrum*/
166  /*-------------------- define kpmdos parameters */
167  double t = evsl_timer();
168  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
169  t = evsl_timer() - t;
170  if (ierr) {
171  printf("kpmdos error %d\n", ierr);
172  return 1;
173  }
174  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n", t);
175  fprintf(fstats, " estimated eig count in interval: %.15e \n", ecount);
176  //-------------------- call splicer to slice the spectrum
177  npts = 10 * ecount;
178  fprintf(fstats, "DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n", Mdeg,
179  nvec, npts);
180  ierr = spslicer(sli, mu, Mdeg, xintv, nslices, npts);
181  if (ierr) {
182  printf("spslicer error %d\n", ierr);
183  return 1;
184  }
185  printf("==================== SLICES FOUND ====================\n");
186  for (j = 0; j < nslices; j++) {
187  printf(" %2d: [% .15e , % .15e]\n", j + 1, sli[j], sli[j + 1]);
188  }
189  //-------------------- # eigs per slice
190  ev_int = (int)(1 + ecount / ((double)nslices));
191  totcnt = 0;
192  //-------------------- For each slice call ChebLanTr
193  for (sl = 0; sl < nslices; sl++) {
194  printf("======================================================\n");
195  int nev2;
196  double *lam, *Y, *res;
197  int *ind;
198  //--------------------
199  a = sli[sl];
200  b = sli[sl + 1];
201  printf(" subinterval: [%.15e , %.15e]\n", a, b);
202  xintv[0] = a;
203  xintv[1] = b;
204  xintv[2] = lmin;
205  xintv[3] = lmax;
206  //-------------------- set up default parameters for pol.
207  set_pol_def(&pol);
208  // can change default values here e.g.
209  pol.damping = 2;
210  pol.thresh_int = 0.8;
211  pol.thresh_ext = 0.2;
212  // pol.max_deg = 300;
213  //-------------------- Now determine polymomial
214  find_pol(xintv, &pol);
215  fprintf(fstats, " polynomial [type = %d], deg %d, bar %e gam %e\n",
216  pol.type, pol.deg, pol.bar, pol.gam);
217  // save_vec(pol.deg+1, pol.mu, "OUT/mu.txt");
218  //-------------------- approximate number of eigenvalues wanted
219  nev = ev_int + 2;
220  //-------------------- Dimension of Krylov subspace and maximal iterations
221  mlan = max(4 * nev, 100);
222  mlan = min(mlan, n);
223  max_its = 3 * mlan;
224  //-------------------- the call ChebLanTr
225  ierr = ChebLanTr(mlan, nev, xintv, max_its, tol, vinit, &pol, &nev2, &lam,
226  &Y, &res, fstats);
227  if (ierr) {
228  printf("ChebLanTr error %d\n", ierr);
229  return 1;
230  }
231 
232  /* sort the eigenvals: ascending order
233  * ind: keep the orginal indices */
234  ind = (int *)malloc(nev2 * sizeof(int));
235  sort_double(nev2, lam, ind);
236  printf(" number of eigenvalues found: %d\n", nev2);
237  /* print eigenvalues */
238  fprintf(fstats, " Eigenvalues in [a, b]\n");
239  fprintf(fstats, " Computed [%d] ||Res||\n", nev2);
240  for (i = 0; i < nev2; i++) {
241  fprintf(fstats, "% .15e %.1e\n", lam[i], res[ind[i]]);
242  }
243  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - "
244  "- - - - - - - - - - - - - - - - - -\n");
245  memcpy(&alleigs[totcnt], lam, nev2 * sizeof(double));
246  totcnt += nev2;
247  counts[sl] = nev2;
248  //-------------------- free allocated space withing this scope
249  if (lam)
250  free(lam);
251  if (Y)
252  free(Y);
253  if (res)
254  free(res);
255  free_pol(&pol);
256  free(ind);
257  } // for (sl=0; sl<nslices; sl++)
258  //-------------------- free other allocated space
259  fprintf(fstats, " --> Total eigenvalues found = %d\n", totcnt);
260  sprintf(path, "OUT/EigsOut_KPM_MMPLanR_(%s_%s)", io.MatNam1, io.MatNam2);
261  FILE *fmtout = fopen(path, "w");
262  if (fmtout) {
263  for (j = 0; j < totcnt; j++)
264  fprintf(fmtout, "%.15e\n", alleigs[j]);
265  fclose(fmtout);
266  }
267  free(vinit);
268  free(sli);
269  free_coo(&Acoo);
270  free_csr(&Acsr);
271  free_coo(&Bcoo);
272  free_csr(&Bcsr);
273  FreeBSolDirectData(Bsol);
274  free(alleigs);
275  free(counts);
276  if (fstats != stdout) {
277  fclose(fstats);
278  }
279  /*-------------------- end matrix loop */
280  }
281  free(mu);
282  if (flog != stdout) {
283  fclose(flog);
284  }
285  fclose(fmat);
286  /*-------------------- finalize EVSL */
287  EVSLFinish();
288  return 0;
289 }
290 
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
double b
Definition: io.h:25
#define MAX_LINE
Definition: io.h:5
int SetGenEig()
Set the problem to generalized eigenvalue problem.
Definition: evsl.c:158
int ChebLanTr(int lanm, int nev, double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats)
Chebyshev polynomial filtering Lanczos process [Thick restart version].
Definition: cheblanTr.c:62
int SetAMatrix(csrMat *A)
Set the matrix A.
Definition: evsl.c:68
int SetupBSolDirect(csrMat *B, void **data)
Setup the B-sol by computing the Cholesky factorization of B.
Definition: evsl_cxsparse.c:57
void rand_double(int n, double *v)
Definition: vect.c:11
int Fmt
Definition: io.h:20
int main()
Definition: MMPLanR.c:18
int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr)
convert coo to csr
Definition: spmat.c:135
int LanTrbounds(int lanm, int maxit, double tol, double *vinit, int bndtype, double *lammin, double *lammax, FILE *fstats)
Lanczos process for eigenvalue bounds [Thick restart version].
Definition: lanTrbounds.c:38
int find_pol(double *intv, polparams *pol)
Sets the values in pol.
Definition: chebpoly.c:361
double gam
Definition: struct.h:64
#define max(a, b)
Definition: MMPLanR.c:10
int SetBMatrix(csrMat *B)
Set the B matrix.
Definition: evsl.c:83
char Fname1[MAX_LINE]
Definition: io.h:17
void sort_double(int n, double *v, int *ind)
Definition: vect.c:92
int EVSLStart()
Initialize evslData.
Definition: evsl.c:27
void FreeBSolDirectData(void *data)
Free solver data.
char Fname2[MAX_LINE]
Definition: io.h:19
void free_csr(csrMat *csr)
memory deallocation for csr matrix
Definition: spmat.c:91
char MatNam2[MaxNamLen]
Definition: io.h:20
double thresh_ext
Definition: struct.h:51
#define HB
Definition: io.h:7
void set_pol_def(polparams *pol)
set default values for polparams struct.
Definition: chebpoly.c:18
double thresh_int
Definition: struct.h:52
int SetLTSol(SolFuncR func, void *data)
Set the solve routine for LT.
Definition: evsl.c:185
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.
#define min(a, b)
Definition: MMPLanR.c:11
Definition: io.h:14
int EVSLFinish()
Finish EVSL.
Definition: evsl.c:48
int get_matrix_info(FILE *fmat, io_t *pio)
Definition: io.c:22
void LTSolDirect(double *b, double *x, void *data)
Solver function of L^{T} x = L^{-T}*b.
int SetBSol(SolFuncR func, void *data)
Set the solve routine and the associated data for B.
Definition: evsl.c:132
double a
Definition: io.h:24
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
int nrows
Definition: struct.h:17
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
Definition: io.c:66
parameters for polynomial filter
Definition: struct.h:45
void BSolDirect(double *b, double *x, void *data)
Solver function of B.
Definitions used for direct solver interface.
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 n_intv
Definition: io.h:23
char MatNam1[MaxNamLen]
Definition: io.h:18