improved lanczos

lanczos.h

@@ -50,7 +50,7 @@ if ( n!= (int)bigmat.ncols()) laerror("non-square matrix in lanczos");
 if(eivals.size()<nroots) laerror("too small eivals dimension in lanczos");
 
 NRVec<T> vec1(n),vec2(n);
-NRVec<T> *v0,*v1;
+NRVec<T> *v0;
 AuxStorage<T> *s0,*s1;
 
 if(incore)
@@ -112,9 +112,17 @@ for(j=1; j<maxkrylov;++j)
 		}
 	else
 		{
-		laerror("zero norm in lanczos");
+		//generate an arbitrary vector and orthonormalize it to all previous v_j
-		//could generate an arbitrary vector and orthonormalize it
+		vec2.randomize(1.);
+		NRVec<T> vec3(n);
+		for(int k=0; k<j; ++k)
+			{
+			T f = vec3.dot(vec2);
+			vec2.axpy(-f,vec3);
+			}
+		vec2.normalize();
 		}
+	//vec2 now stores v_j
 	if(incore) v0[j]=vec2; else s0->put(vec2,j);
 	vec1 *= -beta[j-1];
 	bigmat.gemv(1,vec1,'n',1,vec2);

t.cc

@@ -4650,7 +4650,7 @@ if(n<=20)
 	}
 }
 
-if(1);
+if(0)
 {
 int n,m;
 cin>> n>>m;
@@ -4660,4 +4660,36 @@ NRMat<double> b = explicit_matrix<double,NRMat<double> >(a);
 cout <<"Error = "<<(a-b).norm()<<endl;
 }
 
+if(1)
+{
+int m;
+int which;
+cin>>m >>which;
+cout <<"here\n";
+NRSMat<double> a;
+cin>>a;
+int n=a.nrows();
+NRVec<double> rr(n);
+
+NRSMat<double> aa(a);
+NRMat<double> vv(n,n);
+diagonalize(aa,rr,&vv);
+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);
+cout <<"Iterative energies " <<r;
+cout <<"Eigenvectors compare:\n";
+for(int i=0; i<m; ++i)
+        {
+        cout <<eivecs[i];
+        for(int j=0; j<n;++j) cout <<vv[j][i]<<" ";
+        cout <<endl;
+        }
+
+}
+
+
 }//main