//DIIS convergence acceleration according to Pulay: Chem. Phys. Lett. 73, 393 (1980); J. Comp. Chem. 3,556 (1982) #ifndef _DIIS_H_ #define _DIIS_H_ #include "vec.h" #include "smat.h" #include "mat.h" #include "sparsemat.h" #include "nonclass.h" #include "la_traits.h" #include "auxstorage.h" // Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex fields // actually it can be anything what has operator=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits::normtype and elementtype // and get() and put() if external storage is requested template class DIIS { int dim; int aktdim; bool incore; int cyclicshift; //circular buffer of last dim vectors typedef typename LA_traits::elementtype Te; NRSMat bmat; AuxStorage *st; T *stor; public: DIIS() {dim=0;}; //for array of diis DIIS(const int n, const bool core=1); void setup(const int n, const bool core=1); ~DIIS(); typename LA_traits::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm }; template DIIS::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n) { st=incore?NULL: new AuxStorage; stor= incore? new T[dim] : NULL; bmat= (Te)0; for(int i=1; i void DIIS::setup(const int n, const bool core) { dim=n; incore=core; bmat.resize(n); st=incore?NULL: new AuxStorage; stor= incore? new T[dim] : NULL; bmat= (Te)0; for(int i=1; i DIIS::~DIIS() { if(st) delete st; if(stor) delete[] stor; } template typename LA_traits::normtype DIIS::extrapolate(T &vec) { if(!dim) laerror("attempt to extrapolate from uninitialized DIIS"); //if dim exceeded, shift if(aktdim==dim) { cyclicshift=(cyclicshift+1)%dim; for(int i=1; iput(vec,(aktdim-1+cyclicshift)%dim); if(aktdim==1) return (typename LA_traits::normtype)1000000000; //calculate difference; vec.copyonwrite(); if(incore) vec.axpy((Te)-1,stor[(aktdim-2+cyclicshift)%dim]); else { T tmp=vec; st->get(tmp,(aktdim-2+cyclicshift)%dim); vec.axpy((Te)-1,tmp); } //calculate overlaps of differences (if storage is cheap, they could rather be stored than recomputed) typename LA_traits::normtype norm=vec.norm(); bmat(aktdim-1,aktdim-1)= norm*norm; if(incore) for(int i=1; iget(tmp2,(0+cyclicshift)%dim); for(int i=1; iget(tmp,(i+cyclicshift)%dim); tmp2 -= tmp; bmat(i,aktdim-1)= -vec.dot(tmp2); tmp2=tmp; } } //prepare rhs-solution vector NRVec rhs(dim); rhs= (Te)0; rhs[0]= (Te)-1; //solve for coefficients { NRSMat amat=bmat; linear_solve(amat,rhs,NULL,aktdim); } //build the new linear combination vec.clear(); if(incore) for(int i=1; iget(tmp,(i+cyclicshift)%dim); vec.axpy(rhs[i],tmp); } } return norm; } #endif