*** empty log message ***

davidson.h

@@ -1,12 +1,15 @@
+#ifndef _davidson_h
+#define _davidson_h
 #include "vec.h"
 #include "smat.h"
 #include "mat.h"
 #include "sparsemat.h"
 #include "nonclass.h"
+#include "auxstorage.h"
 
 //Davidson diagonalization of real symmetric matrix (modified Lanczos)
 
-//matrix can be any class which has nrows(), ncols(), diagonalof() and NRVec::gemv() available
+//matrix can be any class which has nrows(), ncols(), diagonalof(), issymmetric(), and gemv() available
 //does not even have to be explicitly stored - direct CI
 
 export template <typename T, typename Matrix>
@@ -17,3 +20,220 @@ extern void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, c
 //@@@options: left eigenvectors by matrix transpose, overridesymmetric (for nrmat)
 //@@@small matrix gdiagonalize - shift complex roots up (option to gdiagonalize?)
 //@@@test gdiagonalize whether it sorts the roots and what for complex ones
+
+
+//Davidson algorithm: J. Comp. Phys. 17:817 (1975) 
+
+//@@@implement left eigenvectors for nonsymmetric case
+
+template <typename T, typename Matrix>
+void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile,
+                int nroots,  const bool verbose, const double eps, 
+                const bool incore, int maxit, const int maxkrylov)
+{
+bool flag=0;
+int n=bigmat.nrows();
+if ( n!= (int)bigmat.ncols()) laerror("non-square matrix in davidson");
+if(eivals.size()<nroots) laerror("too small eivals dimension in davidson");
+
+NRVec<T> vec1(n),vec2(n);
+NRMat<T> smallH(maxkrylov,maxkrylov),smallS(maxkrylov,maxkrylov),smallV;
+NRVec<T> r(maxkrylov);
+NRVec<T> *v0,*v1;
+AuxStorage<T> *s0,*s1;
+
+if(incore)
+	{
+	v0 = new NRVec<T>[maxkrylov];
+	v1 = new NRVec<T>[maxkrylov];
+	}
+else
+	{
+	s0 = new AuxStorage<T>;
+	s1 = new AuxStorage<T>;
+	}
+
+int i,j;
+
+if(maxkrylov<maxit) maxit=maxkrylov;
+if(nroots>=maxkrylov) nroots =maxkrylov-1;
+int  nroot=0;
+int oldnroot;
+smallS=0;
+smallH=0;     
+//guess based on lowest diagonal element of the matrix
+bigmat.diagonalof(vec2);
+vec1=0;
+{T t=1e100; int i,j;  
+for(i=0, j= -1; i<n; ++i) if(vec2[i]<t) {t=vec2[i]; j=i;}
+vec1[j]=1;}
+
+//init Krylov matrices
+bigmat.gemv(0,vec2,'n',1,vec1); //avoid bigmat.operator*(vec), since that needs to allocate another n-sized vector
+smallH(0,0) = vec1*vec2;
+smallS(0,0) = vec1*vec1;
+int krylovsize = 0;
+if(incore) v0[0]=vec1; else s0->put(vec1,0); 
+if(incore) v1[0]=vec2; else s1->put(vec2,0);
+
+
+//iterative Davidson loop
+int it;
+for(it=0; it<maxit; ++it)
+{
+if(it>0) //if this is the first iteration just need to diagonalise the matrix
+	{
+	//update reduced overlap matrix
+	if(incore) v0[krylovsize]=vec1; else s0->put(vec1,krylovsize);
+	for(j=0; j<krylovsize; ++j)
+		{
+		if(!incore) s0->get(vec2,j);
+            	smallS(krylovsize,j) = smallS(j,krylovsize) = vec1*(incore?v0[j]:vec2);
+		}
+        smallS(krylovsize,krylovsize) = vec1*vec1;
+	bigmat.gemv(0,vec2,'n',1,vec1);
+	if(incore) v1[krylovsize]=vec2; else s1->put(vec2,krylovsize);
+
+	//update reduced hamiltonian matrix
+        smallH(krylovsize,krylovsize) = vec1*vec2;
+	for(j=0; j<krylovsize; ++j)
+		{
+		if(!incore) s0->get(vec1,j);
+		smallH(j,krylovsize) = (incore?v0[j]:vec1)*vec2;
+            	if(bigmat.issymmetric()) smallH(krylovsize,j) = smallH(j,krylovsize);
+		}
+	if(!bigmat.issymmetric())
+		{
+		if(!incore) s0->get(vec1,krylovsize);
+		for(j=0; j<krylovsize; ++j)
+			{
+			if(!incore) s1->get(vec2,j);
+			smallH(krylovsize,j) = incore? v1[j]*v0[krylovsize] :vec1*vec2;
+			}
+		}
+	}
+smallV=smallH;
+NRMat<T> smallSwork=smallS;
+if(bigmat.issymmetric()) 
+	diagonalize(smallV,r,1,1,krylovsize+1,&smallSwork,1); //for symmetric matrix they have already been sorted to ascending order in lapack
+else 
+	{
+	NRVec<T> ri(krylovsize+1),beta(krylovsize+1);
+	NRMat<T> scratch;
+	scratch=smallV;
+	gdiagonalize(scratch, r, ri,NULL, &smallV, 1, krylovsize+1, 2, 0, &smallSwork, &beta);
+	for(int i=0; i<=krylovsize; ++i) {r[i]/=beta[i]; ri[i]/=beta[i];}
+	}
+
+T eival_n=r[nroot];
+oldnroot=nroot;
+typename LA_traits<T>::normtype test=abs(smallV(krylovsize,nroot));
+if(test<eps) nroot++;
+if(it==0) nroot = 0;
+for(int iroot=0; iroot<=min(krylovsize,nroots-1); ++iroot)
+	{
+        test = abs(smallV(krylovsize,iroot));
+        if(test>eps) nroot=min(nroot,iroot);
+        if(verbose && iroot<=max(oldnroot,nroot)) 
+		{
+		cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" energy="<<r[iroot]<<"\n";
+		}
+	}
+
+if(verbose && oldnroot!=nroot) cout <<"root no. "<<oldnroot<<" converged\n";
+if (nroot>=nroots) goto converged;
+if (it==maxit-1) break; //not converged
+
+if (krylovsize==maxkrylov) //restart, krylov space exceeded
+	{
+	if(nroot!=0) {flag=1; goto finished;}
+	smallH=0;
+	smallS=0;
+	vec1=0;
+	for(i=0; i<=krylovsize; ++i)
+		{
+		if(!incore) s0->get(vec2,i);
+                vec1.axpy(smallV(i,0),incore?v0[i]:vec2);
+		}
+	s0->put(vec1,0);
+	vec1.normalize();
+        krylovsize = 0;
+	continue;
+	}
+
+//generate the update vector
+vec1=0;
+for(j=0; j<=krylovsize; ++j)
+	{
+	if(!incore) s0->get(vec2,j);
+        vec1.axpy(-r[nroot]*smallV(j,nroot),incore?v0[j]:vec2);
+	if(!incore) s1->get(vec2,j);
+        vec1.axpy(smallV(j,nroot),incore?v1[j]:vec2);
+	}
+bigmat.diagonalof(vec2);
+eival_n = r[nroot];
+for(i=0; i<n; ++i)
+	{
+	T denom = vec2[i] - eival_n;
+        denom = denom<0?-max(0.1,abs(denom)):max(0.1,abs(denom));
+        vec1[i] /= denom;
+	}
+
+//orthogonalise to previous vectors
+vec1.normalize();
+for(j=0; j<=krylovsize; ++j)
+	{
+	typename LA_traits<T>::normtype vnorm;
+        if(!incore) s0->get(vec2,j);
+	do	{
+          	T ab = vec1*(incore?v0[j]:vec2) /smallS(j,j);
+            	vec1.axpy(-ab,incore?v0[j]:vec2);
+		vnorm = vec1.norm();
+		vec1 *= (1./vnorm);
+		} while (vnorm<0.99);
+	}
+
+//here it is possible to apply some purification procedure if the eivector has to fulfill other conditions
+//vec1.normalize(); //after the purification
+
+++krylovsize; //enlarge Krylov space
+}
+flag=1;
+goto finished;
+
+converged:
+AuxStorage<typename LA_traits<T>::elementtype> *ev;
+if(eivecsfile) ev = new AuxStorage<typename LA_traits<T>::elementtype>(eivecsfile);
+if(verbose) cout << "Davidson converged in "<<it<<" iterations.\n";
+for(nroot=0; nroot<nroots; ++nroot)
+	{
+        eivals[nroot]=r[nroot];
+	if(eivecs)
+		{
+		vec1=0;
+		for(j=0; j<=krylovsize; ++j )
+			{
+            		if(!incore) s0->get(vec2,j);
+			vec1.axpy(smallV(j,nroot),incore?v0[j]:vec2);
+			}
+		vec1.normalize();
+        	if(eivecs) eivecs[nroot]|=vec1;
+		if(eivecsfile)
+			{
+			ev->put(vec1,nroot);
+			}
+		}
+	}
+
+if(eivecsfile) delete ev;
+
+finished:
+if(incore) {delete[] v0; delete[] v1;}
+else  {delete s0; delete s1;}
+
+if(flag) laerror("no convergence in davidson");
+}
+
+
+
+#endif