//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"
namespaceLA{
//Pulay memorial book remarks - for numerical stabilization small addition to diagonal (but our experience was opposite)
// 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
// and get() and put() if external storage is requested
template<typenameT,typenameU>
classDIIS
{
intdim;
intaktdim;
boolincore;
intcyclicshift;//circular buffer of last dim vectors
typedeftypenameLA_traits<T>::elementtypeTe;
typedeftypenameLA_traits<U>::elementtypeUe;
typedeftypenameLA_traits<U>::normtypeUn;
NRSMat<Ue>bmat;
AuxStorage<Te>*st;
AuxStorage<Ue>*errst;
T*stor;
U*errstor;
public:
DIIS(){dim=0;st=NULL;stor=NULL;errst=NULL;errstor=NULL;};//for array of diis
//@@@@@@or use DGGLSE, interface that in nonclass.cc, the matrix then has to be without the 0-th dimension - will be a bit confusing
//@@@@@this avoid the normal equations (which worsen the condition number), but requires work with matrices with one LARGE dimension and can be inconvenient
//for several approaches to DIIS solution cf. Shepard et al. Mol. Phys. 105, 2839 (2007)