From ea56e7380d4457af0f257f530cf219778bc5e829 Mon Sep 17 00:00:00 2001
From: jiri <jiri>
Date: Thu, 17 Feb 2005 00:30:09 +0000
Subject: [PATCH] *** empty log message ***

---
 diis.h | 59 +++++++++++++++++++++++++++++++++++++---------------------
 1 file changed, 38 insertions(+), 21 deletions(-)

diff --git a/diis.h b/diis.h
index 668a37c..6feb1bf 100644
--- a/diis.h
+++ b/diis.h
@@ -25,15 +25,15 @@ class DIIS
 public:
 	DIIS(const int n, const bool core=1);
 	~DIIS();
-	Te extrapolate(T &vec); //vec is input/output; returns square residual norm
+	typename LA_traits<T>::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm
 	};
 
 template<typename T>
-DIIS<T>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n+1,n+1)
+DIIS<T>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n)
 {
 st=incore?NULL: new AuxStorage<Te>;
 stor= incore? new T[dim] : NULL;
-bmat= (Te)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Te)-1; 
+bmat= (Te)0; for(int i=1; i<n; ++i) bmat(0,i) = (Te)-1; 
 aktdim=cyclicshift=0;
 
 }
@@ -47,13 +47,13 @@ if(stor) delete[] stor;
 
 
 template<typename T>
-typename DIIS<T>::Te DIIS<T>::extrapolate(T &vec)
+typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec)
 {
 //if dim exceeded, shift 
 if(aktdim==dim)
 	{
 	cyclicshift=(cyclicshift+1)%dim;
-	for(int i=1; i<dim; ++i)
+	for(int i=1; i<dim-1; ++i)
 		for(int j=1; j<=i; ++j)
 			bmat(i,j)=bmat(i+1,j+1);
 	}
@@ -62,48 +62,65 @@ else
 
 //store vector
 if(incore) stor[(aktdim-1+cyclicshift)%dim]=vec;
-else st.put(vec,(aktdim-1+cyclicshift)%dim);
+else st->put(vec,(aktdim-1+cyclicshift)%dim);
 
-//calculate overlaps
-bmat(aktdim,aktdim)= vec.dot(vec);
-if(incore)
-	for(int i=1; i<aktdim; ++i) bmat(i,aktdim)=vec.dot(stor[(i-1+cyclicshift)%dim]);
+if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000;
+
+//calculate difference; 
+vec.copyonwrite();
+if(incore) vec -= stor[(aktdim-2+cyclicshift)%dim];
 else
 	{
-	T tmp=vec; //copy dimensions
-	for(int i=1; i<aktdim; ++i) 
+	T tmp=vec;
+	st->get(tmp,(aktdim-2+cyclicshift)%dim);
+	vec -= tmp;
+	}
+
+//calculate overlaps of differences (if storage is cheap, they could rather be stored than recomputed)
+typename LA_traits<T>::normtype norm=vec.norm();
+bmat(aktdim-1,aktdim-1)= norm*norm;
+if(incore)
+	for(int i=1; i<aktdim-1; ++i) 
+		bmat(i,aktdim-1)=vec.dot(stor[(i+cyclicshift)%dim] - stor[(i-1+cyclicshift)%dim]);
+else
+	{
+	T tmp=vec;
+	T tmp2=vec; //copy dimensions
+	st->get(tmp2,(0+cyclicshift)%dim);
+	for(int i=1; i<aktdim-1; ++i) 
 		{
-		st.get(tmp,(i-1+cyclicshift)%dim);
-		bmat(i,aktdim)=vec.dot(tmp);
+		st->get(tmp,(i+cyclicshift)%dim);
+		tmp2 -= tmp;
+		bmat(i,aktdim-1)= -vec.dot(tmp2);
+		tmp2=tmp;
 		}
 	}
 
 //prepare rhs-solution vector
-NRVec<Te> rhs(dim+1);
+NRVec<Te> rhs(dim);
 rhs= (Te)0; rhs[0]= (Te)-1;
 
 //solve for coefficients
 {
 NRSMat<Te> amat=bmat;
-linear_solve(amat,rhs,NULL,aktdim+1);
+linear_solve(amat,rhs,NULL,aktdim);
 }
 
 //build the new linear combination
-vec.copyonwrite();
-vec *= rhs[aktdim];
+vec = (Te)0;
 if(incore)
-	for(int i=1; i<aktdim; ++i) vec.axpy(rhs[i],stor[(i-1+cyclicshift)%dim]);
+	for(int i=1; i<aktdim; ++i) vec.axpy(rhs[i],stor[(i+cyclicshift)%dim]);
 else
         {
 	T tmp=vec; //copy dimensions
 	for(int i=1; i<aktdim; ++i)
                 {
-                st.get(tmp,(i-1+cyclicshift)%dim);
+                st->get(tmp,(i+cyclicshift)%dim);
 		vec.axpy(rhs[i],tmp);
                 }
         }
 
-return rhs[0];
+return norm;
 }