9 #define max(a, b) ((a) > (b) ? (a) : (b))
10 #define min(a, b) ((a) < (b) ? (a) : (b))
19 int n=0, sl, i, j, mlan, nev, totcnt;
20 double a, b, ecount, xintv[4];
40 FILE *flog = stdout, *fmat = NULL, *fstats = NULL;
45 double *vinit, *xdos, *ydos;
55 if( NULL == ( fmat = fopen(
"matfile",
"r" ) ) ) {
56 fprintf( flog,
"Can't open matfile...\n" );
61 if (NULL == fgets( line,
MAX_LINE, fmat )) {
62 fprintf( flog,
"error in reading matfile...\n" );
65 if( ( numat = atoi( line ) ) <= 0 ) {
66 fprintf( flog,
"Invalid count of matrices...\n" );
70 for(mat = 1; mat <= numat; mat++) {
72 fprintf(flog,
"Invalid format in matfile ...\n");
76 fprintf(flog,
"MATRIX: %s...\n", io.
MatNam);
82 if (stat(
"OUT", &st) == -1) {
88 strcpy( path,
"OUT/MMPLanN_");
90 fstats = fopen(path,
"w");
92 printf(
" failed in opening output file in OUT/\n");
95 fprintf(fstats,
"MATRIX: %s...\n", io.
MatNam);
96 fprintf(fstats,
"Partition the interval of interest [%f,%f] into %d slice (s)\n", a,b,n_intv);
97 counts = malloc(n_intv*
sizeof(
int));
98 sli = malloc((n_intv+1)*
sizeof(
double));
103 fprintf(fstats,
"matrix read successfully\n");
107 fprintf(flog,
"read_coo error = %d\n", ierr);
114 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
120 alleigs = (
double *) malloc(n*
sizeof(
double));
121 vinit = (
double *) malloc(n*
sizeof(
double));
124 ierr =
LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
125 fprintf(fstats,
"Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n", lmin, lmax);
128 fprintf(fstats,
" --> interval: a %9.3e b %9.3e \n",a, b);
130 xintv[0] = a; xintv[1] = b;
131 xintv[2] = lmin; xintv[3] = lmax;
138 xdos = (
double *) malloc((npts)*
sizeof(double));
139 ydos = (
double *) malloc((npts)*
sizeof(double));
142 LanDos(nvec, Mdeg, npts, xdos, ydos, &ecount, xintv);
146 fprintf(fstats,
" Time to build DOS (Landos) was : %10.2f \n",t);
147 fprintf(fstats,
" estimated eig count in interval: %.15e \n",ecount);
149 if ((ecount <0) | (ecount > n)) {
150 printf(
" e-count estimate is incorrect \n ");
154 if ((ecount <1) | (ecount>n))
158 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n",
160 spslicer2(xdos, ydos, n_intv, npts, sli);
165 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npts = %d\n",
167 fprintf(fstats,
"Step 1b: Slices found: \n");
168 for (j=0; j<n_intv;j++) {
169 fprintf(fstats,
" %2d: [%.15e , %.15e]\n", j+1, sli[j],sli[j+1]);
174 nev = (int) (1 + ecount / ((
double) n_intv));
175 nev = (int)(fac*nev);
176 fprintf(fstats,
"Step 2: In each slice compute %d eigenvalues ... \n", nev);
179 for (sl =0; sl<n_intv; sl++) {
180 fprintf(fstats,
"======================================================\n");
181 double *lam, *Y, *res;
187 fprintf(fstats,
" subinterval: [% 12.4e , % 12.4e]\n", a, b);
189 mlan =
max(5*nev,300); mlan =
min(mlan, n);
191 fprintf(fstats,
"Non-Restarted Lanczos with dimension %d\n", mlan);
203 fprintf(fstats,
" polynomial [type = %d], deg %d, bar %e gam %e\n",
206 ierr =
ChebLanNr(xintv, mlan, tol, vinit, &pol, &nevOut, &lam, &Y,
209 printf(
"ChebLanNr error %d\n", ierr);
214 ind = (
int*) malloc(nevOut*
sizeof(
int));
218 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
219 fprintf(fstats,
" Eigenvalues in [%f, %f]\n",a,b);
220 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
221 fprintf(fstats,
" Computed [%d out of %d estimated] ||Res|| ", nevOut, nev);
222 fprintf(fstats,
"\n");
223 for (i=0; i<nevOut; i++) {
224 fprintf(fstats,
" % .15e %.1e", lam[i], res[ind[i]]);
225 fprintf(fstats,
"\n");
227 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
228 memcpy(&alleigs[totcnt],lam,nevOut*
sizeof(
double));
240 fprintf(fstats,
" --> Total eigenvalues found = %d\n",totcnt);
241 sprintf(path,
"OUT/EigsOut_PLanN_%s", io.
MatNam);
242 FILE *fmtout = fopen(path,
"w");
244 for (j=0; j<totcnt; j++)
245 fprintf(fmtout,
"%.15e\n", alleigs[j]);
255 if (fstats != stdout) fclose(fstats);
259 if( flog != stdout ) fclose ( flog );
void free_pol(polparams *pol)
void free_coo(cooMat *coo)
memory deallocation for coo matrix
int get_matrix_info(FILE *fmat, io_t *pio)
int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos, double *neig, const double *const intv)
int SetAMatrix(csrMat *A)
Set the matrix A.
void StatsPrint(FILE *fstats)
void rand_double(int n, double *v)
int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr)
convert coo to csr
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].
int find_pol(double *intv, polparams *pol)
Sets the values in pol.
void sort_double(int n, double *v, int *ind)
void spslicer2(double *xi, double *yi, int n_int, int npts, double *sli)
int EVSLStart()
Initialize evslData.
void free_csr(csrMat *csr)
memory deallocation for csr matrix
void set_pol_def(polparams *pol)
set default values for polparams struct.
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
sparse matrix format: the compressed sparse row (CSR) format, 0-based
This file contains function prototypes and constant definitions for EVSL.
int EVSLFinish()
Finish EVSL.
int ChebLanNr(double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
parameters for polynomial filter
void linspace(double a, double b, int num, double *arr)
sparse matrix format: the coordinate (COO) format, 0-based
double evsl_timer()
evsl timer for mac