19 #ifndef _libint2_src_lib_libint_diis_h_ 20 #define _libint2_src_lib_libint_diis_h_ 22 #include <libint2/cxxstd.h> 23 #if LIBINT2_CPLUSPLUS_STD < 2011 24 # error "libint2/diis.h requires C++11 support" 53 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
54 struct traits<Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >> {
55 typedef Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > D;
56 typedef typename D::Scalar element_type;
59 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
61 dot_product(
const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d1,
62 const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d2) {
63 return d1.cwiseProduct(d2).sum();
66 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols>
68 zero(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d) {
69 d.setZero(d.rows(), d.cols());
72 template<
typename _Scalar,
int _Rows,
int _Cols,
int _Options,
int _MaxRows,
int _MaxCols,
typename Scalar>
74 axpy(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& y,
76 const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& x) {
126 template <
typename D>
155 unsigned int ngrdiis=1,
157 error_(0), errorset_(false),
158 start(strt), ndiis(ndi),
159 iter(0), ngroup(ngr),
181 bool extrapolate_error =
false)
185 const value_type zero_determinant = std::numeric_limits<value_type>::epsilon();
186 const value_type zero_norm = 1.0e-10;
190 const bool do_mixing = (mixing_fraction != 0.0);
191 const value_type scale = 1.0 + damping_factor;
194 if (errors_.size() == ndiis) {
197 if (not x_extrap_.empty()) x_extrap_.pop_front();
198 EigenMatrixX Bcrop = B_.bottomRightCorner(ndiis-1,ndiis-1);
199 Bcrop.conservativeResize(ndiis,ndiis);
205 errors_.push_back(error);
206 const unsigned int nvec = errors_.size();
207 assert(x_.size() == errors_.size());
210 for (
unsigned int i=0; i < nvec-1; i++)
211 B_(i,nvec-1) = B_(nvec-1,i) = dot_product(errors_[i], errors_[nvec-1]);
212 B_(nvec-1,nvec-1) = dot_product(errors_[nvec-1], errors_[nvec-1]);
215 if (not x_extrap_.empty() && do_mixing) {
217 axpy(x, (1.0-mixing_fraction), x_[0]);
218 axpy(x, mixing_fraction, x_extrap_[0]);
221 else if (iter > start && (((iter - start) % ngroup) < ngroupdiis)) {
226 unsigned int nskip = 0;
231 const unsigned int rank = nvec - nskip + 1;
234 EigenMatrixX A(rank, rank);
235 A.col(0).setConstant(-1.0);
236 A.row(0).setConstant(-1.0);
238 EigenVectorX rhs = EigenVectorX::Zero(rank);
241 value_type norm = 1.0;
242 if (B_(nskip,nskip) > zero_norm)
243 norm = 1.0/B_(nskip,nskip);
245 A.block(1, 1, rank-1, rank-1) = B_.block(nskip, nskip, rank-1, rank-1) * norm;
246 A.diagonal() *= scale;
255 std::cout <<
"DIIS: iter=" << iter <<
" nskip=" << nskip <<
" nvec=" << nvec << std::endl;
256 std::cout <<
"DIIS: B=" << B_ << std::endl;
257 std::cout <<
"DIIS: A=" << A << std::endl;
258 std::cout <<
"DIIS: rhs=" << rhs << std::endl;
262 Eigen::ColPivHouseholderQR<EigenMatrixX> A_QR = A.colPivHouseholderQr();
264 absdetA = A_QR.absDeterminant();
270 }
while (absdetA < zero_determinant && nskip < nvec);
273 if (absdetA < zero_determinant) {
274 std::ostringstream oss;
275 oss <<
"DIIS::extrapolate: poorly-conditioned system, |A| = " << absdetA;
276 throw std::domain_error(oss.str());
283 for (
unsigned int k=nskip, kk=1; k < nvec; ++k, ++kk) {
284 if (not do_mixing || x_extrap_.empty()) {
286 axpy(x, c[kk], x_[k]);
287 if (extrapolate_error)
288 axpy(error, c[kk], errors_[k]);
290 axpy(x, c[kk] * (1.0 - mixing_fraction), x_[k]);
291 axpy(x, c[kk] * mixing_fraction, x_extrap_[k]);
298 if (do_mixing) x_extrap_.push_back(x);
305 if (start > iter) start = iter+1;
308 void reinitialize(
const D* data = 0) {
311 const bool do_mixing = (mixing_fraction != 0.0);
312 if (do_mixing) x_extrap_.push_front(*data);
324 unsigned int ngroupdiis;
325 value_type damping_factor;
326 value_type mixing_fraction;
328 typedef Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EigenMatrixX;
329 typedef Eigen::Matrix<value_type, Eigen::Dynamic, 1> EigenVectorX;
334 std::deque<D> errors_;
335 std::deque<D> x_extrap_;
337 void set_error(value_type e) { error_ = e; errorset_ =
true; }
338 value_type error() {
return error_; }
343 B_ = EigenMatrixX::Zero(ndiis,ndiis);
359 #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:304
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
DIIS (``direct inversion of iterative subspace'') extrapolation.
Definition: diis.h:127
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:151
void extrapolate(D &x, D &error, bool extrapolate_error=false)
Definition: diis.h:179