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-03 20:17:57 +00:00
parent e6fd8cc2e8
commit 57971235dd
1 changed files with 34 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
diis.h
View File

@ -9,6 +9,9 @@
#include "la_traits.h" #include "la_traits.h"
#include "auxstorage.h" #include "auxstorage.h"
//Pulay memorial book remarks - for numerical stabilization small addition to diagonal
#define DIISEPS 1e-9
// Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields // Typically, T is some solution vector in form of NRVec, NRMat, or NRSMat over double or complex<double> fields
// actually it can be anything what has operator=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype // actually it can be anything what has operator=(const T&), clear(), dot() , axpy(), norm() and copyonwrite(), and LA_traits<T>::normtype and elementtype
// and get() and put() if external storage is requested // and get() and put() if external storage is requested
@ -22,21 +25,23 @@ class DIIS
int cyclicshift; //circular buffer of last dim vectors int cyclicshift; //circular buffer of last dim vectors
typedef typename LA_traits<T>::elementtype Te; typedef typename LA_traits<T>::elementtype Te;
NRSMat<Te> bmat; NRSMat<Te> bmat;
AuxStorage<Te> *st; AuxStorage<Te> *st, *errst;
T *stor; T *stor, *errstor;
public: public:
DIIS() {dim=0;}; //for array of diis DIIS() {dim=0; st=NULL; stor=NULL;}; //for array of diis
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<T>::normtype extrapolate(T &vec); //vec is input/output; returns square residual norm typename LA_traits<T>::normtype extrapolate(T &vec, const T &errvec); //vec is input/output; returns square residual norm
}; };
template<typename T> template<typename T>
DIIS<T>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n) DIIS<T>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n,n)
{ {
st=incore?NULL: new AuxStorage<Te>; st=incore?NULL: new AuxStorage<Te>;
errst=incore?NULL: new AuxStorage<Te>;
stor= incore? new T[dim] : NULL; stor= incore? new T[dim] : NULL;
errstor= 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; aktdim=cyclicshift=0;
} }
@ -48,7 +53,9 @@ dim=n;
incore=core; incore=core;
bmat.resize(n); bmat.resize(n);
st=incore?NULL: new AuxStorage<Te>; st=incore?NULL: new AuxStorage<Te>;
errst=incore?NULL: new AuxStorage<Te>;
stor= incore? new T[dim] : NULL; stor= incore? new T[dim] : NULL;
errstor= 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; aktdim=cyclicshift=0;
} }
@ -58,12 +65,15 @@ template<typename T>
DIIS<T>::~DIIS() DIIS<T>::~DIIS()
{ {
if(st) delete st; if(st) delete st;
if(errst) delete errst;
if(stor) delete[] stor; if(stor) delete[] stor;
if(errstor) delete[] errstor;
} }
template<typename T> template<typename T>
typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec) typename LA_traits<T>::normtype DIIS<T>::extrapolate(T &vec, const T &errvec)
{ {
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
@ -78,46 +88,44 @@ else
++aktdim; ++aktdim;
//store vector //store vector
if(incore) stor[(aktdim-1+cyclicshift)%dim]=vec; if(incore)
else st->put(vec,(aktdim-1+cyclicshift)%dim); {
stor[(aktdim-1+cyclicshift)%dim]=vec;
errstor[(aktdim-1+cyclicshift)%dim]=errvec;
}
else
{
st->put(vec,(aktdim-1+cyclicshift)%dim);
errst->put(errvec,(aktdim-1+cyclicshift)%dim);
}
if(aktdim==1) return (typename LA_traits<T>::normtype)1000000000; if(aktdim==1) return (typename LA_traits<T>::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) //calculate overlaps of the new error with old ones
typename LA_traits<T>::normtype norm=vec.norm(); typename LA_traits<T>::normtype norm=errvec.norm();
bmat(aktdim-1,aktdim-1)= norm*norm; bmat(aktdim-1,aktdim-1)= norm*norm + DIISEPS;
if(incore) if(incore)
for(int i=1; i<aktdim-1; ++i) for(int i=1; i<aktdim-1; ++i)
bmat(i,aktdim-1)=vec.dot(stor[(i+cyclicshift)%dim] - stor[(i-1+cyclicshift)%dim]); bmat(i,aktdim-1)=errvec.dot(errstor[(i+cyclicshift)%dim]);
else else
{ {
T tmp=vec; T tmp = vec; tmp.copyonwrite(); //copy dimensions
T tmp2=vec; //copy dimensions
st->get(tmp2,(0+cyclicshift)%dim);
for(int i=1; i<aktdim-1; ++i) for(int i=1; i<aktdim-1; ++i)
{ {
st->get(tmp,(i+cyclicshift)%dim); st->get(tmp,(i+cyclicshift)%dim);
tmp2 -= tmp; bmat(i,aktdim-1)= errvec.dot(tmp);
bmat(i,aktdim-1)= -vec.dot(tmp2);
tmp2=tmp;
} }
} }
//prepare rhs-solution vector //prepare rhs-solution vector
NRVec<Te> rhs(dim); NRVec<Te> rhs(dim);
rhs= (Te)0; rhs[0]= (Te)-1; rhs= (Te)0; rhs[0]= (Te)-1;
//solve for coefficients //solve for coefficients
//@@@@@@ implement checking for bad condition number and eliminating old vectors
//@@@ 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);