21 #ifndef _libint2_src_lib_libint_diis_h_ 22 #define _libint2_src_lib_libint_diis_h_ 24 #include <libint2/util/cxxstd.h> 25 #if LIBINT2_CPLUSPLUS_STD < 2011 26 # error "libint2/diis.h requires C++11 support" 31 #pragma GCC diagnostic push 32 #pragma GCC system_header 35 #pragma GCC diagnostic pop 58 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
59 struct traits<Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >> {
60 typedef Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > D;
61 typedef typename D::Scalar element_type;
64 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
66 dot_product(
const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d1,
67 const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d2) {
68 return d1.cwiseProduct(d2).sum();
71 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
73 zero(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d) {
74 d.setZero(d.rows(), d.cols());
77 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols,
typename Scalar>
79 axpy(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& y,
81 const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& x) {
131 template <
typename D>
160 unsigned int ngrdiis=1,
162 error_(0), errorset_(false),
163 start(strt), ndiis(ndi),
164 iter(0), ngroup(ngr),
186 bool extrapolate_error =
false)
188 using namespace ::libint2::diis;
190 const value_type zero_determinant = std::numeric_limits<value_type>::epsilon();
191 const value_type zero_norm = 1.0e-10;
195 const bool do_mixing = (mixing_fraction != 0.0);
196 const value_type scale = 1.0 + damping_factor;
199 if (errors_.size() == ndiis) {
202 if (not x_extrap_.empty()) x_extrap_.pop_front();
203 EigenMatrixX Bcrop = B_.bottomRightCorner(ndiis-1,ndiis-1);
204 Bcrop.conservativeResize(ndiis,ndiis);
210 errors_.push_back(error);
211 const unsigned int nvec = errors_.size();
212 assert(x_.size() == errors_.size());
215 for (
unsigned int i=0; i < nvec-1; i++)
216 B_(i,nvec-1) = B_(nvec-1,i) = dot_product(errors_[i], errors_[nvec-1]);
217 B_(nvec-1,nvec-1) = dot_product(errors_[nvec-1], errors_[nvec-1]);
220 if (not x_extrap_.empty() && do_mixing) {
222 axpy(x, (1.0-mixing_fraction), x_[0]);
223 axpy(x, mixing_fraction, x_extrap_[0]);
226 else if (iter > start && (((iter - start) % ngroup) < ngroupdiis)) {
231 unsigned int nskip = 0;
236 const unsigned int rank = nvec - nskip + 1;
239 EigenMatrixX A(rank, rank);
240 A.col(0).setConstant(-1.0);
241 A.row(0).setConstant(-1.0);
243 EigenVectorX rhs = EigenVectorX::Zero(rank);
246 value_type norm = 1.0;
247 if (B_(nskip,nskip) > zero_norm)
248 norm = 1.0/B_(nskip,nskip);
250 A.block(1, 1, rank-1, rank-1) = B_.block(nskip, nskip, rank-1, rank-1) * norm;
251 A.diagonal() *= scale;
260 std::cout <<
"DIIS: iter=" << iter <<
" nskip=" << nskip <<
" nvec=" << nvec << std::endl;
261 std::cout <<
"DIIS: B=" << B_ << std::endl;
262 std::cout <<
"DIIS: A=" << A << std::endl;
263 std::cout <<
"DIIS: rhs=" << rhs << std::endl;
267 Eigen::ColPivHouseholderQR<EigenMatrixX> A_QR = A.colPivHouseholderQr();
269 absdetA = A_QR.absDeterminant();
275 }
while (absdetA < zero_determinant && nskip < nvec);
278 if (absdetA < zero_determinant) {
279 std::ostringstream oss;
280 oss <<
"DIIS::extrapolate: poorly-conditioned system, |A| = " << absdetA;
281 throw std::domain_error(oss.str());
288 for (
unsigned int k=nskip, kk=1; k < nvec; ++k, ++kk) {
289 if (not do_mixing || x_extrap_.empty()) {
291 axpy(x, c[kk], x_[k]);
292 if (extrapolate_error)
293 axpy(error, c[kk], errors_[k]);
295 axpy(x, c[kk] * (1.0 - mixing_fraction), x_[k]);
296 axpy(x, c[kk] * mixing_fraction, x_extrap_[k]);
303 if (do_mixing) x_extrap_.push_back(x);
310 if (start > iter) start = iter+1;
313 void reinitialize(
const D* data = 0) {
316 const bool do_mixing = (mixing_fraction != 0.0);
317 if (do_mixing) x_extrap_.push_front(*data);
329 unsigned int ngroupdiis;
330 value_type damping_factor;
331 value_type mixing_fraction;
333 typedef Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EigenMatrixX;
334 typedef Eigen::Matrix<value_type, Eigen::Dynamic, 1> EigenVectorX;
339 std::deque<D> errors_;
340 std::deque<D> x_extrap_;
342 void set_error(value_type e) { error_ = e; errorset_ =
true; }
343 value_type error() {
return error_; }
348 B_ = EigenMatrixX::Zero(ndiis,ndiis);
364 #include <libint2/engine.h> void start_extrapolation()
calling this function forces the extrapolation to start upon next call to extrapolate() even if this ...
Definition: diis.h:309
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24
DIIS (`‘direct inversion of iterative subspace’') extrapolation.
Definition: diis.h:132
DIIS(unsigned int strt=1, unsigned int ndi=5, value_type dmp=0, unsigned int ngr=1, unsigned int ngrdiis=1, value_type mf=0)
Constructor.
Definition: diis.h:156
void extrapolate(D &x, D &error, bool extrapolate_error=false)
Definition: diis.h:184