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