10 #define max(a, b) ((a) > (b) ? (a) : (b))
11 #define min(a, b) ((a) < (b) ? (a) : (b))
24 int n = 0, i, j, npts, nslices, nvec, Mdeg, nev, mlan, ev_int, sl,
27 double a, b, lmax, lmin, ecount, tol, *sli, *mu;
40 mu = malloc((Mdeg + 1) *
sizeof(
double));
41 FILE *flog = stdout, *fmat = NULL;
49 printf(
"-----------------------------------------\n");
50 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");
51 printf(
"-----------------------------------------\n");
59 if (NULL == (fmat = fopen(
"matfile",
"r"))) {
60 fprintf(flog,
"Can't open matfile...\n");
65 if (NULL == fgets(line,
MAX_LINE, fmat)) {
66 fprintf(flog,
"error in reading matfile...\n");
69 if ((numat = atoi(line)) <= 0) {
70 fprintf(flog,
"Invalid count of matrices...\n");
73 for (mat = 1; mat <= numat; mat++) {
75 fprintf(flog,
"Invalid format in matfile ...\n");
79 fprintf(flog,
"MATRIX A: %s...\n", io.
MatNam1);
80 fprintf(flog,
"MATRIX B: %s...\n", io.
MatNam2);
86 if (stat(
"OUT", &st) == -1) {
91 strcpy(path,
"OUT/KPM_MMPLanN_");
93 fstats = fopen(path,
"w");
95 printf(
" failed in opening output file in OUT/\n");
98 fprintf(fstats,
"MATRIX A: %s...\n", io.
MatNam1);
99 fprintf(fstats,
"MATRIX B: %s...\n", io.
MatNam2);
100 fprintf(fstats,
"Partition the interval of interest [%f,%f] into %d slices\n",
102 counts = malloc(nslices *
sizeof(
int));
103 sli = malloc((nslices + 1) *
sizeof(
double));
108 fprintf(fstats,
"matrix read successfully\n");
114 fprintf(flog,
"read_coo error for A = %d\n", ierr);
119 fprintf(fstats,
"matrix read successfully\n");
120 if (Bcoo.
nrows != n) {
128 fprintf(flog,
"read_coo error for B = %d\n", ierr);
136 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
139 alleigs = malloc(n *
sizeof(
double));
154 vinit = (
double *)malloc(n *
sizeof(
double));
156 ierr =
LanTrbounds(50, 200, 1e-12, vinit, 1, &lmin, &lmax, fstats);
157 fprintf(fstats,
"Step 0: Eigenvalue bound s for B^{-1}*A: [%.15e, %.15e]\n",
168 ierr =
kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
171 printf(
"kpmdos error %d\n", ierr);
174 fprintf(fstats,
" Time to build DOS (kpmdos) was : %10.2f \n", t);
175 fprintf(fstats,
" estimated eig count in interval: %.15e \n", ecount);
178 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n", Mdeg,
180 ierr =
spslicer(sli, mu, Mdeg, xintv, nslices, npts);
182 printf(
"spslicer error %d\n", ierr);
185 printf(
"==================== SLICES FOUND ====================\n");
186 for (j = 0; j < nslices; j++) {
187 printf(
" %2d: [% .15e , % .15e]\n", j + 1, sli[j], sli[j + 1]);
190 ev_int = (int)(1 + ecount / ((
double)nslices));
193 for (sl = 0; sl < nslices; sl++) {
194 printf(
"======================================================\n");
196 double *lam, *Y, *res;
201 printf(
" subinterval: [%.15e , %.15e]\n", a, b);
215 fprintf(fstats,
" polynomial deg %d, bar %.15e gam %.15e\n", pol.
deg,
221 mlan =
max(5 * nev, 100);
224 ierr =
ChebLanNr(xintv, mlan, tol, vinit, &pol, &nev2, &lam, &Y, &res,
227 printf(
"ChebLanTr error %d\n", ierr);
233 ind = (
int *)malloc(nev2 *
sizeof(
int));
235 printf(
" number of eigenvalues found: %d\n", nev2);
237 fprintf(fstats,
" Eigenvalues in [a, b]\n");
238 fprintf(fstats,
" Computed [%d] ||Res||\n", nev2);
239 for (i = 0; i < nev2; i++) {
240 fprintf(fstats,
"% .15e %.1e\n", lam[i], res[ind[i]]);
242 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - "
243 "- - - - - - - - - - - - - - - - - -\n");
244 memcpy(&alleigs[totcnt], lam, nev2 *
sizeof(
double));
258 fprintf(fstats,
" --> Total eigenvalues found = %d\n", totcnt);
259 sprintf(path,
"OUT/EigsOut_KPM_MMPLanN_(%s_%s)", io.
MatNam1, io.
MatNam2);
260 FILE *fmtout = fopen(path,
"w");
262 for (j = 0; j < totcnt; j++)
263 fprintf(fmtout,
"%.15e\n", alleigs[j]);
275 if (fstats != stdout) {
281 if (flog != stdout) {
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 get_matrix_info(FILE *fmat, io_t *pio)
int SetGenEig()
Set the problem to generalized eigenvalue problem.
int SetAMatrix(csrMat *A)
Set the matrix A.
int SetupBSolDirect(csrMat *B, void **data)
Setup the B-sol by computing the Cholesky factorization of B.
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.
int SetBMatrix(csrMat *B)
Set the B matrix.
void sort_double(int n, double *v, int *ind)
int EVSLStart()
Initialize evslData.
void FreeBSolDirectData(void *data)
Free solver data.
void free_csr(csrMat *csr)
memory deallocation for csr matrix
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)
void set_pol_def(polparams *pol)
set default values for polparams struct.
int SetLTSol(SolFuncR func, void *data)
Set the solve routine for LT.
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.
void LTSolDirect(double *b, double *x, void *data)
Solver function of L^{T} x = L^{-T}*b.
int ChebLanNr(double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats)
int SetBSol(SolFuncR func, void *data)
Set the solve routine and the associated data for B.
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...
parameters for polynomial filter
void BSolDirect(double *b, double *x, void *data)
Solver function of B.
Definitions used for direct solver interface.
sparse matrix format: the coordinate (COO) format, 0-based
double evsl_timer()
evsl timer for mac