From 1cb536421f7451248070c91eb0754c5b13a56ea6 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 9 Dec 2025 16:03:04 +0100
Subject: [PATCH] davidson: added root targeting

---
 davidson.h | 22 ++++++++++++++++++++--
 t.cc       | 27 ++++++++++++++++++++++++++-
 2 files changed, 46 insertions(+), 3 deletions(-)

diff --git a/davidson.h b/davidson.h
index 21669a3..2dde47c 100644
--- a/davidson.h
+++ b/davidson.h
@@ -33,7 +33,7 @@ namespace LA {
 //Note that for efficiency in a direct CI case the diagonalof() should cache its result
 
 
-//works for complex hermitian-only too, but does not converge sowell for more than 1 root
+//works for complex hermitian-only too, but does not converge so well for more than 1 root
 //@@@ for large krylov spaces >200 it can occur 'convergence problem in sygv/syev in diagonalize()'
 //@@@options: left eigenvectors by matrix transpose, overridesymmetric (for nrmat)
 //@@@small matrix gdiagonalize - shift complex roots up (option to gdiagonalize?)
@@ -48,7 +48,7 @@ template <typename T, typename Matrix>
 extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
 		int nroots=1,  const bool verbose=0, const double eps=1e-6,
 	 	const bool incore=1, int maxit=100, const int maxkrylov = 500,
-		void (*initguess)(NRVec<T> &)=NULL)
+		void (*initguess)(NRVec<T> &)=NULL, const typename LA_traits<T>::normtype *target=NULL)
 {
 bool flag=0;
 int n=bigmat.nrows();
@@ -155,6 +155,24 @@ else
 	for(int i=0; i<=krylovsize; ++i) {r[i]/=LA_traits<T>::realpart(beta[i]); ri[i]/=LA_traits<T>::realpart(beta[i]);}
 	}
 
+if(target) //resort eigenpairs by distance from the target
+	{
+	NRVec<typename LA_traits<T>::normtype> key(krylovsize+1);
+	for(int i=0; i<=krylovsize; ++i) key[i] = abs(r[i] - *target);
+	NRPerm<int> p(krylovsize+1);
+	key.sort(0,p);
+	
+	NRVec<typename LA_traits<T>::normtype> rp(maxkrylov);
+	NRMat<T> smallVp(maxkrylov,maxkrylov);
+	for(int i=0; i<=krylovsize; ++i)
+		{
+		rp[i]= r[p[i+1]-1];
+		for(int j=0; j<=krylovsize; ++j) smallVp(j,i) = smallV(j,p[i+1]-1);
+		}
+	r = rp;
+	smallV = smallVp;
+	}
+
 typename LA_traits<T>::normtype eival_n=r[nroot];
 oldnroot=nroot;
 typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));
diff --git a/t.cc b/t.cc
index 5e65b65..b68640e 100644
--- a/t.cc
+++ b/t.cc
@@ -4500,7 +4500,7 @@ cout <<"Error = "<<(xt-y).norm()<<endl;
 
 }
 
-if(1)
+if(0)
 {
 //n must not be too small
 int n,m;
@@ -4537,4 +4537,29 @@ for(int i=0; i<m; ++i)
 
 }
 
+if(1)
+{
+int n,m;
+cin>>n >>m;
+NRSMat<double> a(n);
+NRVec<double> rr(n);
+
+for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
+        {
+        a(i,j)= RANDDOUBLE();
+	if(i==j) a(i,i)+= .5*(i-n*.5);
+        }
+
+NRSMat<double> aa;
+NRMat<double> vv(n,n);
+aa=a; diagonalize(aa,rr,&vv);
+NRVec<double> r(m);
+NRVec<double> *eivecs = new NRVec<double>[m];
+double target= 0;
+cout <<"Exact energies " <<rr<<endl;
+
+davidson(a,r,eivecs,NULL,m,1,1e-6,1,200,300,(void (*)(LA::NRVec<double>&))NULL,&target);
+cout <<"Davidson energies " <<r;
+}
+
 }//main