From 10456fcc71bb4162182767153a40eb25a2b1b0c1 Mon Sep 17 00:00:00 2001
From: Luke Erlandson <LukeAErlandson@gmail.com>
Date: Mon, 26 Nov 2018 21:25:14 -0500
Subject: [PATCH 1/5] Improve const correctness

---
 EVSL_1.1.1/INC/evsl.h               | 44 ++++++++++----------
 EVSL_1.1.1/INC/internal_header.h    | 64 ++++++++++++++---------------
 EVSL_1.1.1/INC/struct.h             |  6 +--
 EVSL_1.1.1/SRC/blas/evsl_blas.h     | 42 +++++++++----------
 EVSL_1.1.1/SRC/cheblanNr.c          | 27 ++++++------
 EVSL_1.1.1/SRC/cheblanTr.c          | 43 ++++++++++---------
 EVSL_1.1.1/SRC/chebpoly.c           | 58 ++++++++++++++------------
 EVSL_1.1.1/SRC/chebsi.c             | 43 ++++++++++---------
 EVSL_1.1.1/SRC/dos_utils.c          | 36 ++++++++--------
 EVSL_1.1.1/SRC/dumps.c              |  8 ++--
 EVSL_1.1.1/SRC/exDOS.c              |  2 +-
 EVSL_1.1.1/SRC/lanTrbounds.c        |  2 +-
 EVSL_1.1.1/SRC/lanbounds.c          |  2 +-
 EVSL_1.1.1/SRC/landos.c             |  4 +-
 EVSL_1.1.1/SRC/landosG.c            |  2 +-
 EVSL_1.1.1/SRC/lapack/evsl_lapack.h |  8 ++--
 EVSL_1.1.1/SRC/misc_la.c            | 60 ++++++++++++++-------------
 EVSL_1.1.1/SRC/ratfilter.c          | 22 +++++-----
 EVSL_1.1.1/SRC/ratlanNr.c           | 21 +++++-----
 EVSL_1.1.1/SRC/ratlanTr.c           | 26 ++++++------
 EVSL_1.1.1/SRC/simpson.c            |  2 +-
 EVSL_1.1.1/SRC/spmat.c              | 24 +++++------
 EVSL_1.1.1/SRC/spslice.c            |  6 +--
 EVSL_1.1.1/SRC/spslice2.c           |  2 +-
 EVSL_1.1.1/SRC/vect.c               |  4 +-
 25 files changed, 289 insertions(+), 269 deletions(-)

diff --git a/EVSL_1.1.1/INC/evsl.h b/EVSL_1.1.1/INC/evsl.h
index 35d2112..3118415 100644
--- a/EVSL_1.1.1/INC/evsl.h
+++ b/EVSL_1.1.1/INC/evsl.h
@@ -76,28 +76,28 @@
 )
 
 /*- - - - - - - - - cheblanNr.c */
-int ChebLanNr(double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
+int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit, polparams const * const pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
 
 /*- - - - - - - - - cheblanTr.c */
-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);
+int ChebLanTr(int lanm, int nev, const double * intv, int maxit, double tol, const double * vinit, polparams const * const pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
 
 /*- - - - - - - - - chebpoly.c */
 void set_pol_def(polparams *pol);
-int find_pol(double *intv, polparams *pol);
+int find_pol(const double *intv, polparams *pol);
 void free_pol(polparams *pol);
 /* int ChebAv(polparams *pol, double *v, double *y, double *w); */
 
 /*- - - - - - - - - chebsi.c */
-int ChebSI(int nev, double *intv, int maxit, double tol, double *vinit, polparams *pol, int *nevo, double **lamo, double **Yo, double **reso, FILE *fstats);
+int ChebSI(int nev, const double *intv, int maxit, double tol, const double *vinit, const polparams *pol, int *nevo, double **lamo, double **Yo, double **reso, FILE *fstats);
 
 /*- - - - - - - - - - dos_utils.c */
 void SetupBPol(int n, int max_deg, double tol, double lmin, double lmax, double (*ffun)(double), BSolDataPol *data);
 void SetupPolRec(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data);
 void SetupPolSqrt(int n, int max_deg, double tol, double lmin, double lmax, BSolDataPol *data);
 void FreeBSolPolData(BSolDataPol *data);
-void BSolPol(double *b, double *x, void *data);
-void extrDiagCsr(csrMat *B, double *d);
-void diagScalCsr(csrMat *A, double *d);
+void BSolPol(const double *b, double *x, void *data);
+void extrDiagCsr(const csrMat *B, double *d);
+void diagScalCsr(csrMat *A, const double *d);
 
 /*- - - - - - - - - evsl_memory.c */
 void *_evsl_Malloc(size_t nbytes);
@@ -106,10 +106,10 @@ void *_evsl_Realloc(void *ptr, size_t nbytes);
 void _evsl_Free(void *ptr);
 
 /*- - - - - - - - - exDos.c */
-int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv);
+int exDOS(const double *vals, int n, int npts, double *x, double *y, const double *intv);
 
 /*- - - - - - - - - lanbounds.c */
-int LanBounds(int msteps, double *v, double *lmin, double *lmax);
+int LanBounds(int msteps, const double *v, double *lmin, double *lmax);
 
 /*- - - - - - - - - - landos.c */
 // Computes the density of states (DOS, or spectral density)
@@ -119,41 +119,41 @@ int LanDos(const int nvec, int msteps, const int npts, double *xdos, double *ydo
 int LanDosG(const int nvec, int msteps, const int npts, double *xdos, double *ydos, double *neig, const double *const intv);
 
 /*- - - - - - - - - lanTrbounds.c */
-int LanTrbounds(int lanm, int maxit, double tol, double *vinit, int bndtype, double *lammin, double *lammax, FILE *fstats);
+int LanTrbounds(int lanm, int maxit, double tol, const double *vinit, int bndtype, double *lammin, double *lammax, FILE *fstats);
 
 /*- - - -- - - - - - misc_la.c */
 int scalEigVec(int n, int nev, double *Y, double* sqrtdiag);
 
 /*- - - - - - - - - ratfilter.c */
 void set_ratf_def(ratparams *rat);
-int find_ratf(double *intv, ratparams *rat);
+int find_ratf(const double *intv, ratparams *rat);
 int set_ratf_solfunc(ratparams *rat, csrMat *A, csrMat *B, SolFuncC *funcs, void **data);
 void free_rat(ratparams *rat);
 
 /*- - - - - - - - - ratlanNr.c */
-int RatLanNr(double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
+int RatLanNr(const double *intv, int maxit, double tol, const double *vinit, ratparams *rat, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
 
 /*- - - - - - - - - ratlanTr.c */
-int RatLanTr(int lanm, int nev, double *intv, int maxit, double tol, double *vinit, ratparams *rat, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
+int RatLanTr(int lanm, int nev, const double *intv, int maxit, double tol, const double *vinit, ratparams *rat, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
 
 /*- - - - - - - - - spmat.c */
-void csr_copy(csrMat *A, csrMat *B, int allocB);
+void csr_copy(const csrMat *A, csrMat *B, int allocB);
 // convert a COO matrix to a CSR matrix
-int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr);
+int cooMat_to_csrMat(int cooidx, const cooMat *coo, csrMat *csr);
 // free a CSR
 void free_csr(csrMat *csr);
 // free a COO
 void free_coo(cooMat *coo);
 // matvec y = A*x
-void matvec_csr(double *x, double *y, void *data);
+void matvec_csr(const double *x, double *y, void *data);
 /* add two csr matrices */
-int matadd(double alp, double bet, csrMat *A, csrMat *B, csrMat *C, int *mapA, int *mapB);
+int matadd(double alp, double bet, const csrMat *A, const csrMat *B, csrMat *C, int *mapA, int *mapB);
 /* sparse identity */
 int speye(int n, csrMat *A);
 /* extract upper triangular part of A */
-void triuCsr(csrMat *A, csrMat *U);
+void triuCsr(const csrMat *A, csrMat *U);
 
-int arrays_copyto_csrMat(int nrow, int ncol, int *ia, int *ja, double *a, csrMat *A);
+int arrays_copyto_csrMat(int nrow, int ncol, const int *ia, const int *ja, const double *a, csrMat *A);
 /*- - - - - - - - - evsl.c */
 /* set an external matvec function */
 int SetAMatvec(int n, MVFunc func, void *data);
@@ -175,11 +175,11 @@ int EVSLFinish();
 void SetDiagScal(double *ds);
 
 /*- - - - - - - - - spslicer.c */
-int spslicer(double *sli, double *mu, int Mdeg, double *intv, int n_int, int npts);
-int kpmdos(int Mdeg, int damping, int nvec, double *ab, double *mu, double *ecnt);
+int spslicer(double *sli, const double *mu, int Mdeg, const double *intv, int n_int, int npts);
+int kpmdos(int Mdeg, int damping, int nvec, const double *ab, double *mu, double *ecnt);
 
 /*- - - - - - - - - - spslicer2.c */
-void spslicer2(double *xi, double *yi, int n_int, int npts, double *sli);
+void spslicer2(const double *xi, double *yi, int n_int, int npts, double *sli);
 
 /*- - - - - - - - - timing.c */
 double evsl_timer();
diff --git a/EVSL_1.1.1/INC/internal_header.h b/EVSL_1.1.1/INC/internal_header.h
index ff0869d..663a334 100644
--- a/EVSL_1.1.1/INC/internal_header.h
+++ b/EVSL_1.1.1/INC/internal_header.h
@@ -29,43 +29,43 @@
 /*- - - - - - - - - chebpoly.c */
 int chebxCoefd(int m, double gam, int damping, double *mu);
 int dampcf(int m, int damping, double *jac);
-int chebxPltd(int m, double *mu, int n, double *xi, double *yi);
-int ChebAv(polparams *pol, double *v, double *y, double *w);
+int chebxPltd(int m, const double *mu, int n, const double *xi, double *yi);
+int ChebAv(const polparams *pol, const double *v, double *y, double *w);
 void chext(polparams *pol, double aIn, double bIn);
 
 /*- - - - - - - - - dos_utils.c */
 int apfun(const double c, const double h, const double* xi, double (*ffun)(double), const int npts, double* yi);
 double rec(const double a);
 double isqrt(const double a);
-int pnav(double *mu, const int m, const double cc, const double dd, double *v, double *y, double *w);
+int pnav(const double *mu, const int m, const double cc, const double dd, const double *v, double *y, double *w);
 int lsPol(double (*ffun)(double), BSolDataPol *pol);
 
 /*- - - - - - - - - dump.c */
-void savemat(csrMat *A, const char *fn);
-void savedensemat(double *A, int lda, int m, int n, const char *fn);
+void savemat(const csrMat *A, const char *fn);
+void savedensemat(const double *A, int lda, int m, int n, const char *fn);
 void save_vec(int n, const double *x, const char fn[]);
 
 /*- - - - - - - - - misc_la.c */
 int SymmTridEig(double *eigVal, double *eigVec, int n, const double *diag, const double *sdiag);
 int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu, int *nevO, const double *diag, const double *sdiag);
-void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam);
-void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, double *w);
-void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q, double *v, double *w);
-void orth(double *V, int n, int k, double *Vo, double *work);
+void SymEigenSolver(int n, const double *A, int lda, double *Q, int ldq, double *lam);
+void CGS_DGKS(int n, int k, int i_max, const double *Q, double *v, double *nrmv, double *w);
+void CGS_DGKS2(int n, int k, int i_max, const double *Z, const double *Q, double *v, double *w);
+void orth(const double *V, int n, int k, double *Vo, double *work);
 
 /*- - - - - - - - - ratfilter.c */
 void contQuad(int method, int n, EVSL_Complex *zk);
-void ratf2p2(int n, int *mulp, EVSL_Complex *zk, EVSL_Complex *alp, int m, double *z, double *x);
+void ratf2p2(int n, const int *mulp, const EVSL_Complex *zk, const EVSL_Complex *alp, int m, const double *z, double *x);
 void pfe2(EVSL_Complex s1, EVSL_Complex s2, int k1, int k2, EVSL_Complex* alp, EVSL_Complex* bet);
 EVSL_Complex integg2(EVSL_Complex s1, EVSL_Complex s2, EVSL_Complex* alp, int k1, EVSL_Complex* bet, int k2, double a, double b);
-void weights(int n, EVSL_Complex* zk, int* pow, double lambda, EVSL_Complex* omega);
-int scaleweigthts(int n, double a, double b, EVSL_Complex *zk, int* pow, EVSL_Complex* omegaM);
+void weights(int n, const EVSL_Complex* zk, const int* pow, double lambda, EVSL_Complex* omega);
+int scaleweigthts(int n, double a, double b, EVSL_Complex *zk, const int* pow, EVSL_Complex* omegaM);
 
 /*- - - - - - - - - ratlanNr.c */
-void RatFiltApply(int n, ratparams *rat, double *b, double *x, double *w3);
+void RatFiltApply(int n, const ratparams *rat, const double *b, double *x, double *w3);
 
 /*- - - - - - - - - - simpson.c */
-void simpson(double *xi, double *yi, int npts);
+void simpson(const double *xi, double *yi, int npts);
 
 /*- - - - - - - - - spmat.c */
 // memory allocation/reallocation for a CSR matrix
@@ -77,19 +77,19 @@ int time_seeder();
 
 /*- - - - - - - - - vect.c */
 void vecset(int n, double t, double *v);
-void vec_perm(int n, int *p, double *x, double *y);
-void vec_iperm(int n, int *p, double *x, double *y);
+void vec_perm(int n, const int *p, const double *x, double *y);
+void vec_iperm(int n, const int *p, const double *x, double *y);
 
 /*------------------- inline functions */
 /**
 * @brief y = A * x
 * This is the matvec function for the matrix A in evsldata
 */
-static inline void matvec_A(double *x, double *y) {
+static inline void matvec_A(const double *x, double *y) {
   CHKERR(!evsldata.Amv);
-  double tms = evsl_timer();
+  const double tms = evsl_timer();
   evsldata.Amv->func(x, y, evsldata.Amv->data);
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_mvA += tme - tms;
   evslstat.n_mvA ++;
 }
@@ -98,11 +98,11 @@ static inline void matvec_A(double *x, double *y) {
 * @brief y = B * x
 * This is the matvec function for the matrix B in evsldata
 */
-static inline void matvec_B(double *x, double *y) {
+static inline void matvec_B(double const * const x, double *y) {
   CHKERR(!evsldata.Bmv);
-  double tms = evsl_timer();
+  const double tms = evsl_timer();
   evsldata.Bmv->func(x, y, evsldata.Bmv->data);
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_mvB += tme - tms;
   evslstat.n_mvB ++;
 }
@@ -111,11 +111,11 @@ static inline void matvec_B(double *x, double *y) {
 * @brief y = B \ x
 * This is the solve function for the matrix B in evsldata
 */
-static inline void solve_B(double *x, double *y) {
+static inline void solve_B(double const * const x, double *y) {
   CHKERR(!evsldata.Bsol);
-  double tms = evsl_timer();
+  const double tms = evsl_timer();
   evsldata.Bsol->func(x, y, evsldata.Bsol->data);
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_svB += tme - tms;
   evslstat.n_svB ++;
 }
@@ -124,11 +124,11 @@ static inline void solve_B(double *x, double *y) {
 * @brief y = LT \ x or y = SQRT(B) \ x
 * This is the solve function for the matrix B in evsldata
 */
-static inline void solve_LT(double *x, double *y) {
+static inline void solve_LT(double const * const x, double *y) {
   CHKERR(!evsldata.LTsol);
-  double tms = evsl_timer();
+  const double tms = evsl_timer();
   evsldata.LTsol->func(x, y, evsldata.LTsol->data);
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_svLT += tme - tms;
   evslstat.n_svLT ++;
 }
@@ -138,11 +138,11 @@ static inline void solve_LT(double *x, double *y) {
 * This is the solve function for the complex shifted matrix
 */
 static inline void solve_ASigB(EVSLASIGMABSol *sol, int n,
-                               double *br, double *bz,
+                               double const * const br, double const * const bz,
                                double *xr, double *xz) {
-  double tms = evsl_timer();
+  const double tms = evsl_timer();
   (sol->func)(n, br, bz, xr, xz, sol->data);
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_svASigB+= tme - tms;
   evslstat.n_svASigB ++;
 }
@@ -150,7 +150,7 @@ static inline void solve_ASigB(EVSLASIGMABSol *sol, int n,
 /**
  * @brief check if an interval is valid
  * */
-static inline int check_intv(double *intv, FILE *fstats) {
+static inline int check_intv(double const * const intv, FILE *fstats) {
   /* intv[4]: ( intv[0], intv[1] ) is the inteval of interest
    *          ( intv[2], intv[3] ) is the spectrum bounds
    * return   0: ok
diff --git a/EVSL_1.1.1/INC/struct.h b/EVSL_1.1.1/INC/struct.h
index 2a76b83..05a1165 100644
--- a/EVSL_1.1.1/INC/struct.h
+++ b/EVSL_1.1.1/INC/struct.h
@@ -80,17 +80,17 @@ typedef struct _polparams {
  * be the  solution (complex vector),  and "data" contains  all the
  * data  needed  by  the  solver.
  */
-typedef void (*SolFuncC)(int n, double *br, double *bz, double *xr, double *xz, void *data);
+typedef void (*SolFuncC)(int n, double const * const br, double const * const bz, double *xr, double *xz, void *data);
 
 /**
  * @brief function prototype for applying the solve B x = b
  */
-typedef void (*SolFuncR)(double *b, double *x, void *data);
+typedef void (*SolFuncR)(double const * const b, double *x, void *data);
 
 /**
  * @brief matvec function prototype
  */
-typedef void (*MVFunc)(double *x, double *y, void *data);
+typedef void (*MVFunc)(const double *x, double *y, void *data);
 
 /*!
  * @brief user-provided function and data for solving (A - SIGMA*B) x = b
diff --git a/EVSL_1.1.1/SRC/blas/evsl_blas.h b/EVSL_1.1.1/SRC/blas/evsl_blas.h
index e8a3a03..3675717 100644
--- a/EVSL_1.1.1/SRC/blas/evsl_blas.h
+++ b/EVSL_1.1.1/SRC/blas/evsl_blas.h
@@ -33,35 +33,35 @@
 #if defined(__cplusplus) && !defined(EVSL_USING_EVSL_BLAS)
 extern "C" {
 #endif
-int evsl_daxpy(EVSL_Int *n, EVSL_Real *da, EVSL_Real *dx, EVSL_Int *incx, EVSL_Real *dy, EVSL_Int *incy);
-int evsl_dcopy(EVSL_Int *n, EVSL_Real *dx, EVSL_Int *incx, EVSL_Real *dy, EVSL_Int *incy);
-EVSL_Real evsl_ddot(EVSL_Int *n, EVSL_Real *dx, EVSL_Int *incx, EVSL_Real *dy, EVSL_Int *incy);
-int evsl_dgemm(const char *transa, const char *transb, EVSL_Int *m, EVSL_Int * n, EVSL_Int *k, EVSL_Real *alpha, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *b, EVSL_Int *ldb, EVSL_Real *beta, EVSL_Real *c__, EVSL_Int *ldc);
-int evsl_dgemv(const char *trans, EVSL_Int *m, EVSL_Int *n, EVSL_Real *alpha, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *x, EVSL_Int *incx, EVSL_Real *beta, EVSL_Real *y, EVSL_Int *incy);
-EVSL_Real evsl_dnrm2(EVSL_Int *n, EVSL_Real *x, EVSL_Int *incx);
-int evsl_dscal(EVSL_Int *n, EVSL_Real *da, EVSL_Real *dx, EVSL_Int *incx);
-int evsl_dcopy(EVSL_Int *n, EVSL_Real *dx, EVSL_Int *incx, EVSL_Real *dy, EVSL_Int *incy);
+int evsl_daxpy(const EVSL_Int *n, const EVSL_Real *da, const EVSL_Real *dx, const EVSL_Int *incx, EVSL_Real *dy, const EVSL_Int *incy);
+int evsl_dcopy(const EVSL_Int *n, const EVSL_Real *dx, const EVSL_Int *incx, EVSL_Real *dy, const EVSL_Int *incy);
+EVSL_Real evsl_ddot(const EVSL_Int *n, const EVSL_Real *dx, const EVSL_Int *incx, const EVSL_Real *dy, const EVSL_Int *incy);
+int evsl_dgemm(const char *transa, const char *transb, const EVSL_Int *m, const EVSL_Int * n, const EVSL_Int *k, const EVSL_Real *alpha, const EVSL_Real *a, const EVSL_Int *lda, const EVSL_Real *b, const EVSL_Int *ldb, const EVSL_Real *beta, EVSL_Real *c__, const EVSL_Int *ldc);
+int evsl_dgemv(const char *trans, const EVSL_Int *m, const EVSL_Int *n, const EVSL_Real *alpha, const EVSL_Real *a, const EVSL_Int *lda, const EVSL_Real *x, const EVSL_Int *incx, const EVSL_Real *beta, EVSL_Real *y, const EVSL_Int *incy);
+EVSL_Real evsl_dnrm2(const EVSL_Int *n, const EVSL_Real *x, const EVSL_Int *incx);
+int evsl_dscal(const EVSL_Int *n, const EVSL_Real *da, EVSL_Real *dx, const EVSL_Int *incx);
+int evsl_dcopy(const EVSL_Int *n, const EVSL_Real *dx, const EVSL_Int *incx, EVSL_Real *dy, const EVSL_Int *incy);
 #if defined(LAPACK_F2C_INCLUDE)
 /* blas needed from evsl lapack only */
-int evsl_dtrmm(const char *side, const char *uplo, const char *transa, const char *diag, EVSL_Int *m, EVSL_Int *n, EVSL_Real *alpha, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *b, EVSL_Int *ldb);
-int evsl_dsyr2k(const char *uplo, const char *trans, EVSL_Int *n, EVSL_Int *k, EVSL_Real *alpha, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *b, EVSL_Int *ldb, EVSL_Real *beta, EVSL_Real *c__, EVSL_Int *ldc);
-int evsl_dsymv(const char *uplo, EVSL_Int *n, EVSL_Real *alpha, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *x, EVSL_Int *incx, EVSL_Real *beta, EVSL_Real *y, EVSL_Int *incy);
-int evsl_dger(EVSL_Int *m, EVSL_Int *n, EVSL_Real *alpha, EVSL_Real *x, EVSL_Int *incx, EVSL_Real *y, EVSL_Int *incy, EVSL_Real *a, EVSL_Int *lda);
+int evsl_dtrmm(const char *side, const char *uplo, const char *transa, const char *diag, const EVSL_Int *m, const EVSL_Int *n, const EVSL_Real *alpha, const EVSL_Real *a, const EVSL_Int *lda, EVSL_Real *b, const EVSL_Int *ldb);
+int evsl_dsyr2k(const char *uplo, const char *trans, const EVSL_Int *n, const EVSL_Int *k, const EVSL_Real *alpha, const EVSL_Real *a, const EVSL_Int *lda, const EVSL_Real *b, const EVSL_Int *ldb, const EVSL_Real *beta, EVSL_Real *c__, const EVSL_Int *ldc);
+int evsl_dsymv(const char *uplo, EVSL_Int *n, const EVSL_Real *alpha, const EVSL_Real *a, const EVSL_Int *lda, const EVSL_Real *x, const EVSL_Int *incx, const EVSL_Real *beta, EVSL_Real *y, const EVSL_Int *incy);
+int evsl_dger(const EVSL_Int *m, const EVSL_Int *n, const EVSL_Real *alpha, const EVSL_Real *x, const EVSL_Int *incx, const EVSL_Real *y, const EVSL_Int *incy, EVSL_Real *a, const EVSL_Int *lda);
 EVSL_Int lsame_(const char *ca, const char *cb);
 int evsl_xerbla(const char *srname, EVSL_Int *info);
-int evsl_zgemm(const char *transa, const char *transb, EVSL_Int *m, EVSL_Int *n, EVSL_Int *k, doublecomplex *alpha, doublecomplex *a, EVSL_Int *lda, doublecomplex *b, EVSL_Int *ldb, doublecomplex *beta, doublecomplex *c__, EVSL_Int *ldc);
-int evsl_zgeru(EVSL_Int *m, EVSL_Int *n, doublecomplex *alpha, doublecomplex *x, EVSL_Int *incx, doublecomplex *y, EVSL_Int *incy, doublecomplex *a, EVSL_Int *lda);
-int evsl_ztrsm(const char *side, const char *uplo, const char *transa, const char *diag, EVSL_Int *m, EVSL_Int *n, doublecomplex *alpha, doublecomplex *a, EVSL_Int *lda, doublecomplex *b, EVSL_Int *ldb);
-EVSL_Int evsl_izamax(EVSL_Int *n, doublecomplex *zx, EVSL_Int *incx);
-int evsl_zswap(EVSL_Int *n, doublecomplex *zx, EVSL_Int *incx, doublecomplex *zy, EVSL_Int *incy);
-int evsl_zscal(EVSL_Int *n, doublecomplex *za, doublecomplex *zx, EVSL_Int *incx);
-int evsl_dtrmv(const char *uplo, const char *trans, const char *diag, EVSL_Int *n, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *x, EVSL_Int *incx);
+int evsl_zgemm(const char *transa, const char *transb, const EVSL_Int *m, const EVSL_Int *n, const EVSL_Int *k, const doublecomplex *alpha, const doublecomplex *a, const EVSL_Int *lda, const doublecomplex *b, const EVSL_Int *ldb, const doublecomplex *beta, doublecomplex *c__, const EVSL_Int *ldc);
+int evsl_zgeru(const EVSL_Int *m, const EVSL_Int *n, const doublecomplex *alpha, const doublecomplex *x, const EVSL_Int *incx, const doublecomplex *y, const EVSL_Int *incy, doublecomplex *a, const EVSL_Int *lda);
+int evsl_ztrsm(const char *side, const char *uplo, const char *transa, const char *diag, const EVSL_Int *m, const EVSL_Int *n, const doublecomplex *alpha, const doublecomplex *a, const EVSL_Int *lda, doublecomplex *b, const EVSL_Int *ldb);
+EVSL_Int evsl_izamax(const EVSL_Int *n, const doublecomplex *zx, const EVSL_Int *incx);
+int evsl_zswap(const EVSL_Int *n, doublecomplex *zx, const EVSL_Int *incx, doublecomplex *zy, const EVSL_Int *incy);
+int evsl_zscal(const EVSL_Int *n, const doublecomplex *za, doublecomplex *zx, const EVSL_Int *incx);
+int evsl_dtrmv(const char *uplo, const char *trans, const char *diag, const EVSL_Int *n, const EVSL_Real *a, const EVSL_Int *lda, EVSL_Real *x, const EVSL_Int *incx);
+// dswap is not specified with correct const due to conflicting decl in dsteqr.c
 int evsl_dswap(EVSL_Int *n, EVSL_Real *dx, EVSL_Int *incx, EVSL_Real *dy, EVSL_Int *incy);
-int evsl_dsyr2(const char *uplo, EVSL_Int *n, EVSL_Real *alpha, EVSL_Real *x, EVSL_Int *incx, EVSL_Real *y, EVSL_Int *incy, EVSL_Real *a, EVSL_Int *lda);
+int evsl_dsyr2(const char *uplo, const EVSL_Int *n, const EVSL_Real *alpha, const EVSL_Real *x, const EVSL_Int *incx, const EVSL_Real *y, const EVSL_Int *incy, EVSL_Real *a, const EVSL_Int *lda);
 #endif
 #if defined(__cplusplus) && !defined(EVSL_USING_EVSL_BLAS)
 }
 #endif
 
 #endif
-
diff --git a/EVSL_1.1.1/SRC/cheblanNr.c b/EVSL_1.1.1/SRC/cheblanNr.c
index 2044590..8d773b4 100644
--- a/EVSL_1.1.1/SRC/cheblanNr.c
+++ b/EVSL_1.1.1/SRC/cheblanNr.c
@@ -54,8 +54,8 @@
  * @warning memory allocation for Wo/lamo/reso within this function
  *
  * ------------------------------------------------------------ */
-int ChebLanNr(double *intv, int maxit, double tol, double *vinit,
-              polparams *pol, int *nevOut, double **lamo, double **Wo,
+int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit,
+              polparams const * const pol, int *nevOut, double **lamo, double **Wo,
               double **reso, FILE *fstats) {
   //-------------------- to report timings/
   double tall, tm1 = 0.0, tt;
@@ -70,15 +70,16 @@ int ChebLanNr(double *intv, int maxit, double tol, double *vinit,
     do_print = 0;
   }
   /*-------------------- frequently used constants  */
-  char cN = 'N';
-  int one = 1;
-  double done = 1.0, dzero = 0.0;
+  const char cN = 'N';
+  const int one = 1;
+  const double done = 1.0;
+  const double dzero = 0.0;
   /*-------------------- Ntest = when to start testing convergence */
-  int Ntest = 30;
+  const int Ntest = 30;
   /*-------------------- how often to test */
-  int cycle = 20;
+  const int cycle = 20;
   /* size of the matrix */
-  int n = evsldata.n;
+  const int n = evsldata.n;
   /* max num of its */
   maxit = evsl_min(n, maxit);
   /*-------------------- Caveat !!!: To prevent integer overflow, we save
@@ -91,17 +92,17 @@ int ChebLanNr(double *intv, int maxit, double tol, double *vinit,
   /*-------------------- polynomial filter  approximates the delta
                          function centered at 'gamB'.
                          bar: a bar value to threshold Ritz values of p(A) */
-  double bar = pol->bar;
-  double gamB = pol->gam;
+  const double bar = pol->bar;
+  const double gamB = pol->gam;
   /*-------------------- interval [aa, bb], used for testing only at end */
   if (check_intv(intv, fstats) < 0) {
     *nevOut = 0;
     *lamo = NULL; *Wo = NULL; *reso = NULL;
     return 0;
   }
-  double aa = intv[0];
-  double bb = intv[1];
-  int deg = pol->deg;
+  const double aa = intv[0];
+  const double bb = intv[1];
+  const int deg = pol->deg;
   if (do_print) {
     fprintf(fstats, " ** Cheb Poly of deg = %d, gam = %.15e, bar: %.15e\n",
             deg, gamB, bar);
diff --git a/EVSL_1.1.1/SRC/cheblanTr.c b/EVSL_1.1.1/SRC/cheblanTr.c
index d1b628b..8bba4c5 100644
--- a/EVSL_1.1.1/SRC/cheblanTr.c
+++ b/EVSL_1.1.1/SRC/cheblanTr.c
@@ -56,8 +56,8 @@
  * @warning memory allocation for W/vals/resW within this function
  *
  **/
-int ChebLanTr(int lanm, int nev, double *intv, int maxit,
-              double tol, double *vinit, polparams *pol, int *nev2,
+int ChebLanTr(int lanm, int nev, const double * intv, int maxit,
+              double tol, const double * vinit, polparams const * const pol, int *nev2,
               double **vals, double **W, double **resW, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
   /*-------------------- for stats */
@@ -71,31 +71,33 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
     do_print = 0;
   }
   /* size of the matrix */
-  int n = evsldata.n;
-  size_t n_l = n;
+  const int n = evsldata.n;
+  const size_t n_l = n;
   /*--------------------- adjust lanm and maxit */
   lanm = evsl_min(lanm, n);
-  int lanm1=lanm+1;
-  size_t lanm1_l = lanm1;
+  const int lanm1=lanm+1;
+  const size_t lanm1_l = lanm1;
   /*  if use full lanczos, should not do more than n iterations */
   if (lanm == n) {
     maxit = evsl_min(maxit, n);
   }
   /*-------------------- this is needed to increment the space when we
                          discover more than nev eigenvalues in interval */
-  double nevInc = 0.2;   /* add 1  + 20% each time it is needed */
+  const double nevInc = 0.2;   /* add 1  + 20% each time it is needed */
   /*-------------------- if we have at least nev/ev_frac good candidate
                          eigenvalues from p(A) == then we restart to lock them in */
   //int evFrac = 2;
   /*--------------------   some constants frequently used */
   /* char cT='T'; */
-  char cN = 'N';
-  int one = 1;
-  double done=1.0, dmone=-1.0, dzero=0.0;
+  const char cN = 'N';
+  const int one = 1;
+  const double done=1.0;
+  const double dmone=-1.0;
+  const double dzero=0.0;
   /*-------------------- Ntest = when to start testing convergence */
   int Ntest = evsl_min(lanm, nev+50);
   /*--------------------   how often to test */
-  int cycle = 50;
+  const int cycle = 50;
   int i, ll, /* count, last_count,*/ jl, last_jl;
   /*-----------------------------------------------------------------------
     -----------------------------------------------------------------------*/
@@ -105,13 +107,14 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
     *vals = NULL; *W = NULL; *resW = NULL;
     return 0;
   }
-  double aa = intv[0];
-  double bb = intv[1];
+  const double aa = intv[0];
+  const double bb = intv[1];
   /*-------------------- Assumption: polynomial pol computed before calling cheblanTr
                          pol.  approximates the delta function centered at 'gamB'
                          bar: a bar value to threshold Ritz values of p(A) */
   //int deg = pol->deg;
-  double gamB=pol->gam, bar=pol->bar;
+  const double gamB=pol->gam;
+  const double bar=pol->bar;
   /*-------------------- gamB must be within [-1, 1] */
   if (gamB > 1.0 || gamB < -1.0) {
     if (fstats) {
@@ -347,7 +350,7 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
       }
       /*-------------------- znew = znew - beta*zold */
       if (zold) {
-        double nbeta = -beta;
+        const double nbeta = -beta;
         evsl_daxpy(&n, &nbeta, zold, &one, znew, &one);
       }
       /*-------------------- alpha = znew'*v */
@@ -356,7 +359,7 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
       T[(k-1)*lanm1_l+(k-1)] = alpha;
       wn += fabs(alpha);
       /*-------------------- znew = znew - alpha*z */
-      double nalpha = -alpha;
+      const double nalpha = -alpha;
       evsl_daxpy(&n, &nalpha, z, &one, znew, &one);
       /*-------------------- FULL reortho to all previous Lan vectors */
       if (ifGenEv) {
@@ -395,7 +398,7 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
           CGS_DGKS2(n, k, NGS_MAX, V, Z, vnew, work);
           matvec_B(vnew, znew);
           beta = sqrt(evsl_ddot(&n, vnew, &one, znew, &one));
-          double ibeta = 1.0 / beta;
+          const double ibeta = 1.0 / beta;
           evsl_dscal(&n, &ibeta, vnew, &one);
           evsl_dscal(&n, &ibeta, znew, &one);
           beta = 0.0;
@@ -405,13 +408,13 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
           /*   vnew = vnew - V(:,1:k)*V(:,1:k)'*vnew */
           /*   beta = norm(w) */
           CGS_DGKS(n, k, NGS_MAX, V, vnew, &beta, work);
-          double ibeta = 1.0 / beta;
+          const double ibeta = 1.0 / beta;
           evsl_dscal(&n, &ibeta, vnew, &one);
           beta = 0.0;
         }
       } else {
         /*---------------------- vnew = vnew / beta */
-        double ibeta = 1.0 / beta;
+        const double ibeta = 1.0 / beta;
         evsl_dscal(&n, &ibeta, vnew, &one);
         if (ifGenEv) {
           /*-------------------- znew = znew / beta */
@@ -422,7 +425,7 @@ int ChebLanTr(int lanm, int nev, double *intv, int maxit,
       T[k*lanm1_l+(k-1)] = beta;
       T[(k-1)*lanm1_l+k] = beta;
       /*---------------------- Restarting test */
-      int k1 = k-trlen-Ntest;
+      const int k1 = k-trlen-Ntest;
       if ( ((k1>=0) && (k1 % cycle == 0)) || (k == lanm) || it == maxit) {
         /*-------------------- solve eigen-problem for T(1:k,1:k)
                                vals in Rval, vecs in EvecT */
diff --git a/EVSL_1.1.1/SRC/chebpoly.c b/EVSL_1.1.1/SRC/chebpoly.c
index b1c3320..47e9740 100644
--- a/EVSL_1.1.1/SRC/chebpoly.c
+++ b/EVSL_1.1.1/SRC/chebpoly.c
@@ -75,7 +75,7 @@ int dampcf(int m, int damping, double *jac) {
  * @param[in] jac vector of damping coefficients
 **/
 
-double dif_eval(double *v, int m, double thc, double *jac){
+double dif_eval(const double * v, int m, double thc, const double * jac){
   double fval = 0.0;
   int j;
   for (j=0; j<=m; j++)
@@ -102,8 +102,10 @@ double dif_eval(double *v, int m, double thc, double *jac){
  *  @param[out] yi       = pn(xi(:) )
  *  @return 0
  **/
-int chebxPltd(int m, double *mu, int npts, double *xi, double *yi) {
-  int k, j, one = 1, n = npts;
+int chebxPltd(int m, const double *mu, int npts, const double *xi, double *yi) {
+  int k, j;
+  const int one = 1;
+  const int n = npts;
   double scal, *vkm1, *vkp1, *vk;
   vkm1 = evsl_Malloc(n, double);
   vkp1 = evsl_Malloc(n, double);
@@ -171,19 +173,20 @@ void chext(polparams *pol, double aIn, double bIn){
   int max_deg = pol->max_deg;
   // int min_deg = pol->min_deg;   NOT used
   double thresh = pol->thresh_ext;
-  double *mu = pol->mu;
+  double * mu = pol->mu;
   //-------------------- local variables
-  double del = 0.1*sqrt((bIn-aIn)*0.5);
+  const double del = 0.1*sqrt((bIn-aIn)*0.5);
   //double del = 0.0;
-  double eps = 1e-13;
+  const double eps = 1e-13;
   //double eps = 0.0;
   double a, b, x, e, c, sigma, sigma1, sigma_new, g0, g1, gnew, s1, s2, s3;
   double *t0, *t1, *tnew; // coef. of three consecutive cheby. expansions
   double bar, gam;
   int mbest=0;
-  int m1 = max_deg+1, j, k;
+  const int m1 = max_deg+1;
+  int j, k;
   //-------------------- local work space
-  int work_size = 3*m1;
+  const int work_size = 3*m1;
   //-------------------- this is for the forced degree case
   if (pol->deg > 0){
     thresh = -1.0;
@@ -264,7 +267,7 @@ void chext(polparams *pol, double aIn, double bIn){
  * @param[in] size size of the arry
  * @return index of smallest entry in array
 **/
-int indexofSmallestElement(double *array, int size){
+int indexofSmallestElement(const double * array, int size){
   int index = 0, i;
   for(i = 1; i < size; i++){
     if(array[i] < array[index])
@@ -285,7 +288,7 @@ int indexofSmallestElement(double *array, int size){
  *
  * @warning mu must be preallocated
 **/
-int rootchb(int m, double *v, double* jac, double tha, double thb, double *mu,
+int rootchb(int m, const double * v, const double * jac, double tha, double thb, double *mu,
 	    double *thcOut){
   int MaxIterBalan = 30;     // max steps in Newton to balance interval
   double tolBal;
@@ -360,9 +363,11 @@ int rootchb(int m, double *v, double* jac, double tha, double thb, double *mu,
  *  ChebAv]
  *
 **/
-int find_pol(double *intv, polparams *pol) {
+int find_pol(const double * intv, polparams *pol) {
   double *mu, *v, *jac, t=0.0, itv[2],  vals[2];
-  int max_deg=pol->max_deg, min_deg=pol->min_deg, damping=pol->damping;
+  int max_deg=pol->max_deg;
+  int min_deg=pol->min_deg;
+  const int damping=pol->damping;
   double tha=0.0, thb=0.0, thc=0.0;
   double gam,  thresh;
   int m, j, nitv, mbest;
@@ -378,17 +383,18 @@ int find_pol(double *intv, polparams *pol) {
   /*-------------------- A parameter for interval check */
   // double IntTol = 2*DBL_EPSILON; // If b*IntTol>1 accept [a b] extreme
   //double IntTol = 0.0005;
-  double IntTol = pol->intvtol;
+  const double IntTol = pol->intvtol;
   double aa, bb;
   aa = evsl_max(intv[0], intv[2]);  bb = evsl_min(intv[1], intv[3]);
   if (intv[0] < intv[2] || intv[1] > intv[3]) {
     fprintf(stdout, " warning [%s (%d)]: interval (%e, %e) is adjusted to (%e, %e)\n",
         __FILE__, __LINE__, intv[0], intv[1], aa, bb);
   }
-  double lmin = intv[2], lmax = intv[3];
+  const double lmin = intv[2];
+  const double lmax = intv[3];
   /*-------------------- cc, rr: center and half-width of [lmin, lmax] */
-  double cc = 0.5 * (lmax + lmin);
-  double dd = 0.5 * (lmax - lmin);
+  const double cc = 0.5 * (lmax + lmin);
+  const double dd = 0.5 * (lmax - lmin);
   pol->cc = cc;
   pol->dd = dd;
   /*-------------------- adjust intervals just in case. */
@@ -514,18 +520,18 @@ void free_pol(polparams *pol) {
  * @param[in, out] w Work vector of length 3*n [allocate before call]
  * @param[in, out] v is untouched
  **/
-int ChebAv(polparams *pol, double *v, double *y, double *w) {
-  double tt = evsl_timer();
+int ChebAv(const polparams *pol, const double *v, double *y, double *w) {
+  const double tt = evsl_timer();
   const int ifGenEv = evsldata.ifGenEv;
-  int n = evsldata.n;
+  const int n = evsldata.n;
   /*-------------------- unpack pol */
-  double *mu = pol->mu;
-  double dd = pol->dd;
-  double cc = pol->cc;
-  double ncc = -cc;
-  int m = pol->deg;
-  int one = 1;
-  double dmone = -1.0;
+  const double *mu = pol->mu;
+  const double dd = pol->dd;
+  const double cc = pol->cc;
+  const double ncc = -cc;
+  const int m = pol->deg;
+  const int one = 1;
+  const double dmone = -1.0;
   /*-------------------- pointers to v_[k-1],v_[k], v_[k+1] from w */
   double *vk   = w;
   double *vkp1 = w+n;
diff --git a/EVSL_1.1.1/SRC/chebsi.c b/EVSL_1.1.1/SRC/chebsi.c
index 84fee83..614919e 100644
--- a/EVSL_1.1.1/SRC/chebsi.c
+++ b/EVSL_1.1.1/SRC/chebsi.c
@@ -36,39 +36,42 @@
  *   @warning Memory allocation for Yo/lamo/reso within this function
  */
 
-int ChebSI(int nev, double *intv, int maxit,
-           double tol, double *vinit, polparams *pol, int *nevo,
+int ChebSI(int nev, const double * intv, int maxit,
+           double tol, const double * vinit, polparams const * const pol, int *nevo,
            double **lamo, double **Yo, double **reso, FILE *fstats) {
   /*-------------------- for stats */
   double tm, tall=0.0, tmv=0.0;
   int icol, nmv = 0;
-  double tolP = tol*0.01;
+  const double tolP = tol*0.01;
   tall = evsl_timer();
   //    int max_deg = pol->max_deg,   min_deg = pol->min_deg;
   /*-------------------   size of A */
-  int n = evsldata.n;
+  const int n = evsldata.n;
   /*--------------------   some constants frequently used */
-  char cT = 'T';
-  char cN = 'N';
-  int one = 1;
-  int nev2 = nev*nev;
-  int nnev = n*nev;
-  double done=1.0,dmone=-1.0,dzero=0.0;
+  const char cT = 'T';
+  const char cN = 'N';
+  const int one = 1;
+  const int nev2 = nev*nev;
+  const int nnev = n*nev;
+  const double done=1.0;
+  const double dmone=-1.0;
+  const double dzero=0.0;
   if (check_intv(intv, fstats) < 0) {
     *nevo = 0;
     return 0;
   }
-  double aa = intv[0];
-  double bb = intv[1];
+  const double aa = intv[0];
+  const double bb = intv[1];
   /*--------------------   frequency of convergence checks (check every cvcheck iterations) */
-  int cvcheck = 5;
+  const int cvcheck = 5;
   /*-------------------   misc */
   int i,j;
   // this code(step 1)is the same as in cheblan; to be moved
   // to a separate file for reuse
   /*-------------------- unpack some values from pol */
-  int deg =pol->deg;
-  double gamB=pol->gam, bar=pol->bar;
+  const int deg =pol->deg;
+  const double gamB=pol->gam;
+  const double bar=pol->bar;
   fprintf(fstats, "Cheb Poly- deg = %d, gam = %.15e, bar: %.15e\n",
       deg, gamB, bar);
   /*-------------------- gamB must be within [-1, 1] */
@@ -166,7 +169,7 @@ int ChebSI(int nev, double *intv, int maxit,
       evsl_dcopy(&nnact, buf, &one, PV+nnlock, &one);
       /*---  Compute active residuals R = PV - V*diag(evalT)    */
       for (i=nlock; i<nev; i++) {
-        double t = -evalT[i];
+        const double t = -evalT[i];
         //evsl_dscal(&n, &t, R+i*n, &one);
         for (j=0; j<n; j++) {
           R[i*n+j] = PV[i*n+j]+t*V[i*n+j];
@@ -178,18 +181,18 @@ int ChebSI(int nev, double *intv, int maxit,
       nact = 0;
       for (i=nlock; i<nev; i++) {
         /*---  Compute norm of R(:,i)   */
-        double resP = sqrt(evsl_ddot(&n, R+i*n, &one, R+i*n, &one));
+        const double resP = sqrt(evsl_ddot(&n, R+i*n, &one, R+i*n, &one));
         if (resP < tolP) {
           /*---  Compute norm of AV(:,i) - V(:,i)*Lambda(i)   */
           tm = evsl_timer();
           matvec_A(V+i*n, buf);
           tmv += evsl_timer() - tm;
           nmv++;
-          double rq = evsl_ddot(&n, V+i*n, &one, buf, &one);  // Rayleigh Quotient for A
+          const double rq = evsl_ddot(&n, V+i*n, &one, buf, &one);  // Rayleigh Quotient for A
           for (j=0; j < n; j++) {
             R[i*n+j] = buf[j] - rq*V[i*n+j];
           }
-          double resA = sqrt(evsl_ddot(&n, R+i*n, &one, R+i*n, &one));
+          const double resA = sqrt(evsl_ddot(&n, R+i*n, &one, R+i*n, &one));
           if (resA < tol) {
             lock_idx[nlock_new] = i;
             res[nlock+nlock_new] = resA;
@@ -248,7 +251,7 @@ int ChebSI(int nev, double *intv, int maxit,
   V_out = evsl_Malloc(nnev, double);
   int idx=0;
   for (i=0; i<nlock;i++) {
-    double t = Lam[i];
+    const double t = Lam[i];
     if ( t >= aa - EVSL_DBL_EPS_MULT * DBL_EPSILON && t <= bb + EVSL_DBL_EPS_MULT * DBL_EPSILON ) {
       Lam_out[idx] = t;
       res_out[idx] = res[i];
diff --git a/EVSL_1.1.1/SRC/dos_utils.c b/EVSL_1.1.1/SRC/dos_utils.c
index b24b67a..7900222 100644
--- a/EVSL_1.1.1/SRC/dos_utils.c
+++ b/EVSL_1.1.1/SRC/dos_utils.c
@@ -12,14 +12,14 @@
  *  @param[in] a Number to take reciprocal of
  *  @return reciprocal of a
  */
-double rec(const double a) { return 1.0 / a; }  // Reciprocal
+double rec(double a) { return 1.0 / a; }  // Reciprocal
 
 /**
  * @brief Inverse  square root
  * @param[in] a number to take inverse square root of
  * @return inverse square root of a
  * */
-double isqrt(const double a) { return 1.0 / sqrt(a); }  // Inverse square root
+double isqrt(double a) { return 1.0 / sqrt(a); }  // Inverse square root
 
 /**
  * Initalize BSolDataPol: use L-S polynomial to approximate a function 'ffun'
@@ -92,7 +92,7 @@ void FreeBSolPolData(BSolDataPol *data) {
  * @param[out] x p(A)b
  * @param[in, out] data Data to be cast to BSolDataPol
  */
-void BSolPol(double *b, double *x, void *data) {
+void BSolPol(const double *b, double *x, void *data) {
   BSolDataPol *pol = (BSolDataPol *)data;
   pnav(pol->mu, pol->deg, pol->cc, pol->dd, b, x, pol->wk);
 }
@@ -114,8 +114,8 @@ void BSolPol(double *b, double *x, void *data) {
  *
  *----------------------------------------------------------------------*/
 
-int apfun(const double c, const double h, const double *const xi,
-          double (*ffun)(double), const int npts, double *yi) {
+int apfun(double c, double h, const double *const xi,
+          double (*ffun)(double), int npts, double *yi) {
   int i = 0;
   for (i = 0; i < npts; i++) {
     yi[i] = ffun(c + h * xi[i]);
@@ -140,19 +140,20 @@ int apfun(const double c, const double h, const double *const xi,
  * @b Workspace
  * @param[in,out] w Work vector of length 3*n [allocate before call
  **/
-int pnav(double *mu, const int m, const double cc, const double dd, double *v,
-         double *y, double *w) {  // Really just ChebAv
-  int n = evsldata.n;
+int pnav(const double * mu, int m, double cc, const double dd,
+         const double * v, double *y, double *w) {  // Really just ChebAv
+  const int n = evsldata.n;
   /*-------------------- pointers to v_[k-1],v_[k], v_[k+1] from w */
   double *vk = w;
   double *vkp1 = w + n;
   double *vkm1 = vkp1 + n;
   /*-------------------- */
-  int k, one=1;
-  double t1 = 1.0 / dd;
-  double t2 = 2.0 / dd;
-  double ncc = -cc;
-  double dmone = -1.0;
+  int k;
+  const int one=1;
+  const double t1 = 1.0 / dd;
+  const double t2 = 2.0 / dd;
+  const double ncc = -cc;
+  const double dmone = -1.0;
   /*-------------------- vk <- v; vkm1 <- zeros(n,1) */
 #if 0
   /* LEAVE HERE IT FOR REFERENCE */
@@ -182,9 +183,9 @@ int pnav(double *mu, const int m, const double cc, const double dd, double *v,
   evsl_dscal(&n, mu, y, &one);
   /*-------------------- degree loop. k IS the degree */
   for (k = 1; k <= m; k++) {
-    double *v_cur = k == 1 ? v : vk;
-    double *v_old = k == 2 ? v : vkm1;
-    double t = k == 1 ? t1 : t2;
+    const double *v_cur = k == 1 ? v : vk;
+    const double *v_old = k == 2 ? v : vkm1;
+    const double t = k == 1 ? t1 : t2;
 
     matvec_B(v_cur, vkp1);
     evsl_daxpy(&n, &ncc, v_cur, &one, vkp1, &one);
@@ -228,7 +229,7 @@ int lsPol(double (*ffun)(double), BSolDataPol *pol) {
   pol->cc = (a + b) / 2;
   pol->dd = (b - a) / 2;
   /*------------------------- Number of points for Gauss-Chebyshev integration */
-  int maxDeg = pol->max_deg;
+  const int maxDeg = pol->max_deg;
   const int npts = maxDeg * 4;
   double *theti;
   theti = evsl_Malloc(npts, double);
@@ -307,4 +308,3 @@ int lsPol(double (*ffun)(double), BSolDataPol *pol) {
 
   return 0;
 }
-
diff --git a/EVSL_1.1.1/SRC/dumps.c b/EVSL_1.1.1/SRC/dumps.c
index 02ea98b..348dd6c 100644
--- a/EVSL_1.1.1/SRC/dumps.c
+++ b/EVSL_1.1.1/SRC/dumps.c
@@ -21,8 +21,8 @@
  * @param[in] a Values
  * @param[in] fn filename
  */
-void save_mtx_basic(int nrow, int ncol, int *ia,
-                    int *ja, double *a, const char *fn) {
+void save_mtx_basic(int nrow, int ncol, const int * ia,
+                    const int * ja, const double * a, const char *fn) {
   int i,j,nnz;
   FILE *fp = fopen(fn, "w");
 
@@ -43,7 +43,7 @@ void save_mtx_basic(int nrow, int ncol, int *ia,
  * @param[in] A csr matrix to save
  * @param[in] fn filename
  */
-void savemat(csrMat *A, const char *fn) {
+void savemat(const csrMat *A, const char *fn) {
   fprintf(stdout, " * saving a matrix into %s\n", fn);
   save_mtx_basic(A->nrows, A->ncols, A->ia, A->ja, A->a, fn);
 }
@@ -74,7 +74,7 @@ void save_vec(int n, const double *x, const char fn[]) {
  * @param[in] n num cols
  * @param[in] fn filename
  */
-void savedensemat(double *A, int lda, int m, int n, const char *fn) {
+void savedensemat(const double *A, int lda, int m, int n, const char *fn) {
   fprintf(stdout, " * saving a matrix into %s\n", fn);
   FILE *fp = fopen(fn, "w");
   int i,j;
diff --git a/EVSL_1.1.1/SRC/exDOS.c b/EVSL_1.1.1/SRC/exDOS.c
index b35a2e7..5236f08 100644
--- a/EVSL_1.1.1/SRC/exDOS.c
+++ b/EVSL_1.1.1/SRC/exDOS.c
@@ -14,7 +14,7 @@
  * @param[out] y y coordinates for dos plot
  * @note Both x and y are expected to be preallocated.
  */
-int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv) {
+int exDOS(const double *vals, int n, int npts, double *x, double *y, const double *intv) {
 /* sigma = coefficient for gaussian
    we select sigma so that at end of subinterval
    we have exp[-0.5 [H/sigma]^2] = 1/K. */
diff --git a/EVSL_1.1.1/SRC/lanTrbounds.c b/EVSL_1.1.1/SRC/lanTrbounds.c
index 227807f..eb01fd0 100644
--- a/EVSL_1.1.1/SRC/lanTrbounds.c
+++ b/EVSL_1.1.1/SRC/lanTrbounds.c
@@ -32,7 +32,7 @@
  * @return Returns 0 on success
  *
  **/
-int LanTrbounds(int lanm, int maxit, double tol, double *vinit,
+int LanTrbounds(int lanm, int maxit, double tol, const double *vinit,
                 int bndtype,
                 double *lammin, double *lammax, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
diff --git a/EVSL_1.1.1/SRC/lanbounds.c b/EVSL_1.1.1/SRC/lanbounds.c
index 5e4882b..ca86a92 100644
--- a/EVSL_1.1.1/SRC/lanbounds.c
+++ b/EVSL_1.1.1/SRC/lanbounds.c
@@ -18,7 +18,7 @@
  *    @param[out] lmax right bound of the desired interval containing all
  *    eigenvalues
  *----------------------------------------------------------------------*/
-int LanBounds(int msteps, double *v, double *lmin, double *lmax) {
+int LanBounds(int msteps, const double *v, double *lmin, double *lmax) {
   const int ifGenEv = evsldata.ifGenEv;
   double *alp, *bet, nbet, nalp, t, *V, *Z;
   int j, one=1, n;
diff --git a/EVSL_1.1.1/SRC/landos.c b/EVSL_1.1.1/SRC/landos.c
index fc16f5a..22e935d 100644
--- a/EVSL_1.1.1/SRC/landos.c
+++ b/EVSL_1.1.1/SRC/landos.c
@@ -34,8 +34,8 @@
  *
  *----------------------------------------------------------------------*/
 
-int LanDos(const int nvec, int msteps, int npts, double *xdos, double *ydos,
-           double *neig, const double *const intv) {
+int LanDos(int nvec, int msteps, int npts, double *xdos, double *ydos,
+           double *neig, const double *intv) {
   //--------------------
   const int ifGenEv = evsldata.ifGenEv;
   if (ifGenEv) {
diff --git a/EVSL_1.1.1/SRC/landosG.c b/EVSL_1.1.1/SRC/landosG.c
index ba323f4..ca090b8 100644
--- a/EVSL_1.1.1/SRC/landosG.c
+++ b/EVSL_1.1.1/SRC/landosG.c
@@ -38,7 +38,7 @@
  *----------------------------------------------------------------------*/
 
 int LanDosG(const int nvec, const int msteps, int npts, double *xdos, double *ydos,
-            double *neig, const double *const intv) {
+            double *neig, const double *intv) {
 
   int i, j, k;
   int maxit = msteps, m;  /* Max number of iterations */
diff --git a/EVSL_1.1.1/SRC/lapack/evsl_lapack.h b/EVSL_1.1.1/SRC/lapack/evsl_lapack.h
index afeb623..1056587 100644
--- a/EVSL_1.1.1/SRC/lapack/evsl_lapack.h
+++ b/EVSL_1.1.1/SRC/lapack/evsl_lapack.h
@@ -13,10 +13,10 @@ extern "C" {
 #endif
 #endif
 
-int evsl_dstev(const char *jobz, EVSL_Int *n, EVSL_Real *d__, EVSL_Real *e, EVSL_Real *z__, EVSL_Int *ldz, EVSL_Real *work, EVSL_Int *info);
-int evsl_dsyev(const char *jobz, const char *uplo, EVSL_Int *n, EVSL_Real *a, EVSL_Int *lda, EVSL_Real *w, EVSL_Real *work, EVSL_Int *lwork, EVSL_Int *info);
-int evsl_dstemr(const char *jobz, const char *range, EVSL_Int *n, EVSL_Real *d__, EVSL_Real *e, EVSL_Real *vl, EVSL_Real *vu, EVSL_Int *il, EVSL_Int *iu, EVSL_Int *m, EVSL_Real *w, EVSL_Real *z__, EVSL_Int *ldz, EVSL_Int *nzc, EVSL_Int *isuppz, EVSL_Int *tryrac, EVSL_Real *work, EVSL_Int *lwork, EVSL_Int *iwork, EVSL_Int *liwork, EVSL_Int *info);
-int evsl_zgesv(EVSL_Int *n, EVSL_Int *nrhs, EVSL_Complex *a, EVSL_Int *lda, EVSL_Int *ipiv, EVSL_Complex *b, EVSL_Int *ldb, EVSL_Int *info);
+int evsl_dstev(const char *jobz, const EVSL_Int *n, EVSL_Real *d__, EVSL_Real *e, EVSL_Real *z__, const EVSL_Int *ldz, EVSL_Real *work, EVSL_Int *info);
+int evsl_dsyev(const char *jobz, const char *uplo, const EVSL_Int *n, EVSL_Real *a, const EVSL_Int *lda, EVSL_Real *w, EVSL_Real *work, const EVSL_Int *lwork, EVSL_Int *info);
+int evsl_dstemr(const char *jobz, const char *range, const EVSL_Int *n, const EVSL_Real *d__, EVSL_Real *e, const EVSL_Real *vl, const EVSL_Real *vu, const EVSL_Int *il, const EVSL_Int *iu, EVSL_Int *m, EVSL_Real *w, EVSL_Real *z__, const EVSL_Int *ldz, const EVSL_Int *nzc, EVSL_Int *isuppz, EVSL_Int *tryrac, EVSL_Real *work, const EVSL_Int *lwork, EVSL_Int *iwork, const EVSL_Int *liwork, EVSL_Int *info);
+int evsl_zgesv(const EVSL_Int *n, const EVSL_Int *nrhs, EVSL_Complex *a, const EVSL_Int *lda, EVSL_Int *ipiv, EVSL_Complex *b, const EVSL_Int *ldb, EVSL_Int *info);
 
 #ifndef EVSL_USING_EVSL_LAPACK
 #ifdef __cplusplus
diff --git a/EVSL_1.1.1/SRC/misc_la.c b/EVSL_1.1.1/SRC/misc_la.c
index a4571a3..03ee17e 100644
--- a/EVSL_1.1.1/SRC/misc_la.c
+++ b/EVSL_1.1.1/SRC/misc_la.c
@@ -34,9 +34,9 @@ int SymmTridEig(double *eigVal, double *eigVec, int n,
                 const double *diag, const double *sdiag) {
   double tms = evsl_timer();
   // compute eigenvalues and eigenvectors or eigvalues only
-  char jobz = eigVec ? 'V' : 'N';
-  int nn = n;
-  int ldz = n;
+  const char jobz = eigVec ? 'V' : 'N';
+  const int nn = n;
+  const int ldz = n;
   int info;  // output flag
   // copy diagonal and subdiagonal elements to alp and bet
   double *alp = eigVal;
@@ -90,9 +90,9 @@ int SymmTridEig(double *eigVal, double *eigVec, int n,
  * ----------------------------------------------------------------------- */
 int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu,
                  int *nevO, const double *diag, const double *sdiag) {
-  double tms = evsl_timer();
-  char jobz = 'V';  // compute eigenvalues and eigenvectors
-  char range = 'V'; // compute eigenvalues in an interval
+  const double tms = evsl_timer();
+  const char jobz = 'V';  // compute eigenvalues and eigenvectors
+  const char range = 'V'; // compute eigenvalues in an interval
 
   // this does not use mwlapack for mex files
   int info;
@@ -102,11 +102,11 @@ int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu,
   isuppz = evsl_Malloc(2*n, int);
   //-------------------- real work array
   double *work;
-  int lwork = 18*n;
+  const int lwork = 18*n;
   work = evsl_Malloc(lwork, double);
   //-------------------- int work array
   int *iwork;
-  int liwork = 10*n;
+  const int liwork = 10*n;
   iwork = evsl_Calloc(liwork, int);
   //-------------------- copy diagonal + subdiagonal elements
   //                     to alp and bet
@@ -157,13 +157,13 @@ int SymmTridEigS(double *eigVal, double *eigVec, int n, double vl, double vu,
  *     @param[in] ldq Leading dimension q
  *     @param[out] lam Eigenvalues
  *- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
-void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam) {
+void SymEigenSolver(int n, const double *A, int lda, double *Q, int ldq, double *lam) {
   double tms = evsl_timer();
   /* compute eigenvalues/vectors of A that n x n, symmetric
    * eigenvalues saved in lam: the eigenvalues in ascending order
    * eigenvectors saved in Q */
-  char jobz='V';/* want eigenvectors */
-  char uplo='U';/* use upper triangular part of the matrix */
+  const char jobz='V';/* want eigenvectors */
+  const char uplo='U';/* use upper triangular part of the matrix */
   /*   copy A to Q */
   int i,j;
   for (i=0; i<n; i++) {
@@ -203,9 +203,9 @@ void SymEigenSolver(int n, double *A, int lda, double *Q, int ldq, double *lam)
  * @param[out] v Output
  * @param[out] w Output
  **/
-void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, double *w) {
-  double tms = evsl_timer();
-  double eta = 1.0 / sqrt(2.0);
+void CGS_DGKS(int n, int k, int i_max, const double *Q, double *v, double *nrmv, double *w) {
+  const double tms = evsl_timer();
+  const double eta = 1.0 / sqrt(2.0);
   int i, one=1;
 #if USE_DGEMV
   char cT = 'T', cN = 'N';
@@ -221,7 +221,7 @@ void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, doubl
 #else
     int j;
     for (j=0; j<k; j++) {
-       double t = -evsl_ddot(&n, &Q[j*n], &one, v, &one);
+       const double t = -evsl_ddot(&n, &Q[j*n], &one, v, &one);
        evsl_daxpy(&n, &t, &Q[j*n], &one, v, &one);
     }
 #endif
@@ -235,7 +235,7 @@ void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, doubl
   if (nrmv) {
     *nrmv = new_nrm;
   }
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_reorth += tme - tms;
 }
 
@@ -251,13 +251,17 @@ void CGS_DGKS(int n, int k, int i_max, double *Q, double *v, double *nrmv, doubl
  * @param[out] v Output
  * @param[out] w Output
  **/
-void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q,
+void CGS_DGKS2(int n, int k, int i_max, const double *Z, const double *Q,
                double *v, double *w) {
-  double tms = evsl_timer();
-  int i, one=1;
+  const double tms = evsl_timer();
+  int i;
+  const int one=1;
 #if USE_DGEMV
-  char cT = 'T', cN = 'N';
-  double done=1.0, dmone=-1.0, dzero=0.0;
+  const char cT = 'T';
+  const char cN = 'N';
+  const double done=1.0;
+  const double dmone=-1.0;
+  const double dzero=0.0;
 #endif
   for (i=0; i<i_max; i++) {
 #if USE_DGEMV
@@ -266,12 +270,12 @@ void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q,
 #else
     int j;
     for (j=0; j<k; j++) {
-       double t = -evsl_ddot(&n, &Q[j*n], &one, v, &one);
+       const double t = -evsl_ddot(&n, &Q[j*n], &one, v, &one);
        evsl_daxpy(&n, &t, &Z[j*n], &one, v, &one);
     }
 #endif
   }
-  double tme = evsl_timer();
+  const double tme = evsl_timer();
   evslstat.t_reorth += tme - tms;
 }
 
@@ -287,17 +291,17 @@ void CGS_DGKS2(int n, int k, int i_max, double *Z, double *Q,
  *
  * @warning Aliasing happens in call to CGS_DGKS
  */
-void orth(double *V, int n, int k, double *Vo, double *work) {
+void orth(const double *V, int n, int k, double *Vo, double *work) {
   int i;
-  int one=1;
-  int nk = n*k;
+  const int one=1;
+  const int nk = n*k;
   evsl_dcopy(&nk, V, &one, Vo, &one);
-  double tt = evsl_ddot(&n, Vo, &one, Vo, &one);
+  const double tt = evsl_ddot(&n, Vo, &one, Vo, &one);
   double nrmv = sqrt(tt);
   double t = 1.0 / nrmv;
   evsl_dscal(&n, &t, Vo, &one);
   for (i = 1; i < k; i++) {
-    int istart = i*n;
+    const int istart = i*n;
     CGS_DGKS(n, i, NGS_MAX, Vo, Vo+istart, &nrmv, work);
     t = 1.0 / nrmv;
     evsl_dscal(&n, &t, Vo+istart, &one);
diff --git a/EVSL_1.1.1/SRC/ratfilter.c b/EVSL_1.1.1/SRC/ratfilter.c
index c268d48..5fb3829 100644
--- a/EVSL_1.1.1/SRC/ratfilter.c
+++ b/EVSL_1.1.1/SRC/ratfilter.c
@@ -62,8 +62,8 @@ void contQuad(int method, int n, EVSL_Complex* zk) {
    * @param[out] xx    : function values at real locations z
    *
    *----------------------------------------------------------------------*/
-void ratf2p2(int n, int *mulp, EVSL_Complex *zk, EVSL_Complex* alp, int m,
-             double *z, double *xx) {
+void ratf2p2(int n, const int *mulp, const EVSL_Complex *zk, const EVSL_Complex* alp, int m,
+             const double *z, double *xx) {
   EVSL_Complex y, x, t;
   int i, j, k, k1, k2;
   for (k2=0; k2<m; k2++) {
@@ -162,13 +162,13 @@ EVSL_Complex integg2(EVSL_Complex s1, EVSL_Complex s2,
  * @param[out] omega LS weight for each pole
  *
  *----------------------------------------------------------------------*/
-void weights(int n, EVSL_Complex* zk, int* mulp, double lambda,
+void weights(int n, const EVSL_Complex* zk, const int* mulp, double lambda,
              EVSL_Complex* omega) {
   int INFO;
-  int nrhs = 1;
+  const int nrhs = 1;
   int *ipiv;
   int m;
-  double mu = 10.0;
+  const double mu = 10.0;
   int i, j, ii, jj, ki, kj, n1, n2, nf=0, k1, k2;
   EVSL_Complex s1, s2, s3, t;
   EVSL_Complex *rhs, *A, *B, *mat, *alp, *bet;
@@ -306,7 +306,7 @@ void weights(int n, EVSL_Complex* zk, int* mulp, double lambda,
    * @param[out] omegaM: multiple LS weights
    *
    *----------------------------------------------------------------------*/
-int scaleweigthts(int n, double a, double b, EVSL_Complex *zk, int* mulp,
+int scaleweigthts(int n, double a, double b, EVSL_Complex *zk, const int* mulp,
                   EVSL_Complex* omegaM) {
   int i, j, k, nf=0;
   double c, h;
@@ -367,7 +367,7 @@ void set_ratf_def(ratparams *rat) {
  *    cc  : ..center of interval\n
  *
  *--------------------------------------------------------------------*/
-int find_ratf(double *intv, ratparams *rat) {
+int find_ratf(const double *intv, ratparams *rat) {
   EVSL_Complex *omega; // weights of the poles
   EVSL_Complex *zk;    // location of the poles
   int *mulp;             // multiplicity of the each pole
@@ -440,16 +440,16 @@ void free_rat(ratparams *rat) {
  * @param[out] x Becomes R(A)b
  *
  */
-void RatFiltApply(int n, ratparams *rat, double *b, double *x, double *w6) {
+void RatFiltApply(int n, const ratparams *rat, const double *b, double *x, double *w6) {
   double tt = evsl_timer();
   const int ifGenEv = evsldata.ifGenEv;
   int jj, kk, k=0, kf;
   int *mulp = rat->mulp;
   int num = rat->num;
   EVSL_Complex *omega = rat->omega;
-  double dtwo = 2.0;
-  double done = 1.0;
-  int one = 1;
+  const double dtwo = 2.0;
+  const double done = 1.0;
+  const int one = 1;
 
   double *xr, *xz, *bz, *br, *yr=NULL, *yz=NULL;
   double zkr, zkc;
diff --git a/EVSL_1.1.1/SRC/ratlanNr.c b/EVSL_1.1.1/SRC/ratlanNr.c
index 019512c..78a315a 100644
--- a/EVSL_1.1.1/SRC/ratlanNr.c
+++ b/EVSL_1.1.1/SRC/ratlanNr.c
@@ -44,7 +44,7 @@
  *  @param[out] fstats   File stream which stats are printed to
  *
  * ------------------------------------------------------------ */
-int RatLanNr(double *intv, int maxit, double tol, double *vinit,
+int RatLanNr(const double *intv, int maxit, double tol, const double *vinit,
              ratparams *rat, int *nevOut, double **lamo, double **Wo,
              double **reso, FILE *fstats) {
   //-------------------- to report timings/
@@ -60,15 +60,16 @@ int RatLanNr(double *intv, int maxit, double tol, double *vinit,
     do_print = 0;
   }
   /*--------------------   frequently used constants  */
-  char cN = 'N';
-  int one = 1;
-  double done = 1.0, dzero = 0.0;
+  const char cN = 'N';
+  const int one = 1;
+  const double done = 1.0;
+  const double dzero = 0.0;
   /*--------------------   Ntest = when to start testing convergence */
-  int Ntest = 30;
+  const int Ntest = 30;
   /*--------------------   how often to test */
-  int cycle = 20;
+  const int cycle = 20;
   /* size of the matrix */
-  int n = evsldata.n;
+  const int n = evsldata.n;
   /* max num of its */
   maxit = evsl_min(n, maxit);
   /*-------------------- Caveat !!!: To prevent integer overflow, we save
@@ -81,14 +82,14 @@ int RatLanNr(double *intv, int maxit, double tol, double *vinit,
   /*-------------------- Rational filter with pole at ((a+b)/2,(b-a)/2) with
     multiplicity pow, bar value equals 1/2        */
   /*-------------------- a, b, used for testing only at end */
-  double bar = 0.5;
+  const double bar = 0.5;
   if (check_intv(intv, fstats) < 0) {
     *nevOut = 0;
     *lamo = NULL; *Wo = NULL; *reso = NULL;
     return 0;
   }
-  double aa = intv[0];
-  double bb = intv[1];
+  const double aa = intv[0];
+  const double bb = intv[1];
   //int deg = rat->pow; // multiplicity of the pole
   /*-----------------------------------------------------------------------*
    * *Non-restarted* Lanczos iteration
diff --git a/EVSL_1.1.1/SRC/ratlanTr.c b/EVSL_1.1.1/SRC/ratlanTr.c
index af871d4..103caeb 100644
--- a/EVSL_1.1.1/SRC/ratlanTr.c
+++ b/EVSL_1.1.1/SRC/ratlanTr.c
@@ -53,8 +53,8 @@
  * @warning memory allocation for W/vals/resW within this function
  *
  **/
-int RatLanTr(int lanm, int nev, double *intv, int maxit,
-             double tol, double *vinit, ratparams *rat, int *nev2,
+int RatLanTr(int lanm, int nev, const double *intv, int maxit,
+             double tol, const double *vinit, ratparams *rat, int *nev2,
              double **vals, double **W, double **resW, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
   /*-------------------- for stats */
@@ -68,11 +68,11 @@ int RatLanTr(int lanm, int nev, double *intv, int maxit,
     do_print = 0;
   }
   /* size of the matrix */
-  int n = evsldata.n;
-  size_t n_l = n;
+  const int n = evsldata.n;
+  const size_t n_l = n;
   /*--------------------- adjust lanm and maxit */
   lanm = evsl_min(lanm, n);
-  int lanm1=lanm+1;
+  const int lanm1=lanm+1;
   size_t lanm1_l = lanm1;
   /*  if use full lanczos, should not do more than n iterations */
   if (lanm == n) {
@@ -86,13 +86,15 @@ int RatLanTr(int lanm, int nev, double *intv, int maxit,
   //int evFrac = 2;
   /*--------------------   some constants frequently used */
   /* char cT='T'; */
-  char cN = 'N';
-  int one = 1;
-  double done=1.0, dmone=-1.0, dzero=0.0;
+  const char cN = 'N';
+  const int one = 1;
+  const double done=1.0;
+  const double dmone=-1.0;
+  const double dzero=0.0;
   /*-------------------- Ntest = when to start testing convergence */
   int Ntest = evsl_min(lanm, nev+50);
   /*--------------------   how often to test */
-  int cycle = 30;
+  const int cycle = 30;
   int i, ll, /* count, last_count,*/ jl, last_jl;
   /*-----------------------------------------------------------------------
     -----------------------------------------------------------------------*/
@@ -102,10 +104,10 @@ int RatLanTr(int lanm, int nev, double *intv, int maxit,
     *vals = NULL; *W = NULL; *resW = NULL;
     return 0;
   }
-  double aa = intv[0];
-  double bb = intv[1];
+  const double aa = intv[0];
+  const double bb = intv[1];
   //int deg = rat->pow;
-  double bar = 0.5; // for the scaled rational filter
+  const double bar = 0.5; // for the scaled rational filter
   /*-----------------------------------------------------------------------*
    * *thick restarted* Lanczos step
    *-----------------------------------------------------------------------*/
diff --git a/EVSL_1.1.1/SRC/simpson.c b/EVSL_1.1.1/SRC/simpson.c
index f04e86e..ff09f4d 100644
--- a/EVSL_1.1.1/SRC/simpson.c
+++ b/EVSL_1.1.1/SRC/simpson.c
@@ -22,7 +22,7 @@
  *    @note In-place version.
  *----------------------------------------------------------------------*/
 
-void simpson(double* xi, double* yi, int npts) {
+void simpson(const double* xi, double* yi, int npts) {
   double tc = 0.0, ti, tm, ysav;
   int i = 1;
   //-------------------- save yi[even_i]
diff --git a/EVSL_1.1.1/SRC/spmat.c b/EVSL_1.1.1/SRC/spmat.c
index 244b0b5..9699f09 100644
--- a/EVSL_1.1.1/SRC/spmat.c
+++ b/EVSL_1.1.1/SRC/spmat.c
@@ -123,7 +123,7 @@ void free_csr(csrMat *csr) {
  *  @param[out] B Destination matrix
  *  @param[in] allocB Whether or not to allocate B
  */
-void csr_copy(csrMat *A, csrMat *B, int allocB) {
+void csr_copy(const csrMat *A, csrMat *B, int allocB) {
   int nrows = A->nrows;
   int ncols = A->ncols;
   int nnz = A->ia[nrows];
@@ -156,7 +156,7 @@ void free_coo(cooMat *coo) {
  * @param[in] coo COO matrix
  * @param[out] csr CSR matrix
  */
-int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr) {
+int cooMat_to_csrMat(int cooidx, const cooMat *coo, csrMat *csr) {
   const int nnz = coo->nnz;
   //printf("@@@@ coo2csr, nnz %d\n", nnz);
   /* allocate memory */
@@ -201,7 +201,7 @@ int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr) {
  * @param[in] a Values
  * @param[out] A output CSR matrix
  */
-int arrays_copyto_csrMat(int nrow, int ncol, int *ia, int *ja, double *a,
+int arrays_copyto_csrMat(int nrow, int ncol, const int *ia, const int *ja, const double *a,
                          csrMat *A) {
   int nnz = ia[nrow];
 
@@ -251,8 +251,8 @@ double dcsrinfnrm(csrMat *A){
  * @param[in] x Input vector
  * @param[out] y Output vector
  */
-void dcsrmv(char trans, int nrow, int ncol, double *a,
-    int *ia, int *ja, double *x, double *y) {
+void dcsrmv(char trans, int nrow, int ncol, const double *a,
+    const int *ia, const int *ja, const double *x, double *y) {
   int  len, jj=nrow;
   if (trans == 'N') {
     //#pragma omp parallel for schedule(guided)
@@ -303,7 +303,7 @@ void dcsrmv(char trans, int nrow, int ncol, double *a,
 * @param[out] y Output vector
 * @param[in] data CSR matrix
 */
-void matvec_csr(double *x, double *y, void *data) {
+void matvec_csr(const double *x, double *y, void *data) {
   csrMat *A = (csrMat *) data;
 #ifdef EVSL_USING_INTEL_MKL
   char cN = 'N';
@@ -322,7 +322,7 @@ void matvec_csr(double *x, double *y, void *data) {
 /** @brief inline function used by matadd
  * insert an element pointed by j of A (times t) to location k in C (row i)
  * */
-inline void matadd_insert(double t, csrMat *A, csrMat *C, int i, int *k,
+inline void matadd_insert(double t, const csrMat *A, csrMat *C, int i, int *k,
                           int *j, int *map) {
   /* if this entry already exists in C:
    * checking if it is the first entry of this row
@@ -363,7 +363,7 @@ inline void matadd_insert(double t, csrMat *A, csrMat *C, int i, int *k,
  * entry in C.ja and C.a for entry in A.ja[i] and A.a[i]
  * @param[out] mapB the same as mapA
  * */
-int matadd(double alp, double bet, csrMat *A, csrMat *B, csrMat *C,
+int matadd(double alp, double bet, const csrMat *A, const csrMat *B, csrMat *C,
            int *mapA, int *mapB) {
   int nnzA, nnzB, i, jA, jB, k;
   /* check dimension */
@@ -432,7 +432,7 @@ int speye(int n, csrMat *A) {
  * @param[in,out] A Coo Matrix to scale
  * @param[in] d The vector that contains d(i)
  */
-void diagScalCoo(cooMat *A, double *d) {
+void diagScalCoo(cooMat *A, const double *d) {
   int i, row, col, nnz = A->nnz;
   /* diagonal scaling for A */
   for (i=0; i<nnz; i++) {
@@ -448,7 +448,7 @@ void diagScalCoo(cooMat *A, double *d) {
  * @param[in,out] A CSR Matrix to scale
  * @param[in] d The vector that contains d(i)
  */
-void diagScalCsr(csrMat *A, double *d) {
+void diagScalCsr(csrMat *A, const double *d) {
   int i, j;
   /* diagonal scaling for A */
   for (i=0; i<A->nrows; i++) {
@@ -464,7 +464,7 @@ void diagScalCsr(csrMat *A, double *d) {
  * @param[in]  B Matrix to extract the diagonal
  * @param[out] d preallocated vector of lengeth B.nrows
  */
-void extrDiagCsr(csrMat *B, double *d) {
+void extrDiagCsr(const csrMat *B, double *d) {
   int i, j, nrows = B->nrows;
   for (i=0; i<nrows; i++) {
     d[i] = 0.0;
@@ -482,7 +482,7 @@ void extrDiagCsr(csrMat *B, double *d) {
  * @param[in]  A Matrix
  * @param[out] U Matrix
  */
-void triuCsr(csrMat *A, csrMat *U) {
+void triuCsr(const csrMat *A, csrMat *U) {
   EVSL_Int i, j;
   EVSL_Unsigned k = 0, nnzu = 0;
 
diff --git a/EVSL_1.1.1/SRC/spslice.c b/EVSL_1.1.1/SRC/spslice.c
index 3f87f58..98d259c 100644
--- a/EVSL_1.1.1/SRC/spslice.c
+++ b/EVSL_1.1.1/SRC/spslice.c
@@ -26,7 +26,7 @@
  * @param[out] ecnt estimated num of eigenvalues in the interval of interest
  *
  *----------------------------------------------------------------------*/
-int kpmdos(int Mdeg, int damping, int nvec, double *intv,
+int kpmdos(int Mdeg, int damping, int nvec, const double *intv,
     double *mu, double *ecnt) {
   /*-------------------- initialize variables */
   const int ifGenEv = evsldata.ifGenEv;
@@ -148,7 +148,7 @@ int kpmdos(int Mdeg, int damping, int nvec, double *intv,
  *  where p(t) is the approximate DOS as given in the KPM method
  *  in the expanded form:  \f$\sum mu_i C_i /\sqrt{1-t^2}\f$
  **/
-void intChx(const int Mdeg, double *mu, const int npts, double *xi, double *yi) {
+void intChx(const int Mdeg, const double *mu, const int npts, const double *xi, double *yi) {
   //
   int ndp1, j, k;
   if(npts <= 0) {
@@ -198,7 +198,7 @@ void intChx(const int Mdeg, double *mu, const int npts, double *xi, double *yi)
  *                of eigenvalues in the interval [a b] (input).
  *
  *----------------------------------------------------------------------*/
-int spslicer(double *sli, double *mu, int Mdeg, double *intv, int n_int, int npts) {
+int spslicer(double *sli, const double *mu, int Mdeg, const double *intv, int n_int, int npts) {
   int ls, ii, err=0;
   double  ctr, wid, aL, bL, target, aa, bb;
 
diff --git a/EVSL_1.1.1/SRC/spslice2.c b/EVSL_1.1.1/SRC/spslice2.c
index fa7b2ac..4041d78 100644
--- a/EVSL_1.1.1/SRC/spslice2.c
+++ b/EVSL_1.1.1/SRC/spslice2.c
@@ -22,7 +22,7 @@
  *
  *----------------------------------------------------------------------*/
 
-void spslicer2(double* xi, double* yi, int n_int, int npts, double* sli) {
+void spslicer2(const double* xi, double* yi, int n_int, int npts, double* sli) {
   /*-------------------- makes a call here to  integration by Simpson */
   double want;
   int k = 0;
diff --git a/EVSL_1.1.1/SRC/vect.c b/EVSL_1.1.1/SRC/vect.c
index 5fd53d9..8a384b0 100644
--- a/EVSL_1.1.1/SRC/vect.c
+++ b/EVSL_1.1.1/SRC/vect.c
@@ -149,7 +149,7 @@ void sort_double(int n, double *v, int *ind) {
  * @param[in] x Input vector
  * @param[out] y Output vector
  * */
-void vec_perm(int n, int *p, double *x, double *y) {
+void vec_perm(int n, const int *p, const double *x, double *y) {
   if (!p) {
     memcpy(y, x, n*sizeof(double));
   } else {
@@ -166,7 +166,7 @@ void vec_perm(int n, int *p, double *x, double *y) {
  * @param[in] p Permutation vector
  * @param[in] x Input vector
  * @param[in] y Output vector */
-void vec_iperm(int n, int *p, double *x, double *y) {
+void vec_iperm(int n, const int *p, const double *x, double *y) {
   if (!p) {
     memcpy(y, x, n*sizeof(double));
   } else {

From 4266a06612d5725fc33a2f5b10fdc148b2a3c331 Mon Sep 17 00:00:00 2001
From: Luke Erlandson <LukeAErlandson@gmail.com>
Date: Mon, 26 Nov 2018 21:40:05 -0500
Subject: [PATCH 2/5] Improve consistency of consts

---
 EVSL_1.1.1/INC/evsl.h            |  4 ++--
 EVSL_1.1.1/INC/internal_header.h | 10 +++++-----
 EVSL_1.1.1/INC/struct.h          |  4 ++--
 EVSL_1.1.1/SRC/cheblanNr.c       |  2 +-
 EVSL_1.1.1/SRC/cheblanTr.c       |  2 +-
 EVSL_1.1.1/SRC/chebsi.c          |  2 +-
 6 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/EVSL_1.1.1/INC/evsl.h b/EVSL_1.1.1/INC/evsl.h
index 3118415..128e8bd 100644
--- a/EVSL_1.1.1/INC/evsl.h
+++ b/EVSL_1.1.1/INC/evsl.h
@@ -76,10 +76,10 @@
 )
 
 /*- - - - - - - - - cheblanNr.c */
-int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit, polparams const * const pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
+int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit, const polparams* pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
 
 /*- - - - - - - - - cheblanTr.c */
-int ChebLanTr(int lanm, int nev, const double * intv, int maxit, double tol, const double * vinit, polparams const * const pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
+int ChebLanTr(int lanm, int nev, const double * intv, int maxit, double tol, const double * vinit, const polparams *pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
 
 /*- - - - - - - - - chebpoly.c */
 void set_pol_def(polparams *pol);
diff --git a/EVSL_1.1.1/INC/internal_header.h b/EVSL_1.1.1/INC/internal_header.h
index 663a334..4f4f3de 100644
--- a/EVSL_1.1.1/INC/internal_header.h
+++ b/EVSL_1.1.1/INC/internal_header.h
@@ -98,7 +98,7 @@ static inline void matvec_A(const double *x, double *y) {
 * @brief y = B * x
 * This is the matvec function for the matrix B in evsldata
 */
-static inline void matvec_B(double const * const x, double *y) {
+static inline void matvec_B(const double *x, double *y) {
   CHKERR(!evsldata.Bmv);
   const double tms = evsl_timer();
   evsldata.Bmv->func(x, y, evsldata.Bmv->data);
@@ -111,7 +111,7 @@ static inline void matvec_B(double const * const x, double *y) {
 * @brief y = B \ x
 * This is the solve function for the matrix B in evsldata
 */
-static inline void solve_B(double const * const x, double *y) {
+static inline void solve_B(const double *x, double *y) {
   CHKERR(!evsldata.Bsol);
   const double tms = evsl_timer();
   evsldata.Bsol->func(x, y, evsldata.Bsol->data);
@@ -124,7 +124,7 @@ static inline void solve_B(double const * const x, double *y) {
 * @brief y = LT \ x or y = SQRT(B) \ x
 * This is the solve function for the matrix B in evsldata
 */
-static inline void solve_LT(double const * const x, double *y) {
+static inline void solve_LT(const double *x, double *y) {
   CHKERR(!evsldata.LTsol);
   const double tms = evsl_timer();
   evsldata.LTsol->func(x, y, evsldata.LTsol->data);
@@ -138,7 +138,7 @@ static inline void solve_LT(double const * const x, double *y) {
 * This is the solve function for the complex shifted matrix
 */
 static inline void solve_ASigB(EVSLASIGMABSol *sol, int n,
-                               double const * const br, double const * const bz,
+                               const double *br, const double *bz,
                                double *xr, double *xz) {
   const double tms = evsl_timer();
   (sol->func)(n, br, bz, xr, xz, sol->data);
@@ -150,7 +150,7 @@ static inline void solve_ASigB(EVSLASIGMABSol *sol, int n,
 /**
  * @brief check if an interval is valid
  * */
-static inline int check_intv(double const * const intv, FILE *fstats) {
+static inline int check_intv(const double *intv, FILE *fstats) {
   /* intv[4]: ( intv[0], intv[1] ) is the inteval of interest
    *          ( intv[2], intv[3] ) is the spectrum bounds
    * return   0: ok
diff --git a/EVSL_1.1.1/INC/struct.h b/EVSL_1.1.1/INC/struct.h
index 05a1165..350e4af 100644
--- a/EVSL_1.1.1/INC/struct.h
+++ b/EVSL_1.1.1/INC/struct.h
@@ -80,12 +80,12 @@ typedef struct _polparams {
  * be the  solution (complex vector),  and "data" contains  all the
  * data  needed  by  the  solver.
  */
-typedef void (*SolFuncC)(int n, double const * const br, double const * const bz, double *xr, double *xz, void *data);
+typedef void (*SolFuncC)(int n, const double *br, const double *bz, double *xr, double *xz, void *data);
 
 /**
  * @brief function prototype for applying the solve B x = b
  */
-typedef void (*SolFuncR)(double const * const b, double *x, void *data);
+typedef void (*SolFuncR)(const double *b, double *x, void *data);
 
 /**
  * @brief matvec function prototype
diff --git a/EVSL_1.1.1/SRC/cheblanNr.c b/EVSL_1.1.1/SRC/cheblanNr.c
index 8d773b4..cac46b4 100644
--- a/EVSL_1.1.1/SRC/cheblanNr.c
+++ b/EVSL_1.1.1/SRC/cheblanNr.c
@@ -55,7 +55,7 @@
  *
  * ------------------------------------------------------------ */
 int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit,
-              polparams const * const pol, int *nevOut, double **lamo, double **Wo,
+              const polparams *pol, int *nevOut, double **lamo, double **Wo,
               double **reso, FILE *fstats) {
   //-------------------- to report timings/
   double tall, tm1 = 0.0, tt;
diff --git a/EVSL_1.1.1/SRC/cheblanTr.c b/EVSL_1.1.1/SRC/cheblanTr.c
index 8bba4c5..bac934a 100644
--- a/EVSL_1.1.1/SRC/cheblanTr.c
+++ b/EVSL_1.1.1/SRC/cheblanTr.c
@@ -57,7 +57,7 @@
  *
  **/
 int ChebLanTr(int lanm, int nev, const double * intv, int maxit,
-              double tol, const double * vinit, polparams const * const pol, int *nev2,
+              double tol, const double * vinit, const polparams* const pol, int *nev2,
               double **vals, double **W, double **resW, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
   /*-------------------- for stats */
diff --git a/EVSL_1.1.1/SRC/chebsi.c b/EVSL_1.1.1/SRC/chebsi.c
index 614919e..2e229a5 100644
--- a/EVSL_1.1.1/SRC/chebsi.c
+++ b/EVSL_1.1.1/SRC/chebsi.c
@@ -37,7 +37,7 @@
  */
 
 int ChebSI(int nev, const double * intv, int maxit,
-           double tol, const double * vinit, polparams const * const pol, int *nevo,
+           double tol, const double * vinit, const polparams *pol, int *nevo,
            double **lamo, double **Yo, double **reso, FILE *fstats) {
   /*-------------------- for stats */
   double tm, tall=0.0, tmv=0.0;

From 7bc8f663278a86972f498716af516355bc030643 Mon Sep 17 00:00:00 2001
From: Luke Erlandson <LukeAErlandson@gmail.com>
Date: Mon, 26 Nov 2018 21:54:09 -0500
Subject: [PATCH 3/5] Further improve consistency of consts x *

---
 EVSL_1.1.1/INC/evsl.h      | 4 ++--
 EVSL_1.1.1/SRC/cheblanNr.c | 2 +-
 EVSL_1.1.1/SRC/cheblanTr.c | 4 ++--
 EVSL_1.1.1/SRC/chebpoly.c  | 8 ++++----
 EVSL_1.1.1/SRC/chebsi.c    | 4 ++--
 EVSL_1.1.1/SRC/dos_utils.c | 4 ++--
 EVSL_1.1.1/SRC/dumps.c     | 4 ++--
 7 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/EVSL_1.1.1/INC/evsl.h b/EVSL_1.1.1/INC/evsl.h
index 128e8bd..5c1ce59 100644
--- a/EVSL_1.1.1/INC/evsl.h
+++ b/EVSL_1.1.1/INC/evsl.h
@@ -76,10 +76,10 @@
 )
 
 /*- - - - - - - - - cheblanNr.c */
-int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit, const polparams* pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
+int ChebLanNr(const double *intv, int maxit, double tol, const double *vinit, const polparams* pol, int *nevOut, double **lamo, double **Wo, double **reso, FILE *fstats);
 
 /*- - - - - - - - - cheblanTr.c */
-int ChebLanTr(int lanm, int nev, const double * intv, int maxit, double tol, const double * vinit, const polparams *pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
+int ChebLanTr(int lanm, int nev, const double *intv, int maxit, double tol, const double *vinit, const polparams *pol, int *nev2, double **vals, double **W, double **resW, FILE *fstats);
 
 /*- - - - - - - - - chebpoly.c */
 void set_pol_def(polparams *pol);
diff --git a/EVSL_1.1.1/SRC/cheblanNr.c b/EVSL_1.1.1/SRC/cheblanNr.c
index cac46b4..59da055 100644
--- a/EVSL_1.1.1/SRC/cheblanNr.c
+++ b/EVSL_1.1.1/SRC/cheblanNr.c
@@ -54,7 +54,7 @@
  * @warning memory allocation for Wo/lamo/reso within this function
  *
  * ------------------------------------------------------------ */
-int ChebLanNr(const double * intv, int maxit, double tol, const double * vinit,
+int ChebLanNr(const double *intv, int maxit, double tol, const double *vinit,
               const polparams *pol, int *nevOut, double **lamo, double **Wo,
               double **reso, FILE *fstats) {
   //-------------------- to report timings/
diff --git a/EVSL_1.1.1/SRC/cheblanTr.c b/EVSL_1.1.1/SRC/cheblanTr.c
index bac934a..bf10b24 100644
--- a/EVSL_1.1.1/SRC/cheblanTr.c
+++ b/EVSL_1.1.1/SRC/cheblanTr.c
@@ -56,8 +56,8 @@
  * @warning memory allocation for W/vals/resW within this function
  *
  **/
-int ChebLanTr(int lanm, int nev, const double * intv, int maxit,
-              double tol, const double * vinit, const polparams* const pol, int *nev2,
+int ChebLanTr(int lanm, int nev, const double *intv, int maxit,
+              double tol, const double *vinit, const polparams* const pol, int *nev2,
               double **vals, double **W, double **resW, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
   /*-------------------- for stats */
diff --git a/EVSL_1.1.1/SRC/chebpoly.c b/EVSL_1.1.1/SRC/chebpoly.c
index 47e9740..32cacb7 100644
--- a/EVSL_1.1.1/SRC/chebpoly.c
+++ b/EVSL_1.1.1/SRC/chebpoly.c
@@ -75,7 +75,7 @@ int dampcf(int m, int damping, double *jac) {
  * @param[in] jac vector of damping coefficients
 **/
 
-double dif_eval(const double * v, int m, double thc, const double * jac){
+double dif_eval(const double *v, int m, double thc, const double *jac){
   double fval = 0.0;
   int j;
   for (j=0; j<=m; j++)
@@ -267,7 +267,7 @@ void chext(polparams *pol, double aIn, double bIn){
  * @param[in] size size of the arry
  * @return index of smallest entry in array
 **/
-int indexofSmallestElement(const double * array, int size){
+int indexofSmallestElement(const double *array, int size){
   int index = 0, i;
   for(i = 1; i < size; i++){
     if(array[i] < array[index])
@@ -288,7 +288,7 @@ int indexofSmallestElement(const double * array, int size){
  *
  * @warning mu must be preallocated
 **/
-int rootchb(int m, const double * v, const double * jac, double tha, double thb, double *mu,
+int rootchb(int m, const double *v, const double *jac, double tha, double thb, double *mu,
 	    double *thcOut){
   int MaxIterBalan = 30;     // max steps in Newton to balance interval
   double tolBal;
@@ -363,7 +363,7 @@ int rootchb(int m, const double * v, const double * jac, double tha, double thb,
  *  ChebAv]
  *
 **/
-int find_pol(const double * intv, polparams *pol) {
+int find_pol(const double *intv, polparams *pol) {
   double *mu, *v, *jac, t=0.0, itv[2],  vals[2];
   int max_deg=pol->max_deg;
   int min_deg=pol->min_deg;
diff --git a/EVSL_1.1.1/SRC/chebsi.c b/EVSL_1.1.1/SRC/chebsi.c
index 2e229a5..adbdf7e 100644
--- a/EVSL_1.1.1/SRC/chebsi.c
+++ b/EVSL_1.1.1/SRC/chebsi.c
@@ -36,8 +36,8 @@
  *   @warning Memory allocation for Yo/lamo/reso within this function
  */
 
-int ChebSI(int nev, const double * intv, int maxit,
-           double tol, const double * vinit, const polparams *pol, int *nevo,
+int ChebSI(int nev, const double *intv, int maxit,
+           double tol, const double *vinit, const polparams *pol, int *nevo,
            double **lamo, double **Yo, double **reso, FILE *fstats) {
   /*-------------------- for stats */
   double tm, tall=0.0, tmv=0.0;
diff --git a/EVSL_1.1.1/SRC/dos_utils.c b/EVSL_1.1.1/SRC/dos_utils.c
index 7900222..421283f 100644
--- a/EVSL_1.1.1/SRC/dos_utils.c
+++ b/EVSL_1.1.1/SRC/dos_utils.c
@@ -140,8 +140,8 @@ int apfun(double c, double h, const double *const xi,
  * @b Workspace
  * @param[in,out] w Work vector of length 3*n [allocate before call
  **/
-int pnav(const double * mu, int m, double cc, const double dd,
-         const double * v, double *y, double *w) {  // Really just ChebAv
+int pnav(const double *mu, int m, double cc, const double dd,
+         const double *v, double *y, double *w) {  // Really just ChebAv
   const int n = evsldata.n;
   /*-------------------- pointers to v_[k-1],v_[k], v_[k+1] from w */
   double *vk = w;
diff --git a/EVSL_1.1.1/SRC/dumps.c b/EVSL_1.1.1/SRC/dumps.c
index 348dd6c..3576820 100644
--- a/EVSL_1.1.1/SRC/dumps.c
+++ b/EVSL_1.1.1/SRC/dumps.c
@@ -21,8 +21,8 @@
  * @param[in] a Values
  * @param[in] fn filename
  */
-void save_mtx_basic(int nrow, int ncol, const int * ia,
-                    const int * ja, const double * a, const char *fn) {
+void save_mtx_basic(int nrow, int ncol, const int *ia,
+                    const int *ja, const double *a, const char *fn) {
   int i,j,nnz;
   FILE *fp = fopen(fn, "w");
 

From 4d30113a967e17feedbbc4c59dee1d70724261af Mon Sep 17 00:00:00 2001
From: Luke Erlandson <LukeAErlandson@gmail.com>
Date: Mon, 26 Nov 2018 22:39:21 -0500
Subject: [PATCH 4/5] Fix a const *

---
 EVSL_1.1.1/SRC/cheblanTr.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/EVSL_1.1.1/SRC/cheblanTr.c b/EVSL_1.1.1/SRC/cheblanTr.c
index bf10b24..5bb1daa 100644
--- a/EVSL_1.1.1/SRC/cheblanTr.c
+++ b/EVSL_1.1.1/SRC/cheblanTr.c
@@ -57,7 +57,7 @@
  *
  **/
 int ChebLanTr(int lanm, int nev, const double *intv, int maxit,
-              double tol, const double *vinit, const polparams* const pol, int *nev2,
+              double tol, const double *vinit, const polparams *pol, int *nev2,
               double **vals, double **W, double **resW, FILE *fstats) {
   const int ifGenEv = evsldata.ifGenEv;
   /*-------------------- for stats */

From 9f3203248c57912e744d45e459b4fc67b963d0fb Mon Sep 17 00:00:00 2001
From: Li <li50@roose.llnl.gov>
Date: Tue, 27 Nov 2018 16:10:24 -0800
Subject: [PATCH 5/5] some more const

---
 EVSL_1.1.1/DIRECT/Interface/evsl_cxsparse.c    | 6 +++---
 EVSL_1.1.1/DIRECT/Interface/evsl_direct.h      | 6 +++---
 EVSL_1.1.1/DIRECT/Interface/evsl_pardiso.c     | 6 +++---
 EVSL_1.1.1/DIRECT/Interface/evsl_suitesparse.c | 6 +++---
 EVSL_1.1.1/TESTS/Landos/LanDos.c               | 2 --
 EVSL_1.1.1/TESTS/PLanN/LapPLanN_MatFree.c      | 4 ++--
 EVSL_1.1.1/TESTS/RLanN/LapRLanN.c              | 1 -
 7 files changed, 14 insertions(+), 17 deletions(-)

diff --git a/EVSL_1.1.1/DIRECT/Interface/evsl_cxsparse.c b/EVSL_1.1.1/DIRECT/Interface/evsl_cxsparse.c
index eb5e34a..bbe05dd 100644
--- a/EVSL_1.1.1/DIRECT/Interface/evsl_cxsparse.c
+++ b/EVSL_1.1.1/DIRECT/Interface/evsl_cxsparse.c
@@ -129,7 +129,7 @@ int SetupBSolDirect(csrMat *B, void **data) {
 /** @brief Solver function of B
  *
  * */
-void BSolDirect(double *b, double *x, void *data) {
+void BSolDirect(const double *b, double *x, void *data) {
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
   cs_dis *S = Bsol_data->S;
   cs_din *N = Bsol_data->N;
@@ -145,7 +145,7 @@ void BSolDirect(double *b, double *x, void *data) {
 /** @brief Solver function of L^{T}
  *  x = L^{-T}*b
  * */
-void LTSolDirect(double *b, double *x, void *data) {
+void LTSolDirect(const double *b, double *x, void *data) {
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
   cs_dis *S = Bsol_data->S;
   cs_din *N = Bsol_data->N;
@@ -320,7 +320,7 @@ int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num,
  * @warning: This function MUST be of this prototype
  *
  *------------------------------------------------------------------*/
-void ASIGMABSolDirect(int n, double *br, double *bi, double *xr,
+void ASIGMABSolDirect(int n, const double *br, const double *bi, double *xr,
                       double *xz, void *data) {
 
   int i;
diff --git a/EVSL_1.1.1/DIRECT/Interface/evsl_direct.h b/EVSL_1.1.1/DIRECT/Interface/evsl_direct.h
index abf970c..4391131 100644
--- a/EVSL_1.1.1/DIRECT/Interface/evsl_direct.h
+++ b/EVSL_1.1.1/DIRECT/Interface/evsl_direct.h
@@ -15,13 +15,13 @@
 
 /* functions for B solve */
 int  SetupBSolDirect(csrMat *B, void **data);
-void BSolDirect(double *b, double *x, void *data);
-void LTSolDirect(double *b, double *x, void *data);
+void BSolDirect(const double *b, double *x, void *data);
+void LTSolDirect(const double *b, double *x, void *data);
 void FreeBSolDirectData(void *data);
 
 /* functions for A-SIGMA*B solve */
 int  SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num, EVSL_Complex *zk, void **data);
-void ASIGMABSolDirect(int n, double *br, double *bi, double *xr, double *xz, void *data);
+void ASIGMABSolDirect(int n, const double *br, const double *bi, double *xr, double *xz, void *data);
 void FreeASIGMABSolDirect(int num, void **data);
 
 /* evsl_direct_f90.c */
diff --git a/EVSL_1.1.1/DIRECT/Interface/evsl_pardiso.c b/EVSL_1.1.1/DIRECT/Interface/evsl_pardiso.c
index 7ca7032..0d4f3b6 100644
--- a/EVSL_1.1.1/DIRECT/Interface/evsl_pardiso.c
+++ b/EVSL_1.1.1/DIRECT/Interface/evsl_pardiso.c
@@ -189,7 +189,7 @@ int SetupBSolDirect(csrMat *B, void **data) {
 /** @brief Solver function of B
  *
  * */
-void BSolDirect(double *b, double *x, void *data) {
+void BSolDirect(const double *b, double *x, void *data) {
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
 
   MKL_INT maxfct = 1;
@@ -243,7 +243,7 @@ void BSolDirect(double *b, double *x, void *data) {
 /** @brief Solver function of L^{T}
  *  x = L^{-T}*b
  * */
-void LTSolDirect(double *b, double *x, void *data) {
+void LTSolDirect(const double *b, double *x, void *data) {
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
 
   MKL_INT maxfct = 1;
@@ -512,7 +512,7 @@ int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num,
  * @warning: This function MUST be of this prototype
  *
  *------------------------------------------------------------------*/
-void ASIGMABSolDirect(int n, double *br, double *bi, double *xr,
+void ASIGMABSolDirect(int n, const double *br, const double *bi, double *xr,
                       double *xz, void *data) {
 
   ASBSolDataDirect *sol_data = (ASBSolDataDirect *) data;
diff --git a/EVSL_1.1.1/DIRECT/Interface/evsl_suitesparse.c b/EVSL_1.1.1/DIRECT/Interface/evsl_suitesparse.c
index bdcc241..e1d9ff5 100644
--- a/EVSL_1.1.1/DIRECT/Interface/evsl_suitesparse.c
+++ b/EVSL_1.1.1/DIRECT/Interface/evsl_suitesparse.c
@@ -112,7 +112,7 @@ cholmod_dense cholmod_X, cholmod_B, *cholmod_Y=NULL, *cholmod_E=NULL,
 /** @brief Solver function of B
  *
  * */
-void BSolDirect(double *b, double *x, void *data) {
+void BSolDirect(const double *b, double *x, void *data) {
   //int n;
 
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
@@ -212,7 +212,7 @@ void FreeBSolDirectData(void *vdata) {
 /** @brief Solver function of L^{T}
  *  x = L^{-T}*b
  * */
-void LTSolDirect(double *b, double *x, void *data) {
+void LTSolDirect(const double *b, double *x, void *data) {
   //int n;
 
   BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
@@ -246,7 +246,7 @@ void LTSolDirect(double *b, double *x, void *data) {
  * @warning: This function MUST be of this prototype
  *
  *------------------------------------------------------------------*/
-void ASIGMABSolDirect(int n, double *br, double *bi, double *xr,
+void ASIGMABSolDirect(int n, const double *br, const double *bi, double *xr,
                       double *xz, void *data) {
   void* Numeric = data;
   double Control[UMFPACK_CONTROL];
diff --git a/EVSL_1.1.1/TESTS/Landos/LanDos.c b/EVSL_1.1.1/TESTS/Landos/LanDos.c
index 8a1753f..f1d00cd 100644
--- a/EVSL_1.1.1/TESTS/Landos/LanDos.c
+++ b/EVSL_1.1.1/TESTS/Landos/LanDos.c
@@ -13,8 +13,6 @@
 #include "io.h"
 #include "lapl.h"
 
-int exDOS(double *vals, int n, int npts, double *x, double *y, double *intv);
-
 /**
  *-----------------------------------------------------------------------
  * Tests landos.c -- Includes graphical comparison of exact DOS and calculated.
diff --git a/EVSL_1.1.1/TESTS/PLanN/LapPLanN_MatFree.c b/EVSL_1.1.1/TESTS/PLanN/LapPLanN_MatFree.c
index a5d9548..5f00d29 100644
--- a/EVSL_1.1.1/TESTS/PLanN/LapPLanN_MatFree.c
+++ b/EVSL_1.1.1/TESTS/PLanN/LapPLanN_MatFree.c
@@ -13,7 +13,7 @@
  * We provide a matvec routine, which only needs the stencil and grid sizes
  * The matvec routine and the associated data will need to be registered */
 /* matvec routine [it must be of this prototype] */
-void Lap2D3DMatvec(double *x, double *y, void *data);
+void Lap2D3DMatvec(const double *x, double *y, void *data);
 /* datatype for performing matvec for Laplacians */
 typedef struct _lapmv_t {
   int nx, ny, nz;
@@ -236,7 +236,7 @@ int main(int argc, char *argv[]) {
 }
 
 /*----------------- external matvec routine provided by users */
-void Lap2D3DMatvec(double *x, double *y, void *data) {
+void Lap2D3DMatvec(const double *x, double *y, void *data) {
   /* y = A * x
    * data: pointer to a struct that contains all needed data
    */
diff --git a/EVSL_1.1.1/TESTS/RLanN/LapRLanN.c b/EVSL_1.1.1/TESTS/RLanN/LapRLanN.c
index 252167b..906bcff 100644
--- a/EVSL_1.1.1/TESTS/RLanN/LapRLanN.c
+++ b/EVSL_1.1.1/TESTS/RLanN/LapRLanN.c
@@ -3,7 +3,6 @@
 #include <string.h>
 #include <sys/stat.h>
 #include <math.h>
-#include <complex.h>
 #include "evsl.h"
 #include "io.h"
 #include "evsl_direct.h"