davidson for complex hermitean

davidson.h

@@ -33,7 +33,7 @@ namespace LA {
 //Note that for efficiency in a direct CI case the diagonalof() should cache its result
 
 
-//@@@should work for complex hermitian-only too, but was not tested yet (maybe somwhere complex conjugation will have to be added)
+//works for complex hermitian-only too, but does not converge sowell 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?)
@@ -158,12 +158,19 @@ else
 typename LA_traits<T>::normtype eival_n=r[nroot];
 oldnroot=nroot;
 typename LA_traits<T>::normtype test=std::abs(smallV(krylovsize,nroot));
+//std::cout <<"NROOT = "<<nroot<<" TEST = "<<test<<std::endl;
+if(test>100.)
+	{
+	std::cout <<"Divergence in Davidson\n";
+	flag=1;
+	goto finished;
+	}
 if(test<eps) nroot++;
 if(it==0) nroot = 0;
 for(int iroot=0; iroot<=std::min(krylovsize,nroots-1); ++iroot)
 	{
         test = std::abs(smallV(krylovsize,iroot));
-        if(test>eps) nroot=std::min(nroot,iroot);
+        if(test>eps) nroot=std::min(nroot,iroot); //return to previous root
         if(verbose && iroot<=std::max(oldnroot,nroot)) 
 		{
 		std::cout <<"Davidson: iter="<<it <<" dim="<<krylovsize<<" root="<<iroot<<" eigenvalue="<<r[iroot]<<"\n";
@@ -229,7 +236,8 @@ for(j=0; j<=krylovsize; ++j)
 	typename LA_traits<T>::normtype vnorm2;
         if(!incore) s0->get(vec2,j);
 	do	{
-          	T ab = vec1*(incore?v0[j]:vec2) /smallS(j,j);
+		 // For COMPLEX the order matters, the first one gets conjuagted, changed vec1 to the right!
+          	T ab = (incore?v0[j]:vec2)*vec1 /smallS(j,j);
             	vec1.axpy(-ab,incore?v0[j]:vec2);
 		vnorm2 = vec1.norm();
 		if(vnorm2==0) 

t.cc

@@ -1137,7 +1137,7 @@ if(0)
 {
 int n,m;
 cin>>n >>m;
-NRSMat<double> a(n,n);
+NRSMat<double> a(n);
 NRVec<double> rr(n);
 
 for(int i=0;i<n;++i) for(int j=0;j<=i;++j)
@@ -4460,7 +4460,7 @@ if(inv==false)
 }
 
 
-if(1)
+if(0)
 {
 //test linear transform
 int r=3;
@@ -4500,4 +4500,41 @@ cout <<"Error = "<<(xt-y).norm()<<endl;
 
 }
 
+if(1)
+{
+//n must not be too small
+int n,m;
+cin>>n >>m;
+NRSMat<complex<double> > a(n);
+a.randomize(.5);
+
+for(int i=0;i<n;++i)
+        {
+	a(i,i)=  RANDDOUBLE() +  .2*(i-n);
+        }
+//cout <<"test matrix = "<<a;
+
+NRVec<double> rr(n);
+NRSMat<complex<double> > aa;
+NRMat<complex<double> > vv(n,n);
+aa=a; diagonalize(aa,rr,&vv);
+cout <<"Exact energies "<<rr;
+
+NRVec<complex<double> > r(m);
+NRVec<complex<double> > *eivecs = new NRVec<complex<double> >[m];
+davidson(a,r,eivecs,NULL,m,true,1e-6,true,10*n,n*10);
+
+cout <<"Davidson energies " <<r;
+cout <<"Exact energies "<<rr;
+
+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

vec.h

@@ -311,10 +311,10 @@ public:
 
 	//! compute the Euclidean inner product (with conjugation in complex case) 
 	inline const T operator*(const NRVec &rhs) const;
-	inline const T dot(const NRVec &rhs) const {return *this * rhs;};
 
 	//! compute the Euclidean inner product (with conjugation in complex case) with a stride-vector
-	inline const T dot(const T *a, const int stride = 1) const;
+	inline const T dot(const T *a, const int stride = 1, bool conjugate=true) const;
+	inline const T dot(const NRVec &rhs, bool conjugate=true) const {return dot(&((*this)[0]),1,conjugate);};
 
 	void gemv(const T beta, const NRMat<T> &a, const char trans, const T alpha, const NRVec &x);
 	void gemv(const T beta, const NRSMat<T> &a, const char trans /**< just for compatibility reasons */, const T alpha, const NRVec &x);
@@ -1073,10 +1073,25 @@ inline const T NRVec<T>::operator*(const NRVec<T> &rhs) const {
 	NOT_GPU(rhs);
 
 	T dot(0);
-	for(register int i=0; i<nn; ++i) dot += v[i]*rhs.v[i];
+	if(LA_traits<T>::is_complex()) for(register int i=0; i<nn; ++i) dot += LA_traits<T>::conjugate(v[i])*rhs.v[i];
+	else for(register int i=0; i<nn; ++i) dot += v[i]*rhs.v[i];
 	return dot;
 }
 
+/***************************************************************************
+ * compute dot product with optional conjugation
+ */
+template<typename T>
+inline const T NRVec<T>::dot(const T *a, const int stride , bool conjugate) const
+{
+	 NOT_GPU(*this);
+	T dot(0);
+        if(LA_traits<T>::is_complex() && conjugate) for(register int i=0; i<nn; ++i) dot += LA_traits<T>::conjugate(v[i])*a[i*stride];
+        else for(register int i=0; i<nn; ++i) dot += v[i]*a[i*stride];
+        return dot; 
+}
+
+
 /***************************************************************************//**
  * indexing operator giving the element at given position with range checking in
  * the DEBUG mode
@@ -1801,9 +1816,9 @@ inline const std::complex<double> NRVec<std::complex<double> >::operator*(const
  * @return \f$\sum_{i=1}^N\vec{x}_{i}\cdot y_{\mathrm{stride}\cdot(i-1) + 1}\f$
  ******************************************************************************/
 template<>
-inline const double NRVec<double>::dot(const double *y, const int stride) const {
+inline const double NRVec<double>::dot(const double *y, const int stride, bool conjugate) const {
 	NOT_GPU(*this);
-	return cblas_ddot(nn, y, stride, v, 1);
+	return cblas_ddot(nn,v, 1,y,stride);
 }
 
 /***************************************************************************//**
@@ -1813,10 +1828,11 @@ inline const double NRVec<double>::dot(const double *y, const int stride) const
  * @return \f$\sum_{i=1}^N\vec{x}_{i}\cdot \overbar{y_{\mathrm{stride}\cdot(i-1) + 1}}\f$
  ******************************************************************************/
 template<>
-inline const std::complex<double> NRVec<std::complex<double> >::dot(const std::complex<double> *y, const int stride) const {
+inline const std::complex<double> NRVec<std::complex<double> >::dot(const std::complex<double> *y, const int stride, bool conjugate) const {
         std::complex<double> dot;
 	NOT_GPU(*this);
-        cblas_zdotc_sub(nn, y, stride, v, 1, &dot);
+        if(conjugate) cblas_zdotc_sub(nn, v,1,y, stride, &dot);
+	else cblas_zdotc_sub(nn, v,1,y, stride, &dot);
         return dot;
 }