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;
21 double a, b, ecount, xintv[4];
39 int max_its, Mdeg, nvec;
41 FILE *flog = stdout, *fmat = NULL, *fstats = NULL;
56 double *mu = malloc((Mdeg+1)*
sizeof(
double));
60 if ( NULL == ( fmat = fopen(
"matfile",
"r" ) ) ) {
61 fprintf( flog,
"Can't open matfile...\n" );
66 if (NULL == fgets( line,
MAX_LINE, fmat )) {
67 fprintf( flog,
"error in reading matfile...\n" );
70 if( ( numat = atoi( line ) ) <= 0 ) {
71 fprintf( flog,
"Invalid count of matrices...\n" );
75 for(mat = 1; mat <= numat; mat++) {
77 fprintf(flog,
"Invalid format in matfile ...\n");
81 fprintf(flog,
"MATRIX: %s...\n", io.
MatNam);
87 if (stat(
"OUT", &st) == -1) {
93 strcpy( path,
"OUT/MMPLanR_");
95 fstats = fopen(path,
"w");
97 printf(
" failed in opening output file in OUT/\n");
100 fprintf(fstats,
"MATRIX: %s...\n", io.
MatNam);
101 fprintf(fstats,
"Partition the interval of interest [%f,%f] into %d slices\n", a,b,n_intv);
102 counts = malloc(n_intv*
sizeof(
int));
103 sli = malloc((n_intv+1)*
sizeof(
double));
108 fprintf(fstats,
"matrix read successfully\n");
113 fprintf(flog,
"read_coo error = %d\n", ierr);
120 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
126 alleigs = (
double *) malloc(n*
sizeof(
double));
127 vinit = (
double *) malloc(n*
sizeof(
double));
130 ierr =
LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
131 fprintf(fstats,
"Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n", lmin, lmax);
134 fprintf(fstats,
" --> interval: a %9.3e b %9.3e \n",a, b);
136 xintv[0] = a; xintv[1] = b;
137 xintv[2] = lmin; xintv[3] = lmax;
143 ierr =
kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
146 printf(
"kpmdos error %d\n", ierr);
149 fprintf(fstats,
" Time to build DOS (kpmdos) was : %10.2f \n",t);
150 fprintf(fstats,
"Step 1a: Estimated eig count in interval - %.15e \n",ecount);
151 if ((ecount <0) | (ecount > n)) {
152 printf(
" e-count estimate is incorrect \n ");
156 if (ecount < 1 || ecount > n)
160 ierr =
spslicer(sli, mu, Mdeg, xintv, n_intv, npts);
162 printf(
"spslicer error %d\n", ierr);
166 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npts = %d\n",
168 fprintf(fstats,
"Step 1b: Slices found: \n");
169 for (j=0; j<n_intv;j++) {
170 fprintf(fstats,
" %2d: [%.15e , %.15e]\n", j+1, sli[j],sli[j+1]);
175 nev = (int) (1 + ecount / ((
double) n_intv));
176 nev = (int)(fac*nev);
177 fprintf(fstats,
"Step 2: In each slice compute %d eigenvalues ... \n", nev);
180 for (sl=0; sl<n_intv; sl++) {
181 fprintf(fstats,
"======================================================\n");
182 double *lam, *Y, *res;
187 fprintf(fstats,
" subinterval: [% 12.4e , % 12.4e]\n", a, b);
189 mlan =
max(4*nev,300); mlan =
min(mlan, n);
202 fprintf(fstats,
" polynomial [type = %d], deg %d, bar %e gam %e\n",
205 ierr =
ChebLanTr(mlan, nev, xintv, max_its, tol, vinit,
206 &pol, &nevOut, &lam, &Y, &res, fstats);
208 printf(
"ChebLanTr error %d\n", ierr);
213 ind = (
int *) malloc(nevOut*
sizeof(
int));
217 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
218 fprintf(fstats,
" Eigenvalues in [%f, %f]\n",a,b);
219 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
220 fprintf(fstats,
" Computed [%d out of %d estimated] ||Res|| ", nevOut, nev);
221 fprintf(fstats,
"\n");
222 for (i=0; i<nevOut; i++) {
223 fprintf(fstats,
" % .15e %.1e", lam[i], res[ind[i]]);
224 fprintf(fstats,
"\n");
226 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
227 memcpy(&alleigs[totcnt],lam,nevOut*
sizeof(
double));
238 fprintf(fstats,
" --> Total eigenvalues found = %d\n",totcnt);
239 sprintf(path,
"OUT/EigsOut_PLanR_%s", io.
MatNam);
240 FILE *fmtout = fopen(path,
"w");
242 for (j=0; j<totcnt; j++)
243 fprintf(fmtout,
"%.15e\n", alleigs[j]);
253 if (fstats != stdout) fclose(fstats);
258 if( flog != stdout ) fclose ( flog );
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...
void free_pol(polparams *pol)
void free_coo(cooMat *coo)
memory deallocation for coo matrix
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].
int SetAMatrix(csrMat *A)
Set the matrix A.
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)
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.
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 get_matrix_info(FILE *fmat, io_t *pio)
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...
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
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