00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef _DIIS_H_
00020 #define _DIIS_H_
00021 #include "vec.h"
00022 #include "smat.h"
00023 #include "mat.h"
00024 #include "sparsemat.h"
00025 #include "nonclass.h"
00026 #include "la_traits.h"
00027 #include "auxstorage.h"
00028
00029 namespace LA {
00030
00031
00032 const double DIISEPS=1e-9;
00033
00034
00035
00036
00037
00038 template<typename T, typename U>
00039 class DIIS
00040 {
00041 int dim;
00042 int aktdim;
00043 bool incore;
00044 int cyclicshift;
00045 typedef typename LA_traits<T>::elementtype Te;
00046 typedef typename LA_traits<U>::elementtype Ue;
00047 NRSMat<Ue> bmat;
00048 AuxStorage<Te> *st;
00049 AuxStorage<Ue> *errst;
00050 T *stor;
00051 U *errstor;
00052 public:
00053 DIIS() {dim=0; st=NULL; stor=NULL; errst=NULL; errstor=NULL;};
00054 DIIS(const int n, const bool core=1);
00055 void setup(const int n, const bool core=1);
00056 ~DIIS();
00057 typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec, bool verbose=false);
00058 };
00059
00060 template<typename T, typename U>
00061 DIIS<T,U>::DIIS(const int n, const bool core) : dim(n), incore(core), bmat(n+1,n+1)
00062 {
00063 st=incore?NULL: new AuxStorage<Te>;
00064 errst=incore?NULL: new AuxStorage<Ue>;
00065 stor= incore? new T[dim] : NULL;
00066 errstor= incore? new U[dim] : NULL;
00067 bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1;
00068 aktdim=cyclicshift=0;
00069 }
00070
00071 template<typename T, typename U>
00072 void DIIS<T,U>::setup(const int n, const bool core)
00073 {
00074 dim=n;
00075 incore=core;
00076 bmat.resize(n+1);
00077 st=incore?NULL: new AuxStorage<Te>;
00078 errst=incore?NULL: new AuxStorage<Ue>;
00079 stor= incore? new T[dim] : NULL;
00080 errstor= incore? new U[dim] : NULL;
00081 bmat= (Ue)0; for(int i=1; i<=n; ++i) bmat(0,i) = (Ue)-1;
00082 aktdim=cyclicshift=0;
00083 }
00084
00085
00086 template<typename T, typename U>
00087 DIIS<T,U>::~DIIS()
00088 {
00089 if(st) delete st;
00090 if(errst) delete errst;
00091 if(stor) delete[] stor;
00092 if(errstor) delete[] errstor;
00093 }
00094
00095
00096
00097 template<typename T, typename U>
00098 typename LA_traits<U>::normtype DIIS<T,U>::extrapolate(T &vec, const U &errvec, bool verbose)
00099 {
00100 if(!dim) laerror("attempt to extrapolate from uninitialized DIIS");
00101
00102 if(aktdim==dim)
00103 {
00104 cyclicshift=(cyclicshift+1)%dim;
00105 for(int i=1; i<dim; ++i)
00106 for(int j=1; j<=i; ++j)
00107 bmat(i,j)=bmat(i+1,j+1);
00108 }
00109 else
00110 ++aktdim;
00111
00112
00113 if(incore)
00114 {
00115 stor[(aktdim-1+cyclicshift)%dim]|=vec;
00116 errstor[(aktdim-1+cyclicshift)%dim]|=errvec;
00117 }
00118 else
00119 {
00120 st->put(vec,(aktdim-1+cyclicshift)%dim);
00121 errst->put(errvec,(aktdim-1+cyclicshift)%dim);
00122 }
00123
00124 if(aktdim==1) return (typename LA_traits<T>::normtype)1;
00125
00126
00127
00128 typename LA_traits<T>::normtype norm=errvec.norm();
00129 bmat(aktdim,aktdim)= norm*norm + DIISEPS;
00130 if(incore)
00131 for(int i=1; i<aktdim; ++i)
00132 bmat(i,aktdim)=errvec.dot(errstor[(i+cyclicshift-1)%dim]);
00133 else
00134 {
00135 U tmp = errvec; tmp.copyonwrite();
00136 for(int i=1; i<aktdim; ++i)
00137 {
00138 errst->get(tmp,(i-1+cyclicshift)%dim);
00139 bmat(i,aktdim)= errvec.dot(tmp);
00140 }
00141 }
00142
00143
00144
00145 NRVec<Ue> rhs(dim+1);
00146 rhs= (Ue)0; rhs[0]= (Ue)-1;
00147
00148
00149
00150
00151 {
00152 NRSMat<Te> amat=bmat;
00153 linear_solve(amat,rhs,NULL,aktdim+1);
00154 }
00155 if(verbose) std::cout <<"DIIS coefficients: "<<rhs<<std::endl;
00156
00157
00158 vec.clear();
00159 if(incore)
00160 for(int i=1; i<=aktdim; ++i)
00161 vec.axpy(rhs[i],stor[(i-1+cyclicshift)%dim]);
00162 else
00163 {
00164 T tmp=vec;
00165 for(int i=1; i<=aktdim; ++i)
00166 {
00167 st->get(tmp,(i-1+cyclicshift)%dim);
00168 vec.axpy(rhs[i],tmp);
00169 }
00170 }
00171
00172 return norm;
00173 }
00174
00175 }
00176
00177 #endif