9 #define max(a, b) ((a) > (b) ? (a) : (b))
10 #define min(a, b) ((a) < (b) ? (a) : (b))
22 int n=0, sl, i, j, nev, totcnt;
24 double a, b, ecount, xintv[4];
43 int max_its, Mdeg, nvec;
45 FILE *flog = stdout, *fmat = NULL, *fstats = NULL;
71 double *mu = malloc((Mdeg+1)*
sizeof(
double));
77 if( NULL == ( fmat = fopen(
"matfile",
"r" ) ) ) {
78 fprintf( flog,
"Can't open matfile...\n" );
83 if (NULL == fgets( line,
MAX_LINE, fmat )) {
84 fprintf( flog,
"error in reading matfile...\n" );
87 if( ( numat = atoi( line ) ) <= 0 ) {
88 fprintf( flog,
"Invalid count of matrices...\n" );
92 for(mat = 1; mat <= numat; mat++) {
94 fprintf(flog,
"Invalid format in matfile ...\n");
98 fprintf(flog,
"MATRIX: %s...\n", io.
MatNam);
103 struct stat st = {0};
104 if (stat(
"OUT", &st) == -1) {
109 strcpy( path,
"OUT/MMPSI_");
112 fstats = fopen(path,
"w");
114 printf(
" failed in opening output file in OUT/\n");
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));
125 fprintf(fstats,
"matrix read successfully\n");
130 fprintf(flog,
"read_coo error = %d\n", ierr);
137 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
143 alleigs = malloc(n*
sizeof(
double));
144 vinit = (
double *) malloc(n*
sizeof(
double));
147 ierr =
LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
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);
152 xintv[0] = a; xintv[1] = b;
153 xintv[2] = lmin; xintv[3] = lmax;
155 ierr =
kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
157 printf(
"kpmdos error %d\n", ierr);
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 ");
176 ierr =
spslicer(sli, mu, Mdeg, xintv, n_intv, npnts);
178 printf(
"spslicer error %d\n", ierr);
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]);
190 nev = (int) (1 + ecount / ((
double) n_intv));
191 nev = (int)(fac*nev);
192 fprintf(fstats,
"Step 2: In each slice compute %d eigenvalues ... \n", nev);
195 for (sl =0; sl<n_intv; sl++){
196 fprintf(fstats,
"======================================================\n");
197 double *lam, *Y, *res;
202 fprintf(fstats,
" subinterval: [% 12.4e , % 12.4e]\n", a, b);
204 fprintf(fstats,
"Filtered subspace iteration with block size %d\n", nev);
205 fprintf(fstats,
"Max steps %d\n", max_its);
210 V0 = (
double *) malloc(n*nev*
sizeof(
double));
211 for (i=0; i<n*nev; i++){
212 V0[i] = rand() / ((double)RAND_MAX);
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);
224 printf(
"ChebSI error %d\n", ierr);
229 ind = (
int*) malloc(nevOut*
sizeof(
int));
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");
242 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
244 fprintf(fstats,
"======================================================\n");
246 memcpy(&alleigs[totcnt],lam,nevOut*
sizeof(
double));
258 fprintf(fstats,
" --> Total eigenvalues found = %d\n",totcnt);
259 sprintf(path,
"OUT/EigsOut_PSI_%s",io.
MatNam);
260 FILE *fmtout = fopen(path,
"w");
262 for (j=0; j<totcnt; j++)
263 fprintf(fmtout,
"%.15e\n", alleigs[j]);
273 if (fstats != stdout) fclose(fstats);
277 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 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.
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.
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 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)
int get_matrix_info(FILE *fmat, io_t *pio)
sparse matrix format: the coordinate (COO) format, 0-based