11 #define max(a, b) ((a) > (b) ? (a) : (b))
12 #define min(a, b) ((a) < (b) ? (a) : (b))
22 int n=0, sl, i, j, nev, totcnt;
23 double a, b, ecount, xintv[4];
40 FILE *flog = stdout, *fmat = NULL;
58 double *mu = malloc((Mdeg+1)*
sizeof(
double));
62 printf(
"-----------------------------------------\n");
63 printf(
"Note: You are using CXSparse for the direct solver. \n We recommend a more performance based direct solver for anything more than basic tests. \n SuiteSparse is supported with a makefile change. \n Using SuiteSparse can result in magnitudes faster times. \n\n");
64 printf(
"-----------------------------------------\n");
67 if( NULL == ( fmat = fopen(
"matfile",
"r" ) ) ) {
68 fprintf( flog,
"Can't open matfile...\n" );
73 if (NULL == fgets( line,
MAX_LINE, fmat )) {
74 fprintf( flog,
"error in reading matfile...\n" );
77 if( ( numat = atoi( line ) ) <= 0 ) {
78 fprintf( flog,
"Invalid count of matrices...\n" );
82 for(mat = 1; mat <= numat; mat++) {
84 fprintf(flog,
"Invalid format in matfile ...\n");
88 fprintf(flog,
"MATRIX: %s...\n", io.
MatNam);
94 if (stat(
"OUT", &st) == -1) {
99 strcpy( path,
"OUT/MMRLanN_");
101 fstats = fopen(path,
"w");
103 printf(
" failed in opening output file in OUT/\n");
106 fprintf(fstats,
"MATRIX: %s...\n", io.
MatNam);
107 fprintf(fstats,
"Partition the interval of interest [%f,%f] into %d slices\n",
109 counts = malloc(n_intv*
sizeof(
int));
110 sli = malloc( (n_intv+1)*
sizeof(
double));
115 fprintf(fstats,
"matrix read successfully\n");
119 fprintf(flog,
"read_coo error = %d\n", ierr);
126 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
132 alleigs = malloc(n*
sizeof(
double));
133 vinit = (
double *) malloc(n*
sizeof(
double));
136 ierr =
LanTrbounds(50, 200, 1e-8, vinit, 1, &lmin, &lmax, fstats);
137 fprintf(fstats,
"Step 0: Eigenvalue bounds for A: [%.15e, %.15e]\n",
141 fprintf(fstats,
" --> interval: a %9.3e b %9.3e \n",a, b);
143 xintv[0] = a; xintv[1] = b;
144 xintv[2] = lmin; xintv[3] = lmax;
150 ierr =
kpmdos(Mdeg, 0, nvec, xintv, mu, &ecount);
153 printf(
"kpmdos error %d\n", ierr);
156 fprintf(fstats,
" Time to build DOS (kpmdos) was : %10.2f \n",t);
157 fprintf(fstats,
"Step 1a: Estimated eig count in interval - %10.2e \n",
159 if (ecount < 0 || ecount > n) {
160 printf(
" e-count estimate is incorrect \n ");
170 ierr =
spslicer(sli, mu, Mdeg, xintv, n_intv, npnts);
172 printf(
"spslicer error %d\n", ierr);
176 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n",
178 fprintf(fstats,
"Step 1b: Slices found: \n");
179 for (j=0; j<n_intv;j++)
180 fprintf(fstats,
"[% 12.4e , % 12.4e]\n", sli[j],sli[j+1]);
184 printf(
"Total number of eigenvalues estimated = %d \n", (
int)(ecount));
185 nev = (int) (1 + ecount / ((
double) n_intv));
186 nev = (int)(fac*nev);
187 max_its =
max(4*nev,300); max_its =
min(max_its, n);
188 fprintf(fstats,
"Step 2: In each slice compute %d eigenvalues ... \n", nev);
191 for (sl=0; sl<n_intv; sl++) {
192 fprintf(fstats,
"======================================================\n");
193 double *lam, *Y, *res;
199 fprintf(fstats,
" subinterval: [% 12.4e , % 12.4e]\n", a, b);
201 fprintf(fstats,
" Max Lanczos steps %d\n", max_its);
207 printf(
"=== Compute the %dth subinterval: [%.4e, %.4e] out of %d ===\n",
219 void **solshiftdata = (
void **) malloc(rat.
num*
sizeof(
void *));
225 ierr =
RatLanNr(intv, max_its, tol, vinit, &rat, &nevOut, &lam,
228 printf(
"RatLanNr error %d\n", ierr);
233 ind = (
int *) malloc(nevOut*
sizeof(
int));
237 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
238 fprintf(fstats,
" Eigenvalues in [%f, %f]\n",
240 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
241 fprintf(fstats,
" Computed [%d out of %d estimated] ||Res|| ",
243 fprintf(fstats,
"\n");
244 for (i=0; i<nevOut; i++) {
245 fprintf(fstats,
" % .15e %.1e", lam[i], res[ind[i]]);
246 fprintf(fstats,
"\n");
248 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n");
249 memcpy(&alleigs[totcnt],lam,nevOut*
sizeof(
double));
263 fprintf(fstats,
" --> Total eigenvalues found = %d\n",totcnt);
264 sprintf(path,
"OUT/EigsOut_RLanN_%s", io.
MatNam);
265 FILE *fmtout = fopen(path,
"w");
267 for (j=0; j<totcnt; j++)
268 fprintf(fmtout,
"%.15e\n", alleigs[j]);
278 if (fstats != stdout) fclose(fstats);
283 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_coo(cooMat *coo)
memory deallocation for coo matrix
void ASIGMABSolDirect(int n, double *br, double *bi, double *xr, double *xz, void *data)
complex linear solver routine passed to evsl
void free_rat(ratparams *rat)
int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num, complex double *zk, void **data)
setup CXsparse solver for A - SIGMA B
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].
void FreeASIGMABSolDirect(int num, void **data)
free the data needed by CXSparse
void sort_double(int n, double *v, int *ind)
int EVSLStart()
Initialize evslData.
void free_csr(csrMat *csr)
memory deallocation for csr matrix
int RatLanNr(double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
Rational filtering Lanczos process [NON-restarted version].
int get_matrix_info(FILE *fmat, io_t *pio)
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...
void linspace(double a, double b, int num, double *arr)
void set_ratf_def(ratparams *rat)
Sets default values for ratparams struct.
int SetASigmaBSol(ratparams *rat, SolFuncC *func, SolFuncC allf, void **data)
Set the solve routine and the associated data for A-SIGMA*B if func == NULL, set all functions to be ...
Definitions used for direct solver interface.
int find_ratf(double *intv, ratparams *rat)
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
sparse matrix format: the coordinate (COO) format, 0-based
double evsl_timer()
evsl timer for mac
parameters for rational filter