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, max_its, ev_int, sl,
27 double a, b, lmax, lmin, ecount, tol, *sli, *mu;
39 csrMat Acsr, Bcsr, Acsr0, Bcsr0;
40 double *sqrtdiag = NULL;
44 mu = malloc((Mdeg + 1) *
sizeof(
double));
45 FILE *flog = stdout, *fmat = NULL;
49 const double tau = 1e-4;
53 printf(
"-----------------------------------------\n");
54 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");
55 printf(
"-----------------------------------------\n");
64 if (NULL == (fmat = fopen(
"matfile",
"r"))) {
65 fprintf(flog,
"Can't open matfile...\n");
70 if (NULL == fgets(line,
MAX_LINE, fmat)) {
71 fprintf(flog,
"error in reading matfile...\n");
74 if ((numat = atoi(line)) <= 0) {
75 fprintf(flog,
"Invalid count of matrices...\n");
78 for (mat = 1; mat <= numat; mat++) {
80 fprintf(flog,
"Invalid format in matfile ...\n");
84 fprintf(flog,
"MATRIX A: %s...\n", io.
MatNam1);
85 fprintf(flog,
"MATRIX B: %s...\n", io.
MatNam2);
91 if (stat(
"OUT", &st) == -1) {
96 strcpy(path,
"OUT/KPM_MMRLanN_");
98 fstats = fopen(path,
"w");
100 printf(
" failed in opening output file in OUT/\n");
103 fprintf(fstats,
"MATRIX A: %s...\n", io.
MatNam1);
104 fprintf(fstats,
"MATRIX B: %s...\n", io.
MatNam2);
105 fprintf(fstats,
"Partition the interval of interest [%f,%f] into %d slices\n",
107 counts = malloc(nslices *
sizeof(
int));
108 sli = malloc((nslices + 1) *
sizeof(
double));
113 fprintf(fstats,
"matrix read successfully\n");
116 fprintf(flog,
"read_coo error for A = %d\n", ierr);
121 fprintf(fstats,
"matrix read successfully\n");
122 if (Bcoo.
nrows != n) {
126 fprintf(flog,
"read_coo error for B = %d\n", ierr);
130 sqrtdiag = (
double *)calloc(n,
sizeof(
double));
134 }
else if (io.
Fmt ==
HB) {
135 fprintf(flog,
"HB FORMAT not supported (yet) * \n");
144 for (i=0; i<n; i++) {
145 sqrtdiag[i] = sqrt(sqrtdiag[i]);
161 vinit = (
double *)malloc(n *
sizeof(
double));
163 ierr =
LanTrbounds(50, 200, 1e-10, vinit, 1, &lmin, &lmax, fstats);
170 printf(
"The degree for LS polynomial approximations to B^{-1} and B^{-1/2} "
171 "are %d and %d\n", Bsol.
deg, Bsqrtsol.
deg);
181 ierr =
LanTrbounds(50, 200, 1e-10, vinit, 1, &lmin, &lmax, fstats);
182 fprintf(fstats,
"Step 0: Eigenvalue bound s for B^{-1}*A: [%.15e, %.15e]\n",
194 ierr =
kpmdos(Mdeg, 1, nvec, xintv, mu, &ecount);
197 printf(
"kpmdos error %d\n", ierr);
200 fprintf(fstats,
" Time to build DOS (kpmdos) was : %10.2f \n", t);
201 fprintf(fstats,
" estimated eig count in interval: %.15e \n", ecount);
204 fprintf(fstats,
"DOS parameters: Mdeg = %d, nvec = %d, npnts = %d\n", Mdeg,
206 ierr =
spslicer(sli, mu, Mdeg, xintv, nslices, npts);
208 printf(
"spslicer error %d\n", ierr);
211 printf(
"==================== SLICES FOUND ====================\n");
212 for (j = 0; j < nslices; j++) {
213 printf(
" %2d: [% .15e , % .15e]\n", j + 1, sli[j], sli[j + 1]);
220 ev_int = (int)(1 + ecount / ((
double)nslices));
222 alleigs = malloc(n *
sizeof(
double));
231 for (sl = 0; sl < nslices; sl++) {
232 printf(
"======================================================\n");
234 double *lam, *Y, *res;
239 printf(
" subinterval: [%.4e , %.4e]\n", a, b);
256 void **solshiftdata = (
void **)malloc(num *
sizeof(
void *));
264 max_its =
max(4 * nev, 300);
265 max_its =
min(max_its, n);
267 ierr =
RatLanNr(intv, max_its, tol, vinit, &rat, &nev2, &lam, &Y, &res,
270 printf(
"RatLanNr error %d\n", ierr);
276 ind = (
int *)malloc(nev2 *
sizeof(
int));
278 printf(
" number of eigenvalues found: %d\n", nev2);
280 fprintf(fstats,
" Eigenvalues in [a, b]\n");
281 fprintf(fstats,
" Computed [%d] ||Res||\n", nev2);
282 for (i = 0; i < nev2; i++) {
283 fprintf(fstats,
"% .15e %.1e\n", lam[i], res[ind[i]]);
285 fprintf(fstats,
"- - - - - - - - - - - - - - - - - - - - - - - - - - - - "
286 "- - - - - - - - - - - - - - - - - -\n");
287 memcpy(&alleigs[totcnt], lam, nev2 *
sizeof(
double));
303 fprintf(fstats,
" --> Total eigenvalues found = %d\n", totcnt);
304 sprintf(path,
"OUT/EigsOut_KPM_MMRLanN_(%s_%s)", io.
MatNam1, io.
MatNam2);
305 FILE *fmtout = fopen(path,
"w");
307 for (j = 0; j < totcnt; j++)
308 fprintf(fmtout,
"%.15e\n", alleigs[j]);
321 if (fstats != stdout) {
327 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_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)
void extrDiagCsr(csrMat *B, double *d)
int SetGenEig()
Set the problem to generalized eigenvalue problem.
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 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
int SetBMatrix(csrMat *B)
Set the B matrix.
void sort_double(int n, double *v, int *ind)
int EVSLStart()
Initialize evslData.
void free_csr(csrMat *csr)
memory deallocation for csr matrix
void diagScalCsr(csrMat *A, double *d)
void SetupPolSqrt(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
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 SetLTSol(SolFuncR func, void *data)
Set the solve routine for LT.
int SetStdEig()
Set the problem to standard eigenvalue problem.
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 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...
int get_matrix_info(FILE *fmat, io_t *pio)
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.
void csr_copy(csrMat *A, csrMat *B, int allocB)
copy a csr matrix A into B alloB: 0: will not allocate memory for B (have been alloced outside) 1: wi...
void BSolPol(double *b, double *x, void *data)
int find_ratf(double *intv, ratparams *rat)
void SetupPolRec(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data)
void FreeBSolPolData(BSolDataPol *data)
sparse matrix format: the coordinate (COO) format, 0-based
double evsl_timer()
evsl timer for mac
parameters for rational filter
int read_coo_MM(const char *matfile, int idxin, int idxout, cooMat *Acoo)