EVSL  1.1.0 EigenValues Slicing Library
KPM_MMRLanN.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <sys/stat.h>
5 #include <math.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 rational filtering.
23  *-------------------------------------------------------------*/
24  int n = 0, i, j, npts, nslices, nvec, Mdeg, nev, 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  /* parameters for rational filter */
34  int num = 1; // number of poles used for each slice
35  int pow = 2; // multiplicity of each pole
36  double beta = 0.01; // beta in the LS approximation
37  /*-------------------- matrices A, B: coo format and csr format */
38  cooMat Acoo, Bcoo;
39  csrMat Acsr, Bcsr, Acsr0, Bcsr0;
40  double *sqrtdiag = NULL;
41  /* slicer parameters */
42  Mdeg = 80;
43  nvec = 40;
44  mu = malloc((Mdeg + 1) * sizeof(double));
45  FILE *flog = stdout, *fmat = NULL;
46  FILE *fstats = NULL;
47  io_t io;
48  const int degB = 200; // Max degree to aproximate B with
49  const double tau = 1e-4; // Tolerance in polynomial approximation
50  int numat, mat;
51  char line[MAX_LINE];
52 #if CXSPARSE == 1
53  printf("-----------------------------------------\n");
54  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");
55  printf("-----------------------------------------\n");
56 #endif
57  /*-------------------- stopping tol */
58  tol = 1e-6;
59  /*-------------------- Polynomial approximation to B and sqrtB*/
60  BSolDataPol Bsol, Bsqrtsol;
61  /*-------------------- start EVSL */
62  EVSLStart();
63  /*------------------ file "matfile" contains paths to matrices */
64  if (NULL == (fmat = fopen("matfile", "r"))) {
65  fprintf(flog, "Can't open matfile...\n");
66  exit(2);
67  }
68  /*-------------------- read number of matrices ..*/
69  memset(line, 0, MAX_LINE);
70  if (NULL == fgets(line, MAX_LINE, fmat)) {
71  fprintf(flog, "error in reading matfile...\n");
72  exit(2);
73  }
74  if ((numat = atoi(line)) <= 0) {
75  fprintf(flog, "Invalid count of matrices...\n");
76  exit(3);
77  }
78  for (mat = 1; mat <= numat; mat++) {
79  if (get_matrix_info(fmat, &io) != 0) {
80  fprintf(flog, "Invalid format in matfile ...\n");
81  exit(5);
82  }
83  /*----------------input matrix and interval information -*/
84  fprintf(flog, "MATRIX A: %s...\n", io.MatNam1);
85  fprintf(flog, "MATRIX B: %s...\n", io.MatNam2);
86  a = io.a; // left endpoint of input interval
87  b = io.b; // right endpoint of input interval
88  nslices = io.n_intv;
89
90  struct stat st = {0}; /* Make sure OUT directory exists */
91  if (stat("OUT", &st) == -1) {
92  mkdir("OUT", 0750);
93  }
94
95  char path[1024]; // path to write the output files
96  strcpy(path, "OUT/KPM_MMRLanN_");
97  strcat(path, io.MatNam1);
98  fstats = fopen(path, "w"); // write all the output to the file io.MatNam
99  if (!fstats) {
100  printf(" failed in opening output file in OUT/\n");
101  fstats = stdout;
102  }
103  fprintf(fstats, "MATRIX A: %s...\n", io.MatNam1);
104  fprintf(fstats, "MATRIX B: %s...\n", io.MatNam2);
105  fprintf(fstats, "Partition the interval of interest [%f,%f] into %d slices\n",
106  a, b, nslices);
107  counts = malloc(nslices * sizeof(int));
108  sli = malloc((nslices + 1) * sizeof(double));
109  /*-------------------- Read matrix - case: COO/MatrixMarket formats */
110  if (io.Fmt > HB) {
111  ierr = read_coo_MM(io.Fname1, 1, 0, &Acoo);
112  if (ierr == 0) {
114  n = Acoo.nrows;
115  } else {
116  fprintf(flog, "read_coo error for A = %d\n", ierr);
117  exit(6);
118  }
119  ierr = read_coo_MM(io.Fname2, 1, 0, &Bcoo);
120  if (ierr == 0) {
122  if (Bcoo.nrows != n) {
123  return 1;
124  }
125  } else {
126  fprintf(flog, "read_coo error for B = %d\n", ierr);
127  exit(6);
128  }
129  /*------------------ diagonal scaling for Acoo and Bcoo */
130  sqrtdiag = (double *)calloc(n, sizeof(double));
131  /*------------------ conversion from COO to CSR format */
132  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
133  ierr = cooMat_to_csrMat(0, &Bcoo, &Bcsr);
134  } else if (io.Fmt == HB) {
135  fprintf(flog, "HB FORMAT not supported (yet) * \n");
136  exit(7);
137  }
138
139  /*-------------------- diagonal scaling for L-S poly. approx.
140  * of B^{-1} and B^{-1/2},
141  * which will be used in the DOS */
142  /*-------------------- sqrt of diag(B) */
143  extrDiagCsr(&Bcsr, sqrtdiag);
144  for (i=0; i<n; i++) {
145  sqrtdiag[i] = sqrt(sqrtdiag[i]);
146  }
147  /*-------------------- backup A and B */
148  csr_copy(&Acsr, &Acsr0, 1); /* 1 stands for memory alloc */
149  csr_copy(&Bcsr, &Bcsr0, 1);
150  /*-------------------- Scale A and B */
151  diagScalCsr(&Acsr, sqrtdiag);
152  diagScalCsr(&Bcsr, sqrtdiag);
153  if (sqrtdiag) {
154  free(sqrtdiag);
155  }
156
157  /*---------------- Set EVSL to solve std eig problem to
158  *---------------- compute the range of the spectrum of B */
159  SetStdEig();
160  SetAMatrix(&Bcsr);
161  vinit = (double *)malloc(n * sizeof(double));
162  rand_double(n, vinit);
163  ierr = LanTrbounds(50, 200, 1e-10, vinit, 1, &lmin, &lmax, fstats);
164  /*-------------------- Use polynomial to solve B and sqrt(B) */
165  /*-------------------- Setup the Bsol and Bsqrtsol struct */
166  SetupPolRec (n, degB, tau, lmin, lmax, &Bsol);
167  SetupPolSqrt(n, degB, tau, lmin, lmax, &Bsqrtsol);
168  SetBSol (BSolPol, (void *)&Bsol);
169  SetLTSol(BSolPol, (void *)&Bsqrtsol);
170  printf("The degree for LS polynomial approximations to B^{-1} and B^{-1/2} "
171  "are %d and %d\n", Bsol.deg, Bsqrtsol.deg);
172  /*-------------------- set the left-hand side matrix A */
173  SetAMatrix(&Acsr);
174  /*-------------------- set the right-hand side matrix B */
175  SetBMatrix(&Bcsr);
176  /*-------------------- for generalized eigenvalue problem */
177  SetGenEig();
178  /*-------------------- step 0: get eigenvalue bounds */
179  //-------------------- initial vector
180  rand_double(n, vinit);
181  ierr = LanTrbounds(50, 200, 1e-10, vinit, 1, &lmin, &lmax, fstats);
182  fprintf(fstats, "Step 0: Eigenvalue bound s for B^{-1}*A: [%.15e, %.15e]\n",
183  lmin, lmax);
184  /*-------------------- interval and eig bounds */
185  xintv[0] = a;
186  xintv[1] = b;
187  xintv[2] = lmin;
188  xintv[3] = lmax;
189  /*-------------------- call kpmdos to get the DOS for dividing the
190  * spectrum*/
191  /*-------------------- define kpmdos parameters */
192  //-------------------- call kpmdos
193  double t = evsl_timer();
194  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
195  t = evsl_timer() - t;
196  if (ierr) {
197  printf("kpmdos error %d\n", ierr);
198  return 1;
199  }
200  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n", t);
201  fprintf(fstats, " estimated eig count in interval: %.15e \n", ecount);
202  //-------------------- call splicer to slice the spectrum
203  npts = 10 * ecount;
204  fprintf(fstats, "DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n", Mdeg,
205  nvec, npts);
206  ierr = spslicer(sli, mu, Mdeg, xintv, nslices, npts);
207  if (ierr) {
208  printf("spslicer error %d\n", ierr);
209  return 1;
210  }
211  printf("==================== SLICES FOUND ====================\n");
212  for (j = 0; j < nslices; j++) {
213  printf(" %2d: [% .15e , % .15e]\n", j + 1, sli[j], sli[j + 1]);
214  }
215  /*-------------------- DONE WITH DOS */
216  FreeBSolPolData(&Bsol);
217  FreeBSolPolData(&Bsqrtsol);
218
219  //-------------------- # eigs per slice
220  ev_int = (int)(1 + ecount / ((double)nslices));
221  totcnt = 0;
222  alleigs = malloc(n * sizeof(double));
223
224  /* recover the original matrices A and B before scaling
225  * Note that B-sol and sqrt(B)-sol will not be needed in RatLan,
226  * so we can recover them */
227  csr_copy(&Acsr0, &Acsr, 0); /* 0 stands for no memory alloc */
228  csr_copy(&Bcsr0, &Bcsr, 0);
229
230  //-------------------- For each slice call RatLanrNr
231  for (sl = 0; sl < nslices; sl++) {
232  printf("======================================================\n");
233  int nev2;
234  double *lam, *Y, *res;
235  int *ind;
236  //--------------------
237  a = sli[sl];
238  b = sli[sl + 1];
239  printf(" subinterval: [%.4e , %.4e]\n", a, b);
240  double intv[4];
241  intv[0] = a;
242  intv[1] = b;
243  intv[2] = lmin;
244  intv[3] = lmax;
245  // find the rational filter on this slice
246  ratparams rat;
247  // setup default parameters for rat
248  set_ratf_def(&rat);
249  // change some default parameters here:
250  rat.pw = pow;
251  rat.num = num;
252  rat.beta = beta;
253  // now determine rational filter
254  find_ratf(intv, &rat);
255  // use direct solver function
256  void **solshiftdata = (void **)malloc(num * sizeof(void *));
257  /*------------ factoring the shifted matrices and store the factors */
258  SetupASIGMABSolDirect(&Acsr, &Bcsr, num, rat.zk, solshiftdata);
259  /*------------ give the data to rat */
260  SetASigmaBSol(&rat, NULL, ASIGMABSolDirect, solshiftdata);
261  //-------------------- approximate number of eigenvalues wanted
262  nev = ev_int + 2;
263  //-------------------- maximal Lanczos iterations
264  max_its = max(4 * nev, 300);
265  max_its = min(max_its, n);
266  //-------------------- RationalLanNr
267  ierr = RatLanNr(intv, max_its, tol, vinit, &rat, &nev2, &lam, &Y, &res,
268  fstats);
269  if (ierr) {
270  printf("RatLanNr error %d\n", ierr);
271  return 1;
272  }
273
274  /* sort the eigenvals: ascending order
275  * ind: keep the orginal indices */
276  ind = (int *)malloc(nev2 * sizeof(int));
277  sort_double(nev2, lam, ind);
278  printf(" number of eigenvalues found: %d\n", nev2);
279  /* print eigenvalues */
280  fprintf(fstats, " Eigenvalues in [a, b]\n");
281  fprintf(fstats, " Computed [%d] ||Res||\n", nev2);
282  for (i = 0; i < nev2; i++) {
283  fprintf(fstats, "% .15e %.1e\n", lam[i], res[ind[i]]);
284  }
285  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - "
286  "- - - - - - - - - - - - - - - - - -\n");
287  memcpy(&alleigs[totcnt], lam, nev2 * sizeof(double));
288  totcnt += nev2;
289  counts[sl] = nev2;
290  //-------------------- free allocated space withing this scope
291  if (lam)
292  free(lam);
293  if (Y)
294  free(Y);
295  if (res)
296  free(res);
297  FreeASIGMABSolDirect(rat.num, solshiftdata);
298  free(solshiftdata);
299  free_rat(&rat);
300  free(ind);
301  } // for (sl=0; sl<nslices; sl++)
302  //-------------------- free other allocated space
303  fprintf(fstats, " --> Total eigenvalues found = %d\n", totcnt);
304  sprintf(path, "OUT/EigsOut_KPM_MMRLanN_(%s_%s)", io.MatNam1, io.MatNam2);
305  FILE *fmtout = fopen(path, "w");
306  if (fmtout) {
307  for (j = 0; j < totcnt; j++)
308  fprintf(fmtout, "%.15e\n", alleigs[j]);
309  fclose(fmtout);
310  }
311  free(vinit);
312  free(sli);
313  free_coo(&Acoo);
314  free_csr(&Acsr);
315  free_coo(&Bcoo);
316  free_csr(&Bcsr);
317  free_csr(&Acsr0);
318  free_csr(&Bcsr0);
319  free(alleigs);
320  free(counts);
321  if (fstats != stdout) {
322  fclose(fstats);
323  }
324  /*-------------------- end matrix loop */
325  }
326  free(mu);
327  if (flog != stdout) {
328  fclose(flog);
329  }
330  fclose(fmat);
331  /*-------------------- finalize EVSL */
332  EVSLFinish();
333  return 0;
334 }
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
double b
Definition: io.h:25
complex double * zk
Definition: struct.h:126
#define MAX_LINE
Definition: io.h:5
int num
Definition: struct.h:114
void free_rat(ratparams *rat)
Definition: ratfilter.c:424
void extrDiagCsr(csrMat *B, double *d)
Definition: spmat.c:406
int SetGenEig()
Set the problem to generalized eigenvalue problem.
Definition: evsl.c:158
int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num, complex double *zk, void **data)
setup CXsparse solver for A - SIGMA B
int SetAMatrix(csrMat *A)
Set the matrix A.
Definition: evsl.c:68
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
void FreeASIGMABSolDirect(int num, void **data)
free the data needed by CXSparse
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
int pw
Definition: struct.h:115
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
#define HB
Definition: io.h:7
void diagScalCsr(csrMat *A, double *d)
Definition: spmat.c:390
void SetupPolSqrt(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
Definition: dos_utils.c:45
int RatLanNr(double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
Rational filtering Lanczos process [NON-restarted version].
Definition: ratlanNr.c:54
int SetLTSol(SolFuncR func, void *data)
Set the solve routine for LT.
Definition: evsl.c:185
int SetStdEig()
Set the problem to standard eigenvalue problem.
Definition: evsl.c:148
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 main()
Definition: KPM_MMRLanN.c:18
int EVSLFinish()
Finish EVSL.
Definition: evsl.c:48
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 get_matrix_info(FILE *fmat, io_t *pio)
Definition: io.c:22
int nrows
Definition: struct.h:17
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.
void csr_copy(csrMat *A, csrMat *B, int allocB)
copy a csr matrix A into B alloB: 0: will not allocate memory for B (have been alloced outside) 1: wi...
Definition: spmat.c:106
void BSolPol(double *b, double *x, void *data)
Definition: dos_utils.c:61
int find_ratf(double *intv, ratparams *rat)
Definition: ratfilter.c:375
void SetupPolRec(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
Definition: dos_utils.c:37
void FreeBSolPolData(BSolDataPol *data)
Definition: dos_utils.c:53
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
#define min(a, b)
Definition: KPM_MMRLanN.c:11
int n_intv
Definition: io.h:23
parameters for rational filter
Definition: struct.h:112
char MatNam1[MaxNamLen]
Definition: io.h:18
#define max(a, b)
Definition: KPM_MMRLanN.c:10
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
Definition: io.c:66