EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
LanDos.c
Go to the documentation of this file.
1 #include <errno.h>
2 #include <fcntl.h>
3 #include <math.h>
4 #include <stdio.h>
5 #include <stdlib.h>
6 #include <string.h>
7 #include <sys/stat.h>
8 #include <sys/types.h>
9 #include <sys/utsname.h>
10 #include <time.h>
11 #include <unistd.h>
12 #include "evsl.h"
13 #include "io.h"
14 
15 /*-------------------- protos */
16 int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv);
17 int findarg(const char *argname, ARG_TYPE type, void *val, int argc,
18  char **argv);
19 
20 #define max(a, b) ((a) > (b) ? (a) : (b))
21 #define min(a, b) ((a) < (b) ? (a) : (b))
22 
31 int readVec(const char *filename, int *npts, double **vec) {
32  int i;
33  FILE *ifp = fopen(filename, "r");
34  if (ifp) {
35  fscanf(ifp, "%i", npts);
36  *vec = (double *)malloc(sizeof(double) * *npts);
37  for (i = 0; i < (*npts); i++) {
38  fscanf(ifp, "%lf", (&(*vec)[i]));
39  }
40  }
41  fclose(ifp);
42  return 0;
43 }
44 
52 int main(int argc, char *argv[]) {
53  srand(time(NULL));
54  const int msteps = 40; /* Steps to perform */
55  const int npts = 200; /* Number of points */
56  const int nvec = 100; /* Number of random vectors to generate */
57  double intv[4];
58 
59  cooMat Acoo;
60  csrMat Acsr;
61 
62  /*-------------------- IO */
63  FILE *flog = stdout, *fmat = NULL, *fstats = NULL;
64  io_t io;
65  int numat, mat, ierr, n=0, graph_exact_dos = 0;
66  char line[MAX_LINE];
67 
68  findarg("graph_exact_dos", INT, &graph_exact_dos, argc, argv);
69  /* initial vector: random */
70  /*-------------------- start EVSL */
71  EVSLStart();
72  //-------------------- interior eigensolver parameters
73  /*------------------ file "matfile" contains paths to matrices */
74  if (NULL == (fmat = fopen("matfile", "r"))) {
75  fprintf(flog, "Can't open matfile...\n");
76  exit(2);
77  }
78  /*-------------------- read number of matrices ..*/
79  memset(line, 0, MAX_LINE);
80  if (NULL == fgets(line, MAX_LINE, fmat)) {
81  fprintf(flog, "error in reading matfile...\n");
82  exit(2);
83  }
84  if ((numat = atoi(line)) <= 0) {
85  fprintf(flog, "Invalid count of matrices...\n");
86  exit(3);
87  }
88  /*-------------------- LOOP through matrices -*/
89  for (mat = 1; mat <= numat; mat++) {
90  if (get_matrix_info(fmat, &io) != 0) {
91  fprintf(flog, "Invalid format in matfile ...\n");
92  exit(5);
93  }
94  /*----------------input matrix and interval information -*/
95  fprintf(flog, "MATRIX: %s...\n", io.MatNam);
96 
97  /*-------------------- path to write the output files*/
98  char path[1024];
99  strcpy(path, "OUT/LanDos_");
100  strcat(path, io.MatNam);
101  strcat(path, ".log");
102  fstats = fopen(path, "w"); // write all the output to the file io.MatNam
103  if (!fstats) {
104  printf(" failed in opening output file in OUT/\n");
105  fstats = stdout;
106  }
107  fprintf(fstats, "MATRIX: %s...\n", io.MatNam);
108  /*-------------------- Read matrix - case: COO/MatrixMarket formats */
109  if (io.Fmt > HB) {
110  ierr = read_coo_MM(io.Fname, 1, 0, &Acoo);
111  if (ierr == 0) {
112  fprintf(fstats, "matrix read successfully\n");
113  // nnz = Acoo.nnz;
114  n = Acoo.nrows;
115  } else {
116  fprintf(flog, "read_coo error = %d\n", ierr);
117  exit(6);
118  }
119  /*-------------------- conversion from COO to CSR format */
120  ierr = cooMat_to_csrMat(0, &Acoo, &Acsr);
121  }
122  if (io.Fmt == HB) {
123  fprintf(flog, "HB FORMAT not supported (yet) * \n");
124  exit(7);
125  }
126 
127  /*-------------------- Define some constants to test with */
128  /*-------------------- Intervals of interest
129  intv[0] and intv[1] are the input interval of interest [a. b]
130  intv[2] and intv[3] are the smallest and largest eigenvalues of (A,B)
131  */
132  int i, ret;
133  double neig; /* Number eigenvalues */
134  /*-------------------- exact histogram and computed DOS */
135  double *xHist = NULL;
136  double *yHist = NULL;
137  if (graph_exact_dos) {
138  xHist = (double *)malloc(npts * sizeof(double)); /*Exact DOS x values */
139  yHist = (double *)malloc(npts * sizeof(double)); /* Exact DOS y values */
140  }
141  double *xdos =
142  (double *)malloc(npts * sizeof(double)); /* Calculated DOS x vals */
143  double *ydos =
144  (double *)malloc(npts * sizeof(double)); /* Calculated DOS y */
145 
146  SetStdEig();
147  SetAMatrix(&Acsr);
148 
149  double *vinit = (double *)malloc(n * sizeof(double));
150  rand_double(n, vinit);
151  double lmin = 0.0, lmax = 0.0;
152  ierr = LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
153 
154  intv[0] = lmin;
155  intv[1] = lmax;
156  /*----------------- get the bounds for (A, B) ---------*/
157  intv[2] = lmin;
158  intv[3] = lmax;
159 
160  double *ev = NULL;
161  int numev;
162  if (graph_exact_dos) {
163  readVec("NM1AB_eigenvalues.dat", &numev, &ev);
164  }
165 
166  /* Calculate computed DOS */
167  ret = LanDos(nvec, msteps, npts, xdos, ydos, &neig, intv);
168  fprintf(stdout, " LanDos ret %d \n", ret);
169 
170  /* Calculate the exact DOS */
171  if (graph_exact_dos) {
172  ret = exDOS(Acoo.vv, Acoo.ncols, npts, xHist, yHist, intv);
173  fprintf(stdout, " exDOS ret %d \n", ret);
174  }
175  free_coo(&Acoo);
176  free_csr(&Acsr);
177  free(vinit);
178  /* --------------------Make OUT dir if it doesn't exist */
179  struct stat st = {0};
180  if (stat("OUT", &st) == -1) {
181  mkdir("OUT", 0750);
182  }
183 
184  /*-------------------- Write to output files */
185  char computed_path[1024];
186  strcpy(computed_path, "OUT/LanDos_Approx_DOS_");
187  strcat(computed_path, io.MatNam);
188  FILE *ofp = fopen(computed_path, "w");
189  for (i = 0; i < npts; i++) {
190  fprintf(ofp, " %10.4f %10.4f\n", xdos[i], ydos[i]);
191  }
192  fclose(ofp);
193 
194  if (graph_exact_dos) {
195  /*-------------------- save exact DOS */
196  strcpy(path, "OUT/LanDosG_Exact_DOS_");
197  strcat(path, io.MatNam);
198  ofp = fopen(path, "w");
199  for (i = 0; i < npts; i++)
200  fprintf(ofp, " %10.4f %10.4f\n", xHist[i], yHist[i]);
201  fclose(ofp);
202  }
203  printf("The data output is located in OUT/ \n");
204  struct utsname buffer;
205  errno = 0;
206  if (uname(&buffer) != 0) {
207  perror("uname");
208  exit(EXIT_FAILURE);
209  }
210  char command[1024];
211  strcpy(command, "gnuplot ");
212  strcat(command, " -e \"filename='");
213  strcat(command, computed_path);
214 
215  if (graph_exact_dos) {
216  strcat(command, "';exact_dos='");
217  strcat(command, path);
218  strcat(command, "'\" tester_ex.gnuplot");
219  ierr = system(command);
220  } else {
221  strcat(command, "'\" tester.gnuplot");
222  ierr = system(command);
223  }
224 
225  if (ierr) {
226  fprintf(stderr,
227  "Error using 'gnuplot < tester.gnuplot', \n"
228  "postscript plot could not be generated \n");
229  } else {
230  printf("A postscript graph has been placed in %s%s\n", computed_path,
231  ".eps");
232  /*-------------------- and gv for visualizing / */
233  if (!strcmp(buffer.sysname, "Linux")) {
234  strcpy(command, "gv ");
235  strcat(command, computed_path);
236  strcat(command, ".eps &");
237  ierr = system(command);
238  if (ierr) {
239  fprintf(stderr, "Error using 'gv %s' \n", command);
240  printf(
241  "To view the postscript graph use a postcript viewer such as "
242  "gv \n");
243  }
244  } else {
245  printf(
246  "To view the postscript graph use a postcript viewer such as "
247  "gv \n");
248  }
249  }
250 
251  if (graph_exact_dos) {
252  free(xHist);
253  free(yHist);
254  }
255  free(xdos);
256  free(ydos);
257  if(ev) {
258  free(ev);
259  }
260  fclose(fstats);
261  }
262  fclose(fmat);
263  EVSLFinish();
264  return 0;
265 }
void free_coo(cooMat *coo)
memory deallocation for coo matrix
Definition: spmat.c:126
int main(int argc, char *argv[])
Definition: LanDos.c:52
#define MAX_LINE
Definition: io.h:5
int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos, double *neig, const double *const intv)
Definition: landos.c:40
ARG_TYPE
Definition: io.h:33
double * vv
Definition: struct.h:22
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
Definition: io.h:34
int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv)
Definition: exDOS.c:22
int EVSLStart()
Initialize evslData.
Definition: evsl.c:27
int findarg(const char *argname, ARG_TYPE type, void *val, int argc, char **argv)
Definition: io.c:171
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 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
int ncols
Definition: struct.h:17
This file contains function prototypes and constant definitions for EVSL.
Definition: io.h:14
int EVSLFinish()
Finish EVSL.
Definition: evsl.c:48
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
Definition: io.c:66
int nrows
Definition: struct.h:17
int readVec(const char *filename, int *npts, double **vec)
Definition: LanDos.c:31
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