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