EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
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 <complex.h>
7 #include "evsl.h"
8 #include "io.h"
9 #include "evsl_direct.h"
10 
11 #define max(a, b) ((a) > (b) ? (a) : (b))
12 #define min(a, b) ((a) < (b) ? (a) : (b))
13 #define TRIV_SLICER 0
14 
15 int main () {
16  int ierr = 0;
17  /*--------------------------------------------------------------
18  * this tests the spectrum slicing idea for a generic matrix
19  * read in matrix format -- using
20  * Thick-Restarted Lanczos with rational filtering.
21  *-------------------------------------------------------------*/
22  int n=0, sl, i, j, nev, totcnt;
23  double a, b, ecount, xintv[4];
24  double lmin, lmax;
25  double *alleigs;
26  int n_intv; // number of subintervals (slices)
27  int npnts; // number of integration points for eigenvalue count
28  /*-------------------- matrix A: coo format and csr format */
29  cooMat Acoo;
30  csrMat Acsr;
31  /* tolerance to accept ev */
32  double tol;
33  /* total #ev computed; the computed eig val/vec */
34  int nevOut;
35  /*-------------------- max Lanczos steps */
36  int max_its;
37  /*-------------------- poly. deg. and # of vectors for kpmdos */
38  int Mdeg, nvec;
39  /*-------------------- IO */
40  FILE *flog = stdout, *fmat = NULL;
41  FILE *fstats = NULL;
42  io_t io;
43  int numat, mat;
44  char line[MAX_LINE];
45  /* initial vector: random */
46  double *vinit;
47  tol = 1e-6;
48  /* parameters for rational filter */
49  int pow = 2; // multiplicity of the pole
50  double beta = 0.01; // beta in the LS approximation
51  /* slicer parameters */
52  npnts = 1000;
53  Mdeg = 300;
54  nvec = 60;
55  /*-------------------- start EVSL */
56  EVSLStart();
57  /* interior eigensolver parameters */
58  double *mu = malloc((Mdeg+1)*sizeof(double)); // coeff. for kpmdos
59  int *counts; // #ev computed in each slice
60  double *sli; // endpoints of partitioned slices
61 #if CXSPARSE == 1
62  printf("-----------------------------------------\n");
63  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");
64  printf("-----------------------------------------\n");
65 #endif
66  /*------------------ file "matfile" contains paths to matrices */
67  if( NULL == ( fmat = fopen( "matfile", "r" ) ) ) {
68  fprintf( flog, "Can't open matfile...\n" );
69  exit(2);
70  }
71  /*-------------------- read number of matrices ..*/
72  memset( line, 0, MAX_LINE );
73  if (NULL == fgets( line, MAX_LINE, fmat )) {
74  fprintf( flog, "error in reading matfile...\n" );
75  exit(2);
76  }
77  if( ( numat = atoi( line ) ) <= 0 ) {
78  fprintf( flog, "Invalid count of matrices...\n" );
79  exit(3);
80  }
81  /*-------------------- LOOP through matrices -*/
82  for(mat = 1; mat <= numat; mat++) {
83  if(get_matrix_info(fmat, &io) != 0) {
84  fprintf(flog, "Invalid format in matfile ...\n");
85  exit(5);
86  }
87  /*----------------input matrix and interval information -*/
88  fprintf(flog, "MATRIX: %s...\n", io.MatNam);
89  a = io.a; // left endpoint of input interval
90  b = io.b; // right endpoint of input interval
91  n_intv = io.n_intv;
92 
93  struct stat st = {0}; /* Make sure OUT directory exists */
94  if (stat("OUT", &st) == -1) {
95  mkdir("OUT", 0750);
96  }
97 
98  char path[1024]; // path to write the output files
99  strcpy( path, "OUT/MMRLanN_");
100  strcat( path, io.MatNam);
101  fstats = fopen(path,"w"); // write all the output to the file io.MatNam
102  if (!fstats) {
103  printf(" failed in opening output file in OUT/\n");
104  fstats = stdout;
105  }
106  fprintf(fstats, "MATRIX: %s...\n", io.MatNam);
107  fprintf(fstats,"Partition the interval of interest [%f,%f] into %d slices\n",
108  a,b,n_intv);
109  counts = malloc(n_intv*sizeof(int));
110  sli = malloc( (n_intv+1)*sizeof(double));
111  /*-------------------- Read matrix - case: COO/MatrixMarket formats */
112  if (io.Fmt > HB){
113  ierr =read_coo_MM(io.Fname, 1, 0, &Acoo);
114  if (ierr == 0) {
115  fprintf(fstats,"matrix read successfully\n");
116  n = Acoo.nrows;
117  }
118  else {
119  fprintf(flog, "read_coo error = %d\n", ierr);
120  exit(6);
121  }
122  /*-------------------- conversion from COO to CSR format */
123  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
124  }
125  if (io.Fmt == HB) {
126  fprintf(flog, "HB FORMAT not supported (yet) * \n");
127  exit(7);
128  }
129  /*-------------------- set the left-hand side matrix A */
130  SetAMatrix(&Acsr);
131  /*-------------------- define parameters for DOS */
132  alleigs = malloc(n*sizeof(double));
133  vinit = (double *) malloc(n*sizeof(double));
134  rand_double(n, vinit);
135  /*-------------------- get lambda_min lambda_max estimates */
136  ierr = LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
137  fprintf(fstats, "Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n",
138  lmin, lmax);
139  /*-------------------- define [a b] now so we can get estimates now
140  on number of eigenvalues in [a b] from kpmdos */
141  fprintf(fstats," --> interval: a %9.3e b %9.3e \n",a, b);
142  /*-------------------- define kpmdos parameters */
143  xintv[0] = a; xintv[1] = b;
144  xintv[2] = lmin; xintv[3] = lmax;
145  //-------------------- trivial slicer or kpmdos?
146  if (TRIV_SLICER) {
147  linspace(a, b, n_intv+1, sli);
148  } else {
149  double t = evsl_timer();
150  ierr = kpmdos(Mdeg, 0, nvec, xintv, mu, &ecount);
151  t = evsl_timer() - t;
152  if (ierr) {
153  printf("kpmdos error %d\n", ierr);
154  return 1;
155  }
156  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n",t);
157  fprintf(fstats, "Step 1a: Estimated eig count in interval - %10.2e \n",
158  ecount);
159  if (ecount < 0 || ecount > n) {
160  printf(" e-count estimate is incorrect \n ");
161  exit(7);
162  }
163  /*-------------------- error in ecount */
164  /*
165  if ((ecount <1) | (ecount>n))
166  exit(6);
167  */
168  /*-------------------- define slicer parameters */
169  npnts = 10 * ecount;
170  ierr = spslicer(sli, mu, Mdeg, xintv, n_intv, npnts);
171  if (ierr) {
172  printf("spslicer error %d\n", ierr);
173  return 1;
174  }
175  }
176  fprintf(fstats,"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n",
177  Mdeg, nvec, npnts);
178  fprintf(fstats, "Step 1b: Slices found: \n");
179  for (j=0; j<n_intv;j++)
180  fprintf(fstats,"[% 12.4e , % 12.4e]\n", sli[j],sli[j+1]);
181  //-------------------- # eigs per slice
182  //-------------------- approximate number of eigenvalues wanted
183  double fac = 1.2; // this factor insures that # of eigs per slice is slightly overestimated
184  printf("Total number of eigenvalues estimated = %d \n", (int)(ecount));
185  nev = (int) (1 + ecount / ((double) n_intv)); // # eigs per slice
186  nev = (int)(fac*nev); // want an overestimate of ev_int
187  max_its = max(4*nev,300); max_its = min(max_its, n);
188  fprintf(fstats, "Step 2: In each slice compute %d eigenvalues ... \n", nev);
189  /*-------------------- MAIN intv LOOP: for each sclice Do: */
190  totcnt = 0;
191  for (sl=0; sl<n_intv; sl++) {
192  fprintf(fstats,"======================================================\n");
193  double *lam, *Y, *res;
194  int *ind;
195  //--------------------
196  a = sli[sl];
197  b = sli[sl+1];
198  StatsReset();
199  fprintf(fstats, " subinterval: [% 12.4e , % 12.4e]\n", a, b);
200  //-------------------- Parameters for RatLanNr
201  fprintf(fstats, " Max Lanczos steps %d\n", max_its);
202  double intv[4]; // endpoints of this slice and spectrum
203  intv[0] = a;
204  intv[1] = b;
205  intv[2] = lmin;
206  intv[3] = lmax;
207  printf("=== Compute the %dth subinterval: [%.4e, %.4e] out of %d ===\n",
208  sl+1, a, b, n_intv);
209  //------------Find the rational filter on this slice
210  ratparams rat;
211  //------------Set up default parameters for rat
212  set_ratf_def(&rat);
213  //------------Change some default values here:
214  rat.beta = beta;
215  rat.pow = pow;
216  //-------------Now determine rational filter
217  find_ratf(intv, &rat);
218  /*------------ use direct solver function */
219  void **solshiftdata = (void **) malloc(rat.num*sizeof(void *));
220  /*------------ factoring the shifted matrices and store the factors */
221  SetupASIGMABSolDirect(&Acsr, NULL, rat.num, rat.zk, solshiftdata);
222  /*------------ give the data to rat */
223  SetASigmaBSol(&rat, NULL, ASIGMABSolDirect, solshiftdata);
224  //-------------------- RationalLanNr
225  ierr = RatLanNr(intv, max_its, tol, vinit, &rat, &nevOut, &lam,
226  &Y, &res, fstats);
227  if (ierr) {
228  printf("RatLanNr error %d\n", ierr);
229  return 1;
230  }
231  /* sort the eigenvals: ascending order
232  * ind: keep the orginal indices */
233  ind = (int *) malloc(nevOut*sizeof(int));
234  sort_double(nevOut, lam, ind);
235 
236  /* print eigenvalues */
237  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
238  fprintf(fstats, " Eigenvalues in [%f, %f]\n",
239  a,b);
240  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
241  fprintf(fstats, " Computed [%d out of %d estimated] ||Res|| ",
242  nevOut, nev);
243  fprintf(fstats, "\n");
244  for (i=0; i<nevOut; i++) {
245  fprintf(fstats, " % .15e %.1e", lam[i], res[ind[i]]);
246  fprintf(fstats,"\n");
247  }
248  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
249  memcpy(&alleigs[totcnt],lam,nevOut*sizeof(double));
250  totcnt += nevOut;
251  counts[sl] = nevOut;
252  StatsPrint(stdout);
253  /*-------------------- free memory within this loop */
254  if (lam) free(lam);
255  if (Y) free(Y);
256  if (res) free(res);
257  free(ind);
258  FreeASIGMABSolDirect(rat.num, solshiftdata);
259  free(solshiftdata);
260  free_rat(&rat);
261  /*-------------------- end slice loop */
262  }
263  fprintf(fstats," --> Total eigenvalues found = %d\n",totcnt);
264  sprintf(path, "OUT/EigsOut_RLanN_%s", io.MatNam);
265  FILE *fmtout = fopen(path,"w");
266  if (fmtout) {
267  for (j=0; j<totcnt; j++)
268  fprintf(fmtout, "%.15e\n", alleigs[j]);
269  fclose(fmtout);
270  }
271  /*-------------------- free memory within this loop */
272  free(vinit);
273  free_coo(&Acoo);
274  free_csr(&Acsr);
275  free(sli);
276  free(counts);
277  free(alleigs);
278  if (fstats != stdout) fclose(fstats);
279  /*-------------------- end matrix loop */
280  }
281  //-------------------- free rest of memory
282  free(mu);
283  if( flog != stdout ) fclose ( flog );
284  fclose( fmat );
285  /*-------------------- finalize EVSL */
286  EVSLFinish();
287  return 0;
288 }
289 
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
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
int pow
Definition: struct.h:123
void StatsPrint(FILE *fstats)
Definition: stats.c:14
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
#define TRIV_SLICER
Definition: MMRLanN.c:13
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 main()
Definition: MMRLanN.c:18
void sort_double(int n, double *v, int *ind)
Definition: vect.c:92
int EVSLStart()
Initialize evslData.
Definition: evsl.c:27
void free_csr(csrMat *csr)
memory deallocation for csr matrix
Definition: spmat.c:91
#define HB
Definition: io.h:7
char Fname[MAX_LINE]
Definition: io.h:17
char MatNam[MaxNamLen]
Definition: io.h:18
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 get_matrix_info(FILE *fmat, io_t *pio)
Definition: io.c:22
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
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 nrows
Definition: struct.h:17
void StatsReset()
Definition: stats.c:79
void linspace(double a, double b, int num, double *arr)
Definition: vect.c:54
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
#define min(a, b)
Definition: MMRLanN.c:12
#define max(a, b)
Definition: MMRLanN.c:11
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
Definition: io.c:66
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
int n_intv
Definition: io.h:23
parameters for rational filter
Definition: struct.h:112