ergo
PuriStepInfo.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
36 #ifndef MAT_PURISTEPINFO
37 #define MAT_PURISTEPINFO
38 #include <math.h>
39 #include <iomanip>
40 #include "PuriStepInfoDebug.h"
41 #define __ID__ "$Id: PuriStepInfo.h 4437 2012-07-05 09:01:18Z elias $"
42 namespace mat {
43  /* IDEA: Adjust threshold (looser) in later iteration if threshold in early
44  * iteration was tighter than expected.
45  */
46 
54  template<typename Treal, typename Tvector, typename TdebugPolicy>
55  class PuriStepInfo : public PuriStepInfoDebug<Treal, TdebugPolicy> {
56  public:
57  explicit PuriStepInfo(int nn = -1, int noc = -1, Treal eigvalConvCrit = 0.0)
58  : n(nn),nocc(noc), traceX(-1.0), traceX2(-1.0), poly(-1),
59  chosenThresh(0.0), actualThresh(0.0),
61  eigInterval(-1.0,2.0), delta(-1.0), correctOccupation(0),
62  XmX2EuclNorm(0.0, 10.0), eigVecPtr(0),
64  n0(0.0), n1(0.0),
65  homo(0.0, 1.0), lumo(0.0, 1.0), eigConvCrit(eigvalConvCrit),
66  nnzX(0), nnzX2(0), eigAccLoss(0),
68  timeXX2Write(0), timeXX2Read(0) {}
69 
71  delete eigVecPtr;
72  }
73  bool converged() const {
74  bool homoLumo_converged = (1 - homo.low() < eigConvCrit &&
75  lumo.upp() < eigConvCrit);
76  bool XmX2norm_converged = getXmX2EuclNorm().upp() < eigConvCrit;
77  // FIXME: Possibly, propagating XmX2norm between the
78  // iterations can give more accurate values since this norm is not
79  // computed accurately in each iteration.
80  return homoLumo_converged || XmX2norm_converged;
81  }
82 
83  inline void setChosenThresh(Treal const thr) {chosenThresh = thr;}
84  inline Treal getChosenThresh() const { return chosenThresh;}
85  inline void setActualThresh(Treal const thr) {actualThresh = thr;}
86  inline Treal getActualThresh() const { return actualThresh;}
87  inline void setEstimatedStepsLeft(int const stepsleft) {
88  estimatedStepsLeft = stepsleft;
89  }
90  inline int getEstimatedStepsLeft() const { return estimatedStepsLeft;}
91  inline void setTraceX(Treal const trX) {traceX = trX;}
92  inline void setTraceX2(Treal const trX2) {traceX2 = trX2;}
93  inline Treal getTraceX() const {return traceX;}
94  inline Treal getTraceX2() const {return traceX2;}
95  void setPoly();
96  inline int getPoly() const {return poly;}
98  inline void setXmX2EuclNorm(Interval<Treal> const XmX2eucl) {
99  XmX2EuclNorm = XmX2eucl;
100  }
101 
107  void setEigVecPtr(Tvector * eigVecPtr_) {
108  delete eigVecPtr;
109  eigVecPtr = eigVecPtr_;
110  }
111  Tvector const * const getEigVecPtr() const {return eigVecPtr;}
112 
113  void improveHomoLumo(Interval<Treal> const homoInt,
114  Interval<Treal> const lumoInt);
115  inline Interval<Treal> const & getHomo() const {
116  return homo;
117  }
118  inline Interval<Treal> const & getLumo() const {
119  return lumo;
120  }
121  inline Interval<Treal> const & getEigInterval() const {
122  return eigInterval;
123  }
128  inline int getCorrectOccupation() const {return correctOccupation;}
129  Treal subspaceError() const;
133  void improveEigInterval(Interval<Treal> const eInt);
134 
135  inline void setNnzX(size_t const nzX) {nnzX = nzX;}
136  inline size_t getNnzX() const {return nnzX;}
137  inline void setNnzX2(size_t const nzX2) {nnzX2 = nzX2;}
138  inline size_t getNnzX2() const {return nnzX2;}
139  void computeEigAccLoss();
140  inline Treal getEigAccLoss() const {return eigAccLoss;}
141 
142  inline Treal getN0() const {return n0;}
143  inline Treal getN1() const {return n1;}
144  inline Treal getDelta() const {return delta;}
145 
146 
147  inline void checkIntervals(const char* descriptionString) const {
149  checkIntervals(eigInterval, homo, lumo, XmX2EuclNorm, descriptionString);
150  }
151  template<typename Tmatrix>
152  inline void computeExactValues(Tmatrix const & X, Tmatrix const & X2) {
154  computeExactValues(X, X2, n, nocc);
155  }
156 
157  /* Functions to set and get timings and mem usage info for the different steps: */
160  }
162  memUsageInXmX2Diff = memUsage;
163  }
164  void setTimeThresh(float t) {timeThresh = t;}
165  void setTimeSquare(float t) {timeSquare = t;}
166  void setTimeXmX2Norm(float t) {timeXmX2Norm = t;}
167  void setTimeTotal(float t) {timeTotal = t;}
168  void setTimeXX2Write(float t) {timeXX2Write = t;}
169  void setTimeXX2Read(float t) {timeXX2Read = t;}
170 
173  float getTimeThresh() {return timeThresh;}
174  float getTimeSquare() {return timeSquare;}
175  float getTimeXmX2Norm() {return timeXmX2Norm;}
176  float getTimeTotal() {return timeTotal;}
177  float getTimeXX2Write() {return timeXX2Write;}
178  float getTimeXX2Read() {return timeXX2Read;}
179 
180  inline bool homoIsAccuratelyKnown
181  (Treal accuracyLimit
184  ) const {
185  return homo.length() < accuracyLimit;
186  }
187  inline bool lumoIsAccuratelyKnown
188  (Treal accuracyLimit
191  ) const {
192  return lumo.length() < accuracyLimit;
193  }
194 
195  inline bool getLumoWasComputed() {return lumoWasComputed;}
196  inline bool getHomoWasComputed() {return homoWasComputed;}
197 
198  protected:
202  void computen0n1();
203 
204  int n;
205  int nocc;
206  /* FIXME?: Upper and lower bound on traces? */
207  Treal traceX;
208  Treal traceX2;
209  int poly;
210  Treal chosenThresh;
211  Treal actualThresh;
218  Treal delta;
230  Tvector * eigVecPtr;
235  Treal n0;
238  Treal n1;
244  Treal eigConvCrit;
247  size_t nnzX;
248  size_t nnzX2;
254  Treal eigAccLoss;
255 
256  /* Variables for the time and mem usage of various operations */
259  float timeThresh;
260  float timeSquare;
262  float timeTotal;
264  float timeXX2Read;
265  private:
266  };
267 
268 #if 1
269  template<typename Treal, typename Tvector, typename TdebugPolicy>
270  std::ostream& operator<<(std::ostream& s,
272  s<<" Trace X: "<<psi.getTraceX()
273  <<" Trace X2: "<<psi.getTraceX2()
274  <<" poly: "<<psi.getPoly()
275  <<" ||E||_2: "<<psi.getActualThresh()
276  <<" Eiginterval: "<<psi.getEigInterval()
277  <<" correctOccupation: "<<psi.getCorrectOccupation()
278  <<std::endl
279  <<" XmX2EuclNorm: "<<psi.getXmX2EuclNorm()
280  <<" n0: "<<psi.getN0()
281  <<" n1: "<<psi.getN1()
282  <<" homo: "<<psi.getHomo()
283  <<" lumo: "<<psi.getLumo()
284  <<" nnzX: "<<psi.getNnzX()
285  <<" nnzX2: "<<psi.getNnzX2();
286  return s;
287  }
288 
289 #endif
290 
291  template<typename Treal, typename Tvector, typename TdebugPolicy>
293  if (Interval<Treal>::intersect(homo,lumo).empty()) {
294  ASSERTALWAYS(homo.low() > lumo.upp());
295  if (homo.low() + lumo.upp() > 1)
296  poly = 0; /* x*x */
297  else
298  poly = 1; /* 2*x - x*x */
299  }
300  else {
301  if (template_blas_fabs(traceX2 - nocc) < template_blas_fabs(2 * traceX - traceX2 - nocc))
302  poly = 0; /* x*x */
303  else
304  poly = 1; /* 2*x - x*x */
305  }
306  }
307 
308  template<typename Treal, typename Tvector, typename TdebugPolicy>
311  Interval<Treal> const lumoInt) {
312  checkIntervals("PuriStepInfo::improveHomoLumo A0");
313  homo.intersect(homoInt);
314  checkIntervals("PuriStepInfo::improveHomoLumo A1");
315  lumo.intersect(lumoInt);
316  checkIntervals("PuriStepInfo::improveHomoLumo A2");
317  ASSERTALWAYS(!homo.empty());
318  ASSERTALWAYS(!lumo.empty());
319  if (homo.low() > 0.5 && lumo.upp() < 0.5)
320  this->setCorrectOccupation();
321 
322  if (correctOccupation && 1.0 - XmX2EuclNorm.upp() * 4.0 > 0) {
323  Interval<Treal> tmp =
324  sqrtInt((XmX2EuclNorm * (Treal)(-4.0)) + (Treal)1.0);
325  ASSERTALWAYS(tmp.length() > 0);
326 
327  ASSERTALWAYS(!homo.empty());
328  ASSERTALWAYS(!lumo.empty());
329 
330  homo.intersect(Interval<Treal>((1.0 + tmp.low()) / 2.0, homo.upp()));
331 
332  ASSERTALWAYS(!homo.empty());
333  ASSERTALWAYS(!lumo.empty());
334 
335  lumo.intersect(Interval<Treal>(lumo.low(), (1.0 - tmp.low()) / 2.0));
336  checkIntervals("PuriStepInfo::improveHomoLumo B");
337  ASSERTALWAYS(!homo.empty());
338  ASSERTALWAYS(!lumo.empty());
339  if (!eigInterval.cover((1 + template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2) &&
340  !eigInterval.cover((1 - template_blas_sqrt(1.0 + 4.0 * XmX2EuclNorm.low())) / 2)){
341  /* Either homo lies in homoTmp or lumo lies in lumoTmp. */
342  Interval<Treal> homoTmp =
343  (tmp + (Treal)1.0) / (Treal)2.0;
344  Interval<Treal> lumoTmp =
345  ((tmp * (Treal)(-1.0)) + (Treal)1.0) / (Treal)2.0;
346  ASSERTALWAYS(!(Interval<Treal>::intersect(homo, homoTmp).empty() &&
347  Interval<Treal>::intersect(lumo, lumoTmp).empty()));
348  if (Interval<Treal>::intersect(homo, homoTmp).empty()) {
349  // ok, lumo was computed in this iteration
350  if (eigVecPtr)
351  lumoWasComputed = 1;
352  lumo.intersect(lumoTmp);
353  }
354  if (Interval<Treal>::intersect(lumo, lumoTmp).empty()) {
355  // ok, homo was computed in this iteration
356  if (eigVecPtr)
357  homoWasComputed = 1;
358  homo.intersect(homoTmp);
359  }
360  checkIntervals("PuriStepInfo::improveHomoLumo C");
361  }
362  }
363  }
364 
365  template<typename Treal, typename Tvector, typename TdebugPolicy>
368  Interval<Treal> zeroOneInt(0.0,1.0);
369 
370  next.checkIntervals("PuriStepInfo::exchangeInfoWithNext A");
371 
372  /* Improve homo/lumo and eig-bounds for next */
373  Interval<Treal> homoForNext(homo);
374  Interval<Treal> lumoForNext(lumo);
375  Interval<Treal> eigIntForNext(eigInterval);
376  homoForNext.puriStep(poly);
377  lumoForNext.puriStep(poly);
378  eigIntForNext.puriStep(poly);
379  ASSERTALWAYS(!homoForNext.empty());
380  ASSERTALWAYS(!lumoForNext.empty());
381  ASSERTALWAYS(!eigIntForNext.empty());
382  /* Increase intervals because relative precision in matrix-matrix
383  * multiplication can result in loss of accuracy in eigenvalues
384  */
385  homoForNext.increase(eigAccLoss);
386  lumoForNext.increase(eigAccLoss);
387  eigIntForNext.increase(eigAccLoss);
388  /* Increase intervals because of thresholding */
389  homoForNext.increase(next.actualThresh);
390  lumoForNext.increase(next.actualThresh);
391  homoForNext.intersect(zeroOneInt);
392  lumoForNext.intersect(zeroOneInt);
393  eigIntForNext.increase(next.actualThresh);
394  next.improveEigInterval(eigIntForNext);
395  next.improveHomoLumo(homoForNext, lumoForNext);
396 
397  /* Improve homo/lumo for this */
398  /* FIXME: Consider improving also eigInterval from next
399  * This could possibly only be done in one end of the interval since
400  * for example information about negative eigenvalues is lost in case
401  * of an x*x step.
402  */
403  Interval<Treal> homoTmp(next.homo);
404  Interval<Treal> lumoTmp(next.lumo);
405  ASSERTALWAYS(!homoTmp.empty());
406  ASSERTALWAYS(!lumoTmp.empty());
407  /* Increase intervals because of thresholding. */
408  homoTmp.increase(next.actualThresh);
409  lumoTmp.increase(next.actualThresh);
410  homoTmp.intersect(zeroOneInt);
411  lumoTmp.intersect(zeroOneInt);
412  homoTmp.invPuriStep(poly);
413  lumoTmp.invPuriStep(poly);
414  /* Increase intervals because relative precision in matrix-matrix
415  * multiplication can result in loss of accuracy in eigenvalues
416  */
417  homoTmp.increase(eigAccLoss);
418  lumoTmp.increase(eigAccLoss);
419 
420  this->improveHomoLumo(homoTmp, lumoTmp);
421  }
422 
423  template<typename Treal, typename Tvector, typename TdebugPolicy>
425  Interval<Treal> gap = Interval<Treal>(lumo.upp(),homo.low());
426  if (actualThresh >= gap.length())
427  return 1.0; /* 1.0 means no accuracy. */
428  else {
429  Treal error = actualThresh / (gap.length() - actualThresh);
430  return error < 1.0 ? error : (Treal)1.0;
431  }
432  }
433 
434  template<typename Treal, typename Tvector, typename TdebugPolicy>
437  eigInterval.intersect(eInt);
438  Treal delta1 = eigInterval.upp() - 1;
439  Treal delta0 = -eigInterval.low();
440  delta = delta1 > delta0 ? delta1 : delta0;
441  this->computen0n1();
442  }
443 
444 
445  template<typename Treal, typename Tvector, typename TdebugPolicy>
447  Treal beta = 0.5;
448  /* Increase beta if possible */
449  if (XmX2EuclNorm.upp() < 1 / (Treal)4)
450  beta = (1 + template_blas_sqrt(1 - 4 * XmX2EuclNorm.upp())) / 2;
451  n1 = (traceX2 - delta * (1 - beta) * n -
452  (1 - delta - beta) * traceX) /
453  ((1 + 2 * delta) * (delta + beta));
454  n0 = (traceX2 + beta * (1 + delta) * n -
455  (1 + delta + beta) * traceX) /
456  ((1 + 2 * delta) * (delta + beta));
457  if (n1 > nocc -1 &&
458  n0 > n - nocc - 1)
459  correctOccupation = 1;
460  }
461 
466  template<typename Treal, typename Tvector, typename TdebugPolicy>
468  Treal nnzPerRowX = nnzX / (Treal)n;
469  Treal maxAbsErrPerElX2 = getRelPrecision<Treal>() * nnzPerRowX;
470  /* mah = max(abs(h_ij)) \approx relPrec * (nnz(X) / n)
471  * e is the exact eigenvalue of X^2, e' is the eigenvalue of the
472  * computed X^2.
473  * | e - e' | <= || H ||_2 <= || H ||_F = sqrt( sum h_ij^2 ) <=
474  * sqrt( mah^2 * nnz(X2))
475  */
476  eigAccLoss = maxAbsErrPerElX2 * template_blas_sqrt((Treal)nnzX2);
477  ASSERTALWAYS(eigAccLoss >= 0);
478  }
479 
480 
481 
482 } /* end namespace mat */
483 #undef __ID__
484 #endif