jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2007-06-06 15:14:19 +00:00
parent a80f82b6ee
commit 18a55f7ede
1 changed files with 22 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
diis.h
View File

@ -35,17 +35,17 @@ public:
DIIS(const int n, const bool core=1); DIIS(const int n, const bool core=1);
void setup(const int n, const bool core=1); void setup(const int n, const bool core=1);
~DIIS(); ~DIIS();
typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec); //vec is input/output; returns square residual norm typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec, bool verbose=false); //vec is input/output; returns square residual norm
}; };
template<typename T, typename U> template<typename T, typename U>
DIIS<T,U>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n) DIIS<T,U>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n+1,n+1)
{ {
st=incore?NULL: new AuxStorage<Te>; st=incore?NULL: new AuxStorage<Te>;
errst=incore?NULL: new AuxStorage<Ue>; errst=incore?NULL: new AuxStorage<Ue>;
stor= incore? new T[dim] : NULL; stor= incore? new T[dim] : NULL;
errstor= incore? new U[dim] : NULL; errstor= incore? new U[dim] : NULL;
bmat= (Ue)0; for(int i=1; i<n; ++i) bmat(0,i) = (Ue)-1; bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1;
aktdim=cyclicshift=0; aktdim=cyclicshift=0;
} }
@ -54,12 +54,12 @@ void DIIS<T,U>::setup(const int n, const bool core)
{ {
dim=n; dim=n;
incore=core; incore=core;
bmat.resize(n); bmat.resize(n+1);
st=incore?NULL: new AuxStorage<Te>; st=incore?NULL: new AuxStorage<Te>;
errst=incore?NULL: new AuxStorage<Ue>; errst=incore?NULL: new AuxStorage<Ue>;
stor= incore? new T[dim] : NULL; stor= incore? new T[dim] : NULL;
errstor= incore? new U[dim] : NULL; errstor= incore? new U[dim] : NULL;
bmat= (Ue)0; for(int i=1; i<n; ++i) bmat(0,i) = (Ue)-1; bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1;
aktdim=cyclicshift=0; aktdim=cyclicshift=0;
} }
@ -76,14 +76,14 @@ if(errstor) delete[] errstor;
template<typename T, typename U> template<typename T, typename U>
typename LA_traits<U>::normtype DIIS<T,U>::extrapolate(T &vec, const U &errvec) typename LA_traits<U>::normtype DIIS<T,U>::extrapolate(T &vec, const U &errvec, bool verbose)
{ {
if(!dim) laerror("attempt to extrapolate from uninitialized DIIS"); if(!dim) laerror("attempt to extrapolate from uninitialized DIIS");
//if dim exceeded, shift //if dim exceeded, shift
if(aktdim==dim) if(aktdim==dim)
{ {
cyclicshift=(cyclicshift+1)%dim; cyclicshift=(cyclicshift+1)%dim;
for(int i=1; i<dim-1; ++i) for(int i=1; i<dim; ++i)
for(int j=1; j<=i; ++j) for(int j=1; j<=i; ++j)
bmat(i,j)=bmat(i+1,j+1); bmat(i,j)=bmat(i+1,j+1);
} }
@ -93,8 +93,8 @@ else
//store vector //store vector
if(incore) if(incore)
{ {
stor[(aktdim-1+cyclicshift)%dim]=vec; stor[(aktdim-1+cyclicshift)%dim]|=vec;
errstor[(aktdim-1+cyclicshift)%dim]=errvec; errstor[(aktdim-1+cyclicshift)%dim]|=errvec;
} }
else else
{ {
@ -107,23 +107,23 @@ if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000;
//calculate overlaps of the new error with old ones //calculate overlaps of the new error with old ones
typename LA_traits<T>::normtype norm=errvec.norm(); typename LA_traits<T>::normtype norm=errvec.norm();
bmat(aktdim-1,aktdim-1)= norm*norm + DIISEPS; bmat(aktdim,aktdim)= norm*norm + DIISEPS;
if(incore) if(incore)
for(int i=1; i<aktdim-1; ++i) for(int i=1; i<aktdim; ++i)
bmat(i,aktdim-1)=errvec.dot(errstor[(i+cyclicshift)%dim]); bmat(i,aktdim)=errvec.dot(errstor[(i+cyclicshift-1)%dim]);
else else
{ {
U tmp = errvec; tmp.copyonwrite(); //copy dimensions U tmp = errvec; tmp.copyonwrite(); //copy dimensions
for(int i=1; i<aktdim-1; ++i) for(int i=1; i<aktdim; ++i)
{ {
errst->get(tmp,(i+cyclicshift)%dim); errst->get(tmp,(i-1+cyclicshift)%dim);
bmat(i,aktdim-1)= errvec.dot(tmp); bmat(i,aktdim)= errvec.dot(tmp);
} }
} }
//prepare rhs-solution vector //prepare rhs-solution vector
NRVec<Ue> rhs(dim); NRVec<Ue> rhs(dim+1);
rhs= (Ue)0; rhs[0]= (Ue)-1; rhs= (Ue)0; rhs[0]= (Ue)-1;
//solve for coefficients //solve for coefficients
@ -131,19 +131,21 @@ rhs= (Ue)0; rhs[0]= (Ue)-1;
//@@@ explicit solution - cf. remarks in Pulay memorial book //@@@ explicit solution - cf. remarks in Pulay memorial book
{ {
NRSMat<Te> amat=bmat; NRSMat<Te> amat=bmat;
linear_solve(amat,rhs,NULL,aktdim); linear_solve(amat,rhs,NULL,aktdim+1);
} }
if(verbose) cout <<"DIIS coefficients: "<<rhs<<endl;
//build the new linear combination //build the new linear combination
vec.clear(); vec.clear();
if(incore) if(incore)
for(int i=1; i<aktdim; ++i) vec.axpy(rhs[i],stor[(i+cyclicshift)%dim]); for(int i=1; i<=aktdim; ++i)
vec.axpy(rhs[i],stor[(i-1+cyclicshift)%dim]);
else else
{ {
T tmp=vec; //copy dimensions T tmp=vec; //copy dimensions
for(int i=1; i<aktdim; ++i) for(int i=1; i<=aktdim; ++i)
{ {
st->get(tmp,(i+cyclicshift)%dim); st->get(tmp,(i-1+cyclicshift)%dim);
vec.axpy(rhs[i],tmp); vec.axpy(rhs[i],tmp);
} }
} }