diff --git a/davidson.h b/davidson.h
index 5e46254..5840b3f 100644
--- a/davidson.h
+++ b/davidson.h
@@ -4,11 +4,12 @@
 #include "sparsemat.h"
 #include "nonclass.h"
 
-//Davidson diagonalization of real symmetric matrix
+//Davidson diagonalization of real symmetric matrix (modified Lanczos)
 
 //matrix can be any class which has nrows(), ncols(), diagonalof() and NRVec::gemv() available
 //does not even have to be explicitly stored - direct CI
-//n<0 highest eigenvalues, n>0 lowest eigenvalues
 
 export template <typename T, typename Matrix>
-extern void davidson(const Matrix &bigmat, NRVec<T> *eivecs /*input-output*/, NRVec<T> &eivals, int nroots=1, const double accur=1e-6, int nits=50, const int ndvdmx = 500, const bool incore=1);
+extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, void *eivecs, 
+		int nroots=1,  const bool verbose=0, const double eps=1e-6,
+	 	const bool incore=1, int maxit=100, const int maxkrylov = 500);
diff --git a/nonclass.cc b/nonclass.cc
index ef87064..a5e5e60 100644
--- a/nonclass.cc
+++ b/nonclass.cc
@@ -111,6 +111,7 @@ void laerror(const char *s1, const char *s2, const char *s3, const char *s4)
     if(s4) std::cerr << s4;
   }
   std::cerr << endl;
+//@@@throw
   exit(1);
 }
 
@@ -119,17 +120,14 @@ void laerror(const char *s1, const char *s2, const char *s3, const char *s4)
 //////////////////////
 
 // A will be overwritten, B will contain the solutions, A is nxn, B is rhs x n
-void linear_solve(NRMat<double> &A, NRMat<double> *B, double *det)
+static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const int ldb, double *det)
 {
 	int r, *ipiv;
 	
 	if (A.nrows() != A.ncols()) laerror("linear_solve() call for non-square matrix");
-	if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()");
 	A.copyonwrite();
-	if (B) B->copyonwrite();
 	ipiv = new int[A.nrows()];
-	r = clapack_dgesv(CblasRowMajor, A.nrows(), B ? B->nrows() : 0, A[0], A.ncols(),
-			ipiv, B ? B[0] : (double *)0, B ? B->ncols() : A.nrows());
+	r = clapack_dgesv(CblasRowMajor, A.nrows(), nrhs, A[0], A.ncols(), ipiv, B , ldb);
 	if (r < 0) {
 		delete[] ipiv;
 		laerror("illegal argument in lapack_gesv");
@@ -144,28 +142,36 @@ void linear_solve(NRMat<double> &A, NRMat<double> *B, double *det)
 	if (r>0 && B) laerror("singular matrix in lapack_gesv");
 }
 
+void linear_solve(NRMat<double> &A, NRMat<double> *B, double *det)
+{
+if (B && A.nrows() != B->ncols()) laerror("incompatible matrices in linear_solve()");
+if(B) B->copyonwrite();
+linear_solve_do(A,B?(*B)[0]:NULL,B?B->nrows() : 0,  B?B->ncols():A.nrows(), det);
+}
 
-// Next routines are not available in clapack, fotran ones will b used with an
+void linear_solve(NRMat<double> &A, NRVec<double> &B, double *det)
+{
+if(A.nrows() != B.size()) laerror("incompatible matrices in linear_solve()");
+B.copyonwrite();
+linear_solve_do(A,&B[0],1,A.nrows(),det);
+}
+
+
+// Next routines are not available in clapack, fotran ones will be used with an
 // additional swap/transpose of outputs when needed
 
 extern "C" void FORNAME(dspsv)(const char *UPLO, const int *N, const int *NRHS,
 		double *AP, int *IPIV, double *B, const int *LDB, int *INFO);
 
-void linear_solve(NRSMat<double> &a, NRMat<double> *b, double *det)
+static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const int ldb, double *det)
 {
 	int r, *ipiv;
 	if (det) cerr << "@@@ sign of the determinant not implemented correctly yet\n";
-	if (b && a.nrows() != b->ncols())
-		laerror("incompatible matrices in symmetric linear_solve()");
 	a.copyonwrite();
-	if (b) b->copyonwrite();
 	ipiv = new int[a.nrows()];
 	char U = 'U';
 	int n = a.nrows();
-	int nrhs = 0;
-	if (b) nrhs = b->nrows();
-	int ldb = b ? b->ncols() : a.nrows();
-	FORNAME(dspsv)(&U, &n, &nrhs, a, ipiv, b?(*b)[0]:0, &ldb,&r);
+	FORNAME(dspsv)(&U, &n, &nrhs, a, ipiv, b, &ldb,&r);
 	if (r < 0) {
 		delete[] ipiv;
 		laerror("illegal argument in spsv() call of linear_solve()");
@@ -181,6 +187,23 @@ void linear_solve(NRSMat<double> &a, NRMat<double> *b, double *det)
 }
 
 
+void linear_solve(NRSMat<double> &a, NRMat<double> *B, double *det)
+{
+	if (B && a.nrows() != B->ncols())
+		laerror("incompatible matrices in symmetric linear_solve()");
+	if (B) B->copyonwrite();
+linear_solve_do(a,B?(*B)[0]:NULL,B?B->nrows() : 0, B?B->ncols():a.nrows(),det);
+}
+
+
+void linear_solve(NRSMat<double> &a, NRVec<double> &B, double *det)
+{
+	if (a.nrows() != B.size())
+		laerror("incompatible matrices in symmetric linear_solve()");
+	B.copyonwrite();
+linear_solve_do(a,&B[0],1,a.nrows(),det);
+}
+
 extern "C" void FORNAME(dsyev)(const char *JOBZ, const char *UPLO, const int *N,
 		double *A, const int *LDA, double *W, double *WORK, const int *LWORK, int *INFO);
 
@@ -547,56 +570,8 @@ if(n==0) n=m;
 if(n<0 || n>m) laerror("actual dimension in gendiagonalize out of range");
 
 //transform the problem to usual diagonalization
-/*
-c
-c     this routine is a translation of the algol procedure reduc1,
-c     num. math. 11, 99-110(1968) by martin and wilkinson.
-c     handbook for auto. comp., vol.ii-linear algebra, 303-314(1971).
-c
-c     this routine reduces the generalized symmetric eigenproblem
-c     ax=(lambda)bx, where b is positive definite, to the standard
-c     symmetric eigenproblem using the cholesky factorization of b.
-c
-c     on input
-c
-c        nm must be set to the row dimension of two-dimensional
-c          array parameters as declared in the calling program
-c          dimension statement.
-c
-c        n is the order of the matrices a and b.  if the cholesky
-c          factor l of b is already available, n should be prefixed
-c          with a minus sign.
-c
-c        a and b contain the real symmetric input matrices.  only the
-c          full upper triangles of the matrices need be supplied.  if
-c          n is negative, the strict lower triangle of b contains,
-c          instead, the strict lower triangle of its cholesky factor l.
-c
-c        dl contains, if n is negative, the diagonal elements of l.
-c
-c     on output
-c
-c        a contains in its full lower triangle the full lower triangle
-c          of the symmetric matrix derived from the reduction to the
-c          standard form.  the strict upper triangle of a is unaltered.
-c
-c        b contains in its strict lower triangle the strict lower
-c          triangle of its cholesky factor l.  the full upper
-c          triangle of b is unaltered.
-c
-c        dl contains the diagonal elements of l.
-c
-c        ierr is set to
-c          zero       for normal return,
-c          7*n+1      if b is not positive definite.
-c
-c     questions and comments should be directed to burton s. garbow,
-c     mathematics and computer science div, argonne national laboratory
-c
-c     this version dated august 1983.
-c
-*/
-//    .......... form l in the arrays b and dl ..........
+
+//cholesky decompose in b and dl
 for(i=0; i<n; ++i)
 {
 for(j=i; j<n; ++j)
@@ -612,8 +587,7 @@ for(j=i; j<n; ++j)
 	}
 }
 
-//     .......... form the transpose of the upper triangle of inv(l)*a
-//                in the lower triangle of the array a ..........
+// form the transpose of the upper triangle of inv(l)*a in the lower triangle of a
 for(i=0; i<n; ++i)
 {
 for(j=i; j<n ; ++j)
@@ -623,7 +597,7 @@ for(j=i; j<n ; ++j)
 	}
 }
 
-//     .......... pre-multiply by inv(l) and overwrite ..........
+//pre-multiply by l^-1
 for(j=0; j<n ; ++j)
 {
 	for(i=j;i<n;++i)
@@ -641,38 +615,6 @@ for(i=1;i<n;++i) for(j=0; j<i; ++j) a(j,i)=a(i,j);
 diagonalize(a,w,1,1,n);
 
 //transform the eigenvectors back
-/*
-c
-c     this routine is a translation of the algol procedure rebaka,
-c     num. math. 11, 99-110(1968) by martin and wilkinson.
-c     handbook for auto. comp., vol.ii-linear algebra, 303-314(1971).
-c
-c     this routine forms the eigenvectors of a generalized
-c     symmetric eigensystem by back transforming those of the
-c     derived symmetric matrix determined by  reduc.
-c
-c     on input
-c
-c        n is the order of the matrix system.
-c
-c        b contains information about the similarity transformation
-c          (cholesky decomposition) used in the reduction by  reduc
-c          in its strict lower triangle.
-c
-c        dl contains further information about the transformation.
-c
-c        a contains the eigenvectors to be back transformed
-c          in its  columns
-c
-c     on output
-c
-c        a contains the transformed eigenvectors
-c          in its columns
-c
-c     questions and comments should be directed to burton s. garbow,
-c     mathematics and computer science div, argonne national laboratory
-c
-*/
 for(j=0; j<n; ++j)//eigenvector loop
 {
 for(int i=n-1; i>=0; --i)//component loop