EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
MMPLanR_omp.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 <omp.h>
7 #include "evsl.h"
8 #include "io.h"
9 
10 #define max(a, b) ((a) > (b) ? (a) : (b))
11 #define min(a, b) ((a) < (b) ? (a) : (b))
12 #define TRIV_SLICER 0
13 
14 int main () {
15  int ierr = 0;
16  /*--------------------------------------------------------------
17  * this tests the spectrum slicing idea for a generic matrix
18  * read in sparse matrix format
19  * Uses:
20  * Thick-restart Lanczos with pol. filtering
21  * this driver implements thread parallelism across slices/
22  *-------------------------------------------------------------*/
23  int n=0, sl, i, j, mlan, nev, totcnt;
24  //int nnz;
25  double a, b, ecount, xintv[4];
26  double lmin, lmax;
27  double *alleigs;
28  int n_intv; // number of subintervals (slices)
29  int npts; // number of integration points for eigenvalue count
30  //--------------------related to kpmdos
31  //int Mdeg = 180, nvec = 100;
32  /*-------------------- matrix A: coo format and csr format */
33  cooMat Acoo;
34  csrMat Acsr;
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  // double *lam=NULL, *Y=NULL, *res=NULL ;
41  /*-------------------- mim/max degree of polynomial, max Lanczos steps */
42  int max_its, Mdeg, nvec;
43  /*-------------------- IO */
44  FILE *flog = stdout, *fmat = NULL, *fstats = NULL;
45  io_t io;
46  int numat, mat;
47  char line[MAX_LINE];
48  /* initial vector: random */
49  double *vinit;
50  tol = 1e-8;
51  //-------------------- slicer parameters
52  Mdeg = 300;
53  nvec = 60;
54  /*-------------------- start EVSL */
55  EVSLStart();
56  //-------------------- interior eigensolver parameters
57  double *mu = malloc((Mdeg+1)*sizeof(double));
58  int *counts;
59  double *sli;
60  /*------------------ file "matfile" contains paths to matrices */
61  if( NULL == ( fmat = fopen( "matfile", "r" ) ) ) {
62  fprintf( flog, "Can't open matfile...\n" );
63  exit(2);
64  }
65  /*-------------------- read number of matrices ..*/
66  memset( line, 0, MAX_LINE );
67  if (NULL == fgets( line, MAX_LINE, fmat )) {
68  fprintf( flog, "error in reading matfile...\n" );
69  exit(2);
70  }
71  if( ( numat = atoi( line ) ) <= 0 ) {
72  fprintf( flog, "Invalid count of matrices...\n" );
73  exit(3);
74  }
75  /*-------------------- LOOP through matrices -*/
76  for(mat = 1; mat <= numat; mat++) {
77  if(get_matrix_info(fmat, &io) != 0) {
78  fprintf(flog, "Invalid format in matfile ...\n");
79  exit(5);
80  }
81  /*----------------input matrix and interval information -*/
82  fprintf(flog, "MATRIX: %s...\n", io.MatNam);
83  a = io.a;
84  b = io.b;
85  n_intv = io.n_intv;
86 
87  struct stat st = {0}; /* Make sure OUT directory exists */
88  if (stat("OUT", &st) == -1) {
89  mkdir("OUT", 0750);
90  }
91 
92  /*-------------------- path to write the output files*/
93  char path[1024];
94  strcpy( path, "OUT/MMPLanR_OMP_");
95  strcat( path, io.MatNam);
96  fstats = fopen(path,"w"); // write all the output to the file io.MatNam
97  if (!fstats) {
98  printf(" failed in opening output file in OUT/\n");
99  fstats = stdout;
100  }
101  fprintf(fstats, "MATRIX: %s...\n", io.MatNam);
102  fprintf(fstats,"Partition the interval of interest [%f,%f] into %d slices\n", a,b,n_intv);
103  counts = malloc(n_intv*sizeof(int));
104  sli = malloc((n_intv+1)*sizeof(double));
105  /*-------------------- Read matrix - case: COO/MatrixMarket formats */
106  if (io.Fmt > HB) {
107  ierr =read_coo_MM(io.Fname, 1, 0, &Acoo);
108  if (ierr == 0) {
109  fprintf(fstats,"matrix read successfully\n");
110  //nnz = Acoo.nnz;
111  n = Acoo.nrows;
112  }
113  else {
114  fprintf(flog, "read_coo error = %d\n", ierr);
115  exit(6);
116  }
117  /*-------------------- conversion from COO to CSR format */
118  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
119  }
120  if (io.Fmt == HB) {
121  fprintf(flog, "HB FORMAT not supported (yet) * \n");
122  exit(7);
123  }
124  /*-------------------- set the left-hand side matrix A */
125  SetAMatrix(&Acsr);
126  /*-------------------- define ChebLanTr parameters */
127  alleigs = (double *) malloc(n*sizeof(double));
128  vinit = (double *) malloc(n*sizeof(double));
129  rand_double(n, vinit);
130  /*-------------------- get lambda_min lambda_max estimates */
131  ierr = LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
132  fprintf(fstats, "Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n", lmin, lmax);
133  /*-------------------- define [a b] now so we can get estimates now
134  on number of eigenvalues in [a b] from kpmdos */
135  fprintf(fstats," --> interval: a %9.3e b %9.3e \n",a, b);
136  /*-------------------- define kpmdos parameters */
137  xintv[0] = a; xintv[1] = b;
138  xintv[2] = lmin; xintv[3] = lmax;
139  //-------------------- kpmdos or triv_slicer
140  if (TRIV_SLICER) {
141  linspace(a, b, n_intv+1, sli);
142  } else {
143  double t = evsl_timer();
144  ierr = kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
145  t = evsl_timer() - t;
146  if (ierr) {
147  printf("kpmdos error %d\n", ierr);
148  return 1;
149  }
150  fprintf(fstats, " Time to build DOS (kpmdos) was : %10.2f \n",t);
151  fprintf(fstats, "Step 1a: Estimated eig count in interval - %.15e \n",ecount);
152  if ((ecount <0) | (ecount > n)) {
153  printf(" e-count estimate is incorrect \n ");
154  exit(7);
155  }
156  /*-------------------- error in ecount */
157  /*
158  if (ecount < 1 || ecount > n)
159  exit(6);
160  */
161  /*-------------------- define slicer parameters */
162  npts = 10 * ecount;
163  ierr = spslicer(sli, mu, Mdeg, xintv, n_intv, npts);
164  if (ierr) {
165  printf("spslicer error %d\n", ierr);
166  return 1;
167  }
168  }
169  fprintf(fstats,"DOS parameters: Mdeg = %d, nvec = %d, npts = %d\n",
170  Mdeg, nvec, npts);
171  fprintf(fstats, "Step 1b: Slices found: \n");
172  for (j=0; j<n_intv;j++) {
173  fprintf(fstats, " %2d: [%.15e , %.15e]\n", j+1, sli[j],sli[j+1]);
174  }
175  //-------------------- # eigs per slice
176  //-------------------- approximate number of eigenvalues wanted
177  double fac = 1.2; // this factor insures that # of eigs per slice is slightly overestimated
178  nev = (int) (1 + ecount / ((double) n_intv)); // # eigs per slice
179  nev = (int)(fac*nev); // want an overestimate of ev_int
180  fprintf(fstats,"Step 2: In each slice compute %d eigenvalues ... \n", nev);
181  /*-------------------- MAIN intv LOOP: for each sclice Do: */
182  double tsolve = evsl_timer();
183  totcnt = 0;
184  mlan = max(4*nev,300); mlan = min(n, mlan);
185  max_its = 3*mlan;
186  // shared array for storing eigenvalues and residuals for each slice
187  double *lam_global, *res_global;
188  lam_global = (double *) malloc(nev*n_intv*sizeof(double));
189  res_global = (double *) malloc(nev*n_intv*sizeof(double));
190 
191 #pragma omp parallel for private(sl)
192  for (sl=0; sl<n_intv; sl++) {
193  int nevOut, j;
194  double *lam, *Y, *res;
195  double *res_sorted;
196  int *ind;
197  double intv[4];
198  polparams pol;
199  /*
200  int tmp1 = omp_get_thread_num()+1;
201  int tmp2 = omp_get_num_threads();
202  #pragma omp critical
203  fprintf(flog,"Thread %d out of %d: \n", tmp1, tmp2);
204  */
205  //-------------------- ChebLanTr
206  intv[0] = sli[sl];
207  intv[1] = sli[sl+1];
208  intv[2] = lmin;
209  intv[3] = lmax;
210 
211  set_pol_def(&pol);
212  find_pol(intv, &pol);
213 
214  fprintf(fstats, " polynomial [type = %d], deg %d, bar %e gam %e\n",
215  pol.type, pol.deg, pol.bar, pol.gam);
216 
217  ierr = ChebLanTr(mlan, nev, intv, max_its, tol, vinit,
218  &pol, &nevOut, &lam, &Y, &res, NULL);
219  if (ierr) {
220  printf("ChebLanTr error %d\n", ierr);
221  }
222  /* sort the eigenvals: ascending order
223  * ind: keep the orginal indices */
224  ind = (int *) malloc(nevOut*sizeof(int));
225  sort_double(nevOut, lam, ind);
226  // sort residuals accordingly
227  res_sorted = (double *) malloc(nevOut*sizeof(double));
228  for (j=0; j<nevOut; j++){
229  res_sorted[j] = res[ind[j]];
230  }
231 
232  /* copy result to global arrays */
233  memcpy(&lam_global[sl*nev],lam,nevOut*sizeof(double));
234  memcpy(&res_global[sl*nev],res_sorted,nevOut*sizeof(double));
235  counts[sl] = nevOut;
236  //-------------------- free memory within this loop
237  if (lam) free(lam);
238  if (Y) free(Y);
239  if (res) free(res);
240  free(res_sorted);
241  free(ind);
242  free_pol(&pol);
243  /*-------------------- end slice loop */
244  } // end of parallel for-loop
245 
246  tsolve = evsl_timer() - tsolve;
247 
248  // Now output result and move computed eigenvalues to alleigs
249  for (sl=0; sl<n_intv; sl++) {
250  //--------------------
251  a = sli[sl];
252  b = sli[sl+1];
253 
254  fprintf(fstats,"======================================================\n");
255  fprintf(fstats, " subinterval: [% 12.4e , % 12.4e ]\n", a, b);
256  fprintf(fstats, "Thick Restarted Lanczos with dimension %d\n", mlan);
257  fprintf(fstats, "Max Lanczos steps %d\n", max_its);
258 
259  /* print eigenvalues */
260  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
261  fprintf(fstats, " Eigenvalues in [%f, %f]\n",a,b);
262  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
263  fprintf(fstats, " Computed [%d out of %d estimated] ||Res|| ", counts[sl], nev);
264  fprintf(fstats, "\n");
265  for (i=0; i<counts[sl]; i++) {
266  fprintf(fstats, " % .15e %.1e",
267  lam_global[sl*nev+i], res_global[sl*nev+i]);
268  fprintf(fstats,"\n");
269  }
270  fprintf(fstats, "- - - - - - - - - - - - - - - - - - - - - - - - - -\n");
271 
272  memcpy(&alleigs[totcnt],&lam_global[sl*nev],counts[sl]*sizeof(double));
273  totcnt += counts[sl];
274  }
275  fprintf(fstats, "Solution time : %.2f\n", tsolve);
276  fprintf(fstats," --> Total eigenvalues found = %d\n",totcnt);
277  //FILE *fmtout = fopen("EigsOut","w");
278  //for (j=0; j<totcnt; j++)
279  // fprintf(fmtout, "%22.15e\n", alleigs[j]);
280  //fclose(fmtout);
281  /*-------------------- free memory */
282  free(counts);
283  free(sli);
284  free(lam_global);
285  free(res_global);
286  free(vinit);
287  free_coo(&Acoo);
288  free_csr(&Acsr);
289  free(alleigs);
290  if (fstats != stdout) fclose(fstats);
291  /*-------------------- end matrix loop */
292  }
293  //-------------------- free rest of memory
294  free(mu);
295  if( flog != stdout ) fclose ( flog );
296  fclose( fmat );
297  /*-------------------- finalize EVSL */
298  EVSLFinish();
299  return 0;
300 }
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 main()
Definition: MMPLanR_omp.c:14
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 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
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
#define TRIV_SLICER
Definition: MMPLanR_omp.c:12
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 max(a, b)
Definition: MMPLanR_omp.c:10
#define HB
Definition: io.h:7
char Fname[MAX_LINE]
Definition: io.h:17
char MatNam[MaxNamLen]
Definition: io.h:18
void set_pol_def(polparams *pol)
set default values for polparams struct.
Definition: chebpoly.c:18
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.
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 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 evsl_timer()
evsl timer for mac
Definition: mactime.c:14
double bar
Definition: struct.h:65
#define min(a, b)
Definition: MMPLanR_omp.c:11
int n_intv
Definition: io.h:23