*** empty log message ***
This commit is contained in:
jiri
parent
820d92b85f
commit
7fd5c953ff
354
diis.h
354
diis.h
@ -1,177 +1,177 @@
|
|||||||
/*
|
/*
|
||||||
LA: linear algebra C++ interface library
|
LA: linear algebra C++ interface library
|
||||||
Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
|
Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
|
||||||
|
|
||||||
This program is free software: you can redistribute it and/or modify
|
This program is free software: you can redistribute it and/or modify
|
||||||
it under the terms of the GNU General Public License as published by
|
it under the terms of the GNU General Public License as published by
|
||||||
the Free Software Foundation, either version 3 of the License, or
|
the Free Software Foundation, either version 3 of the License, or
|
||||||
(at your option) any later version.
|
(at your option) any later version.
|
||||||
|
|
||||||
This program is distributed in the hope that it will be useful,
|
This program is distributed in the hope that it will be useful,
|
||||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||||
GNU General Public License for more details.
|
GNU General Public License for more details.
|
||||||
|
|
||||||
You should have received a copy of the GNU General Public License
|
You should have received a copy of the GNU General Public License
|
||||||
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||||
*/
|
*/
|
||||||
//DIIS convergence acceleration according to Pulay: Chem. Phys. Lett. 73, 393 (1980); J. Comp. Chem. 3,556 (1982)
|
//DIIS convergence acceleration according to Pulay: Chem. Phys. Lett. 73, 393 (1980); J. Comp. Chem. 3,556 (1982)
|
||||||
#ifndef _DIIS_H_
|
#ifndef _DIIS_H_
|
||||||
#define _DIIS_H_
|
#define _DIIS_H_
|
||||||
#include "vec.h"
|
#include "vec.h"
|
||||||
#include "smat.h"
|
#include "smat.h"
|
||||||
#include "mat.h"
|
#include "mat.h"
|
||||||
#include "sparsemat.h"
|
#include "sparsemat.h"
|
||||||
#include "nonclass.h"
|
#include "nonclass.h"
|
||||||
#include "la_traits.h"
|
#include "la_traits.h"
|
||||||
#include "auxstorage.h"
|
#include "auxstorage.h"
|
||||||
|
|
||||||
namespace LA {
|
namespace LA {
|
||||||
|
|
||||||
//Pulay memorial book remarks - for numerical stabilization small addition to diagonal
|
//Pulay memorial book remarks - for numerical stabilization small addition to diagonal (but our experience was opposite)
|
||||||
const double 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
|
|
||||||
|
template<typename T, typename U>
|
||||||
template<typename T, typename U>
|
class DIIS
|
||||||
class DIIS
|
{
|
||||||
{
|
int dim;
|
||||||
int dim;
|
int aktdim;
|
||||||
int aktdim;
|
bool incore;
|
||||||
bool incore;
|
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;
|
typedef typename LA_traits<U>::elementtype Ue;
|
||||||
typedef typename LA_traits<U>::elementtype Ue;
|
typedef typename LA_traits<U>::normtype Un;
|
||||||
NRSMat<Ue> bmat;
|
NRSMat<Ue> bmat;
|
||||||
AuxStorage<Te> *st;
|
AuxStorage<Te> *st;
|
||||||
AuxStorage<Ue> *errst;
|
AuxStorage<Ue> *errst;
|
||||||
T *stor;
|
T *stor;
|
||||||
U *errstor;
|
U *errstor;
|
||||||
public:
|
public:
|
||||||
DIIS() {dim=0; st=NULL; stor=NULL; errst=NULL; errstor=NULL;}; //for array of diis
|
DIIS() {dim=0; st=NULL; stor=NULL; errst=NULL; errstor=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<U>::normtype extrapolate(T &vec, const U &errvec, bool verbose=false); //vec is input/output; returns square residual norm
|
typename LA_traits<U>::normtype extrapolate(T &vec, const U &errvec, bool verbose=false, const Un diiseps=0); //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+1,n+1)
|
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;
|
||||||
}
|
}
|
||||||
|
|
||||||
template<typename T, typename U>
|
template<typename T, typename U>
|
||||||
void DIIS<T,U>::setup(const int n, const bool core)
|
void DIIS<T,U>::setup(const int n, const bool core)
|
||||||
{
|
{
|
||||||
dim=n;
|
dim=n;
|
||||||
incore=core;
|
incore=core;
|
||||||
bmat.resize(n+1);
|
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;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
template<typename T, typename U>
|
template<typename T, typename U>
|
||||||
DIIS<T,U>::~DIIS()
|
DIIS<T,U>::~DIIS()
|
||||||
{
|
{
|
||||||
if(st) delete st;
|
if(st) delete st;
|
||||||
if(errst) delete errst;
|
if(errst) delete errst;
|
||||||
if(stor) delete[] stor;
|
if(stor) delete[] stor;
|
||||||
if(errstor) delete[] errstor;
|
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, bool verbose, const Un diiseps)
|
||||||
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; ++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);
|
}
|
||||||
}
|
else
|
||||||
else
|
++aktdim;
|
||||||
++aktdim;
|
|
||||||
|
//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
|
{
|
||||||
{
|
st->put(vec,(aktdim-1+cyclicshift)%dim);
|
||||||
st->put(vec,(aktdim-1+cyclicshift)%dim);
|
errst->put(errvec,(aktdim-1+cyclicshift)%dim);
|
||||||
errst->put(errvec,(aktdim-1+cyclicshift)%dim);
|
}
|
||||||
}
|
|
||||||
|
if(aktdim==1) return (typename LA_traits<T>::normtype)1;
|
||||||
if(aktdim==1) return (typename LA_traits<T>::normtype)1;
|
|
||||||
|
//calculate overlaps of the new error with old ones
|
||||||
|
typename LA_traits<T>::normtype norm=errvec.norm();
|
||||||
//calculate overlaps of the new error with old ones
|
bmat(aktdim,aktdim) = norm*norm;
|
||||||
typename LA_traits<T>::normtype norm=errvec.norm();
|
// LV
|
||||||
bmat(aktdim,aktdim)= norm*norm + DIISEPS;
|
bmat(aktdim,aktdim) += diiseps;
|
||||||
if(incore)
|
|
||||||
for(int i=1; i<aktdim; ++i)
|
if(incore)
|
||||||
bmat(i,aktdim)=errvec.dot(errstor[(i+cyclicshift-1)%dim]);
|
for(int i=1; i<aktdim; ++i)
|
||||||
else
|
bmat(i,aktdim)=errvec.dot(errstor[(i+cyclicshift-1)%dim]);
|
||||||
{
|
else
|
||||||
U tmp = errvec; tmp.copyonwrite(); //copy dimensions
|
{
|
||||||
for(int i=1; i<aktdim; ++i)
|
U tmp = errvec; tmp.copyonwrite(); //copy dimensions
|
||||||
{
|
for(int i=1; i<aktdim; ++i)
|
||||||
errst->get(tmp,(i-1+cyclicshift)%dim);
|
{
|
||||||
bmat(i,aktdim)= errvec.dot(tmp);
|
errst->get(tmp,(i-1+cyclicshift)%dim);
|
||||||
}
|
bmat(i,aktdim)= errvec.dot(tmp);
|
||||||
}
|
}
|
||||||
|
}
|
||||||
|
|
||||||
//prepare rhs-solution vector
|
//prepare rhs-solution vector
|
||||||
NRVec<Ue> rhs(dim+1);
|
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
|
||||||
//@@@@@@ implement checking for bad condition number and eliminating old vectors
|
//@@@@@@ implement checking for bad condition number and eliminating old vectors
|
||||||
//@@@ 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+1);
|
linear_solve(amat,rhs,NULL,aktdim+1);
|
||||||
}
|
}
|
||||||
if(verbose) std::cout <<"DIIS coefficients: "<<rhs<<std::endl;
|
if(verbose) std::cout <<"DIIS coefficients: "<<rhs<<std::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)
|
for(int i=1; i<=aktdim; ++i)
|
||||||
vec.axpy(rhs[i],stor[(i-1+cyclicshift)%dim]);
|
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-1+cyclicshift)%dim);
|
st->get(tmp,(i-1+cyclicshift)%dim);
|
||||||
vec.axpy(rhs[i],tmp);
|
vec.axpy(rhs[i],tmp);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
return norm;
|
return norm;
|
||||||
}
|
}
|
||||||
|
|
||||||
}//namespace
|
}//namespace
|
||||||
|
|
||||||
#endif
|
#endif
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user