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