LIBINT  2.6.0
diis.h
1 /*
2  * Copyright (C) 2004-2019 Edward F. Valeev
3  *
4  * This file is part of Libint.
5  *
6  * Libint is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Libint is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_lib_libint_diis_h_
22 #define _libint2_src_lib_libint_diis_h_
23 
24 #include <libint2/util/cxxstd.h>
25 #if LIBINT2_CPLUSPLUS_STD < 2011
26 # error "libint2/diis.h requires C++11 support"
27 #endif
28 
29 #include <deque>
30 
31 #pragma GCC diagnostic push
32 #pragma GCC system_header
33 #include <Eigen/Core>
34 #include <Eigen/QR>
35 #pragma GCC diagnostic pop
36 
37 namespace libint2 {
38 
39  namespace diis {
40 
41  template <typename D>
42  struct traits;
43 
44  /*
45  template <typename D>
46  typename traits<D>::element_type
47  dot_product(const D& d1, const D& d2);
48 
49  template <typename D>
50  void
51  zero(D& d);
52 
53  template <typename D, typename Scalar>
54  void
55  axpy(const D& y, Scalar a, const D& x);
56  */
57 
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;
62  };
63 
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();
69  }
70 
71  template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
72  void
73  zero(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& d) {
74  d.setZero(d.rows(), d.cols());
75  }
76 
77  template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols, typename Scalar>
78  void
79  axpy(Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& y,
80  Scalar a,
81  const Eigen::Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols >& x) {
82  y += a*x;
83  }
84 
85  }
86 
88 
131  template <typename D>
132  class DIIS {
133  public:
134  typedef typename diis::traits<D>::element_type value_type;
135 
137 
156  DIIS(unsigned int strt=1,
157  unsigned int ndi=5,
158  value_type dmp =0,
159  unsigned int ngr=1,
160  unsigned int ngrdiis=1,
161  value_type mf=0) :
162  error_(0), errorset_(false),
163  start(strt), ndiis(ndi),
164  iter(0), ngroup(ngr),
165  ngroupdiis(ngr),
166  damping_factor(dmp),
167  mixing_fraction(mf)
168  {
169  init();
170  }
171  ~DIIS() {
172  x_.clear();
173  errors_.clear();
174  x_extrap_.clear();
175  }
176 
184  void extrapolate(D& x,
185  D& error,
186  bool extrapolate_error = false)
187  {
188  using namespace ::libint2::diis;
189 
190  const value_type zero_determinant = std::numeric_limits<value_type>::epsilon();
191  const value_type zero_norm = 1.0e-10;
192 
193  iter++;
194 
195  const bool do_mixing = (mixing_fraction != 0.0);
196  const value_type scale = 1.0 + damping_factor;
197 
198  // if have ndiis vectors
199  if (errors_.size() == ndiis) { // holding max # of vectors already? drop the least recent {x, error} pair
200  x_.pop_front();
201  errors_.pop_front();
202  if (not x_extrap_.empty()) x_extrap_.pop_front();
203  EigenMatrixX Bcrop = B_.bottomRightCorner(ndiis-1,ndiis-1);
204  Bcrop.conservativeResize(ndiis,ndiis);
205  B_ = Bcrop;
206  }
207 
208  // push {x, error} to the set
209  x_.push_back(x);
210  errors_.push_back(error);
211  const unsigned int nvec = errors_.size();
212  assert(x_.size() == errors_.size());
213 
214  // and compute the most recent elements of B, B(i,j) = <ei|ej>
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]);
218 
219  if (iter == 1) { // the first iteration
220  if (not x_extrap_.empty() && do_mixing) {
221  zero(x);
222  axpy(x, (1.0-mixing_fraction), x_[0]);
223  axpy(x, mixing_fraction, x_extrap_[0]);
224  }
225  }
226  else if (iter > start && (((iter - start) % ngroup) < ngroupdiis)) { // not the first iteration and need to extrapolate?
227 
228  EigenVectorX c;
229 
230  value_type absdetA;
231  unsigned int nskip = 0; // how many oldest vectors to skip for the sake of conditioning?
232  // try zero
233  // skip oldest vectors until found a numerically stable system
234  do {
235 
236  const unsigned int rank = nvec - nskip + 1; // size of matrix A
237 
238  // set up the DIIS linear system: A c = rhs
239  EigenMatrixX A(rank, rank);
240  A.col(0).setConstant(-1.0);
241  A.row(0).setConstant(-1.0);
242  A(0,0) = 0.0;
243  EigenVectorX rhs = EigenVectorX::Zero(rank);
244  rhs[0] = -1.0;
245 
246  value_type norm = 1.0;
247  if (B_(nskip,nskip) > zero_norm)
248  norm = 1.0/B_(nskip,nskip);
249 
250  A.block(1, 1, rank-1, rank-1) = B_.block(nskip, nskip, rank-1, rank-1) * norm;
251  A.diagonal() *= scale;
252  //for (unsigned int i=1; i < rank ; i++) {
253  // for (unsigned int j=1; j <= i ; j++) {
254  // A(i, j) = A(j, i) = B_(i+nskip-1, j+nskip-1) * norm;
255  // if (i==j) A(i, j) *= scale;
256  // }
257  //}
258 
259 #if 0
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;
264 #endif
265 
266  // finally, solve the DIIS linear system
267  Eigen::ColPivHouseholderQR<EigenMatrixX> A_QR = A.colPivHouseholderQr();
268  c = A_QR.solve(rhs);
269  absdetA = A_QR.absDeterminant();
270 
271  //std::cout << "DIIS: |A|=" << absdetA << " sol=" << c << std::endl;
272 
273  ++nskip;
274 
275  } while (absdetA < zero_determinant && nskip < nvec); // while (system is poorly conditioned)
276 
277  // failed?
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());
282  }
283  --nskip; // undo the last ++ :-(
284 
285  {
286 
287  zero(x);
288  for (unsigned int k=nskip, kk=1; k < nvec; ++k, ++kk) {
289  if (not do_mixing || x_extrap_.empty()) {
290  //std::cout << "contrib " << k << " c=" << c[kk] << ":" << std::endl << x_[k] << std::endl;
291  axpy(x, c[kk], x_[k]);
292  if (extrapolate_error)
293  axpy(error, c[kk], errors_[k]);
294  } else {
295  axpy(x, c[kk] * (1.0 - mixing_fraction), x_[k]);
296  axpy(x, c[kk] * mixing_fraction, x_extrap_[k]);
297  }
298  }
299  }
300  } // do DIIS
301 
302  // only need to keep extrapolated x if doing mixing
303  if (do_mixing) x_extrap_.push_back(x);
304  }
305 
310  if (start > iter) start = iter+1;
311  }
312 
313  void reinitialize(const D* data = 0) {
314  iter=0;
315  if (data) {
316  const bool do_mixing = (mixing_fraction != 0.0);
317  if (do_mixing) x_extrap_.push_front(*data);
318  }
319  }
320 
321  private:
322  value_type error_;
323  bool errorset_;
324 
325  unsigned int start;
326  unsigned int ndiis;
327  unsigned int iter;
328  unsigned int ngroup;
329  unsigned int ngroupdiis;
330  value_type damping_factor;
331  value_type mixing_fraction;
332 
333  typedef Eigen::Matrix<value_type, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EigenMatrixX;
334  typedef Eigen::Matrix<value_type, Eigen::Dynamic, 1> EigenVectorX;
335 
336  EigenMatrixX B_;
337 
338  std::deque<D> x_;
339  std::deque<D> errors_;
340  std::deque<D> x_extrap_;
341 
342  void set_error(value_type e) { error_ = e; errorset_ = true; }
343  value_type error() { return error_; }
344 
345  void init() {
346  iter = 0;
347 
348  B_ = EigenMatrixX::Zero(ndiis,ndiis);
349 
350  x_.clear();
351  errors_.clear();
352  x_extrap_.clear();
353  //x_.resize(ndiis);
354  //errors_.resize(ndiis);
355  // x_extrap_ is bigger than the other because
356  // it must hold data associated with the next iteration
357  //x_extrap_.resize(diis+1);
358  }
359 
360  }; // class DIIS
361 
362 } // namespace libint2
363 
364 #include <libint2/engine.h>
365 
366 #endif /* _libint2_src_lib_libint_diis_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
Definition: diis.h:42
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