From 36983222e803c61c56f49d943a76ccb8dc8eddd0 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Tue, 10 Feb 2026 15:18:56 +0100
Subject: [PATCH] improevd lanczos

---
 lanczos.h | 16 +++++++++++++---
 t.cc      |  9 +++++++--
 2 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/lanczos.h b/lanczos.h
index d7ab205..b4b756c 100644
--- a/lanczos.h
+++ b/lanczos.h
@@ -114,11 +114,21 @@ for(j=1; j<maxkrylov;++j)
 		{
 		//generate an arbitrary vector and orthonormalize it to all previous v_j
 		vec2.randomize(1.);
-		NRVec<T> vec3(n);
 		for(int k=0; k<j; ++k)
 			{
-			T f = vec3.dot(vec2);
-			vec2.axpy(-f,vec3);
+			T f;
+			if(incore) 
+				{
+				f = v0[k].dot(vec2);
+				vec2.axpy(-f,v0[k]);
+				}
+			else
+				{
+				NRVec<T> vec3(n);
+				s0->get(vec3,k);
+                                f = vec3.dot(vec2);
+                                vec2.axpy(-f,vec3);
+				}
 			}
 		vec2.normalize();
 		}
diff --git a/t.cc b/t.cc
index d023ad5..c835e60 100644
--- a/t.cc
+++ b/t.cc
@@ -61,6 +61,11 @@ inline int randind(const int n)
 
 complex<double> mycident (const complex<double>&x) {return x;}
 
+void randomguess(NRVec<double> &v)
+{
+v.randomize(1);
+}
+
 void printme(const NRPerm<int> &p)
 {
 PERM_RANK_TYPE rank=p.rank();
@@ -4678,8 +4683,8 @@ NRVec<double> r(m);
 NRVec<double> *eivecs = new NRVec<double>[m];
 cout <<"Exact energies " <<rr<<endl;
 
-if(which) lanczos(a,r,eivecs,NULL,m,1,1e-6,1,200,300);
-else davidson(a,r,eivecs,NULL,m,1,1e-6,1,200,300);
+if(which) lanczos(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
+else davidson(a,r,eivecs,NULL,m,1,1e-6,1,500,500,randomguess);
 cout <<"Iterative energies " <<r;
 cout <<"Eigenvectors compare:\n";
 for(int i=0; i<m; ++i)