ergo
PuriInfo.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_PURIINFO
37 #define MAT_PURIINFO
38 #include <math.h>
39 #include <iomanip>
40 #include "PuriStepInfo.h"
41 #define __ID__ "$Id$"
42 namespace mat {
46  template<typename Treal, typename Tvector, typename TdebugPolicy>
47  class PuriInfo : public TdebugPolicy {
48  public:
49  PuriInfo(int nn, int noc,
50  Interval<Treal> eigFInt,
51  Interval<Treal> hoF,
52  Interval<Treal> luF,
53  Treal toleratedEigenvalError,
54  Treal toleratedSubspaceError,
55  normType normForTruncation_,
56  int maxS = 100)
57  : n(nn),nocc(noc), step(new PuriStepInfo<Treal, Tvector, TdebugPolicy>[maxS]),
59  eigFInterval(eigFInt), homoF(hoF), lumoF(luF),
60  tolSubspaceError(toleratedSubspaceError),
61  tolEigenvalError(toleratedEigenvalError),
62  normForTruncation(normForTruncation_) {
63  for (int ind = 0; ind < maxSteps; ++ind)
64  step[ind] =
66  }
67  virtual ~PuriInfo() {
68  delete[] step;
69  }
80  void improveInfo();
82  return eigFInterval;
83  }
85  nSteps++;
87  return step[nSteps - 1];
88  }
90  assert(ind >= 0);
91  assert(ind < nSteps);
92  return step[ind];
93  }
94 
95 
97  inline Treal subspaceError() const {
98  return subspaceError(nSteps);
99  }
105  void estimateStepsLeft(int& stepsLeft, Treal& thresh) const;
106  Treal getOptimalThresh() const;
107  Treal getThreshIncreasingGap(Interval<Treal> const & middleGap) const;
108  bool ShouldComputeXmX2EuclNormAccurately(Treal & howAccurate) const;
112  inline Interval<Treal> getHomoF() const {return homoF;}
116  inline Interval<Treal> getLumoF() const {return lumoF;}
117  inline int getMaxSteps() const {return maxSteps;}
118  inline int getNSteps() const {return nSteps;}
119  /* Returns a vector or 0 and 1 giving the used poly sequence. */
120  void getPolys(std::vector<int> & resultVector) {
121  resultVector.resize(nSteps);
122  for(int i = 0; i < nSteps; i++)
123  resultVector[i] = step[i].getPoly();
124  }
125  /* Returns a vector of chosen threshold values. */
126  void getThreshValues(std::vector<Treal> & resultVector) {
127  resultVector.resize(nSteps);
128  for(int i = 0; i < nSteps; i++)
129  resultVector[i] = step[i].getChosenThresh();
130  }
131  /* Tries to improve the homoF and lumoF eigenvalues.
132  * Called after the purification process.
133  */
136  void improveHomoLumoF();
137 
138  inline bool converged() {return step[nSteps - 1].converged();}
139 
140  double getAccumulatedTimeSquare() const;
141  double getAccumulatedTimeThresh() const;
142  double getAccumulatedTimeXmX2Norm() const;
143  double getAccumulatedTimeTotal() const;
144 
145  void mTimings(std::ostream & file) const;
146  void mInfo(std::ostream & file) const;
147  void mMemUsage(std::ostream & file) const;
152  bool homoIsAccuratelyKnown() const {
153  bool res = false;
154  for(int ind = 0; ind < nSteps; ++ind)
155  res = res || step[ind].homoIsAccuratelyKnown(tolSubspaceError / 100);
156  return res;
157  }
162  bool lumoIsAccuratelyKnown() const {
163  bool res = false;
164  for(int ind = 0; ind < nSteps; ++ind)
165  res = res || step[ind].lumoIsAccuratelyKnown(tolSubspaceError / 100);
166  return res;
167  }
168 
170  return normForTruncation;
171  }
172 
173  void getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
174  Tvector & eigVecHOMO) const;
175 
176  protected:
177  int n;
178  int nocc;
181  int maxSteps;
183  int nSteps;
203  Treal const tolSubspaceError;
207  Treal const tolEigenvalError;
214  Treal subspaceError(int end) const;
215 
216  private:
217  };
218 
219 
220  template<typename Treal, typename Tvector, typename TdebugPolicy>
222  if ( step[nSteps-1].getCorrectOccupation() )
223  return;
224  step[nSteps-1].setCorrectOccupation();
226  return;
227  }
228 
229  template<typename Treal, typename Tvector, typename TdebugPolicy>
231  if (step[nSteps-1].getCorrectOccupation()) {
232  Treal distance;
233  Interval<Treal> middleInt;
234  Interval<Treal> zeroOneInt(0.0,1.0);
235  for (int ind = 0; ind < nSteps; ind++) {
236  Treal XmX2Eucl = step[ind].getXmX2EuclNorm().upp();
237  if ( XmX2Eucl < 1 / (Treal)4) {
238  distance = (1 - template_blas_sqrt(1 - 4 * XmX2Eucl)) / 2;
239  middleInt = Interval<Treal>(distance, 1.0 - distance);
240  int i = ind;
241  while (!middleInt.empty() && i < nSteps - 1) {
242  middleInt.puriStep(step[i].getPoly());
243  middleInt.decrease(step[i+1].getActualThresh());
244  /* Make sure we stay in [0, 1] */
245  middleInt.intersect(zeroOneInt);
246  ++i;
247  } /* end while */
248  if (middleInt.cover(0.5))
249  step[ind].setCorrectOccupation();
250  } /* end if */
251  } /* end for */
252  }
253  }
254 
255  template<typename Treal, typename Tvector, typename TdebugPolicy>
257  for (int ind = nSteps - 2; ind >= 0; ind--) {
258  step[ind].exchangeInfoWithNext(step[ind + 1]);
259  }
260  for (int ind = 0; ind < nSteps - 1; ind++) {
261  step[ind].exchangeInfoWithNext(step[ind + 1]);
262  }
263  }
264 
265  template<typename Treal, typename Tvector, typename TdebugPolicy>
267  Treal error = 0;
268  for (int ind = 0; ind < end; ind++) {
269  error += step[ind].subspaceError();
270  }
271  return error;
272  }
273 
274  template<typename Treal, typename Tvector, typename TdebugPolicy>
276  estimateStepsLeft(int& stepsLeft, Treal& thresh) const {
277  stepsLeft = -1;
278  thresh = 0;
279  Interval<Treal> initialGap;
280  Interval<Treal> gap;
281  Treal tolError = tolSubspaceError - subspaceError(nSteps - 1);
282  Treal initialAltThresh = 1;
283  if (tolError <= 0.0)
284  return;
285 
286  /* Compute number of steps needed to converge. */
287  /* Compute initial interval */
288  /* nSteps == 0 means that a Purification object has not been
289  * associated to this instance yet.
290  * nSteps == 1 means that a Purification object has been associated
291  * to this instance but no steps have been taken yet.
292  */
293  Treal lastThresh = 0;
294  if (nSteps == 0 || nSteps == 1) {
295  Treal lmax = eigFInterval.upp();
296  Treal lmin = eigFInterval.low();
297  /* Compute transformed homo and lumo eigenvalues. */
298  initialGap = Interval<Treal>((lumoF.low() - lmax) / (lmin - lmax),
299  (homoF.upp() - lmax) / (lmin - lmax));
300  }
301  else {
302  initialGap = Interval<Treal>(step[nSteps - 2].getLumo().upp(),
303  step[nSteps - 2].getHomo().low());
304  initialAltThresh = getThreshIncreasingGap(initialGap);
305  lastThresh = step[nSteps - 2].getChosenThresh();
306  initialAltThresh = initialAltThresh > lastThresh / 10 ?
307  initialAltThresh : lastThresh / 10;
308  initialGap.puriStep(step[nSteps - 2].getPoly());
309  }
310  if (initialGap.empty())
311  return;
312 #if 0
313  /* Already converged? */
314  if (1.0 - initialGap.upp() < tolEigenvalError &&
315  initialGap.low() < tolEigenvalError) {
316  stepsLeft = 0;
317  thresh = 0;
318  return;
319  }
320 #endif
321  Treal thetaPerStep = 0.0; /* Tolerated subspace error per step. */
322  Treal altThresh = 0; /* Maximum threshold that guarantees increased
323  * gap.
324  */
325  Treal currThresh = 0; /* Chosen threshold. */
326  int stepsLeftPrev = -2;
327 
328  while (stepsLeft > stepsLeftPrev) {
329  gap = initialGap;
330  altThresh = initialAltThresh;
331  currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
332  currThresh = currThresh < altThresh ? currThresh : altThresh;
333  gap.decrease(currThresh);
334  lastThresh = currThresh;
335  stepsLeftPrev = stepsLeft;
336  stepsLeft = 0;
337  while (!gap.empty() &&
338  (1.0 - gap.upp() > tolEigenvalError ||
339  gap.low() > tolEigenvalError)) {
340  altThresh = getThreshIncreasingGap(gap);
341  altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
342 
343  gap.puriStep(gap.upp() + gap.low() < 1);
344 
345  currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
346  currThresh = currThresh < altThresh ? currThresh : altThresh;
347  gap.decrease(currThresh);
348  lastThresh = currThresh;
349  ++stepsLeft;
350  }
351  thetaPerStep = tolError / (stepsLeft + 1);
352  }
353 
354  /* Compute optimal threshold for convergence of subspace. */
355  thetaPerStep = tolError / (stepsLeftPrev + 1);
356  thresh = thetaPerStep * initialGap.length() / (1 + thetaPerStep);
357  return;
358  }
359 
360 
361  /* FIXME !! */
362  template<typename Treal, typename Tvector, typename TdebugPolicy>
364  int stepsLeft = -1;
365  Treal thresh = 0.0;
366  estimateStepsLeft(stepsLeft, thresh);
367  if (nSteps > 0)
368  step[nSteps - 1].setEstimatedStepsLeft(stepsLeft);
369  if (stepsLeft < 0)
370  thresh = tolSubspaceError / 100; /* No reason */
371  if (nSteps > 1) {
372  Interval<Treal> middleGap =
373  Interval<Treal>(step[nSteps - 2].getLumo().upp(),
374  step[nSteps - 2].getHomo().low());
375  if (!middleGap.empty()) {
376  /* Get thresh that definitely gives increasing gap. */
377  Treal altThresh = getThreshIncreasingGap(middleGap);
378  /* Allow thresh to decrease only to a ten times smaller threshold */
379  Treal lastThresh = step[nSteps - 2].getChosenThresh();
380  altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
381  thresh = thresh < altThresh ? thresh : altThresh;
382 #if 0
383  Interval<Treal> xmX2Eucl = step[nSteps - 2].getXmX2EuclNorm();
384  Treal tmp = 1 - xmX2Eucl.upp() * 4.0;
385  if (tmp > 0) {
386  Treal altThresh = 0.25 * (1 - template_blas_sqrt(tmp)) / 2;
387  thresh = thresh < altThresh ? thresh : altThresh;
388  }
389 #endif
390  }
391  }
392  //std::cout<<"nsteps left, thresh: "<<stepsLeft<<" , "<<thresh<<std::endl;
393  ASSERTALWAYS(thresh >= 0.0);
394  return thresh;
395  }
396 
397  template<typename Treal, typename Tvector, typename TdebugPolicy>
399  getThreshIncreasingGap(Interval<Treal> const & middleGap) const {
400  Treal x = middleGap.midPoint();
401  Treal delta = middleGap.upp() - x;
402  Treal thresh;
403 #if 0
404  thresh = delta * template_blas_fabs(2 * x - 1) / 10;
405 #else
406  /* After a purification step, truncation is done.
407  * The threshold t (measured by the Euclidean norm of the error matrix)
408  * has to satisfy t <= | 1 - 2x | * d / 2 where x and d are the midpoint
409  * and the length of the HOMO-LUMO gap respectively. In this way the
410  * gap is guaranteed to increase to the next step. Note that x = 0.5 and
411  * d = 0 gives zero threshold. x -> 0.5 as the process converges.
412  */
413  if (delta > 0.4)
414  /* This choice gives much better estimation of the remaining
415  * number of iterations.
416  */
417  /* Close to convergence we chose quadratical dependence on the
418  * midpoint since we expect the midpoint to go to 0.5 quadratically
419  * at convergence.
420  */
421  thresh = delta * template_blas_fabs(2 * x - 1) * template_blas_fabs(2 * x - 1);
422  else
423  thresh = delta * template_blas_fabs(2 * x - 1) / 2;
424  // thresh = thresh > 1e-7 ? thresh : 1e-7;
425 #endif
426  /**************************************************************
427  ELIAS NOTE 2010-11-18: I got assertion failure when using gcc
428  in Fedora 14 [g++ (GCC) 4.5.1 20100924 (Red Hat 4.5.1-4)] and
429  it turned out to be because the fabs call returned zero while
430  thresh was something like 1e-30 for single precision. Therefore
431  I added std::numeric_limits<Treal>::epsilon() in the assertion,
432  which seemed to solve the problem.
433  ***************************************************************/
434  ASSERTALWAYS(thresh <= delta * template_blas_fabs(2 * x - 1) + std::numeric_limits<Treal>::epsilon());
435  ASSERTALWAYS(thresh >= 0.0);
436  return thresh;
437  }
438 
439  template<typename Treal, typename Tvector, typename TdebugPolicy>
441  ShouldComputeXmX2EuclNormAccurately(Treal & howAccurate) const {
442 
443  if (nSteps == 0 || nSteps == 1) {
444  howAccurate = 0;
445  return false;
446  }
447  Treal ep = 0.207106781186547;
448  /* = (sqrt(2) - 1) / 2
449  * This value is obtained by noting that x = 1 / sqrt(2) -> x^2 = 1 / 2.
450  * If an interval ]1 - 1 / sqrt(2), 1 / sqrt(2)[ is empty from eigenvalues,
451  * and if the occ. count is correct, then the occupation
452  * count is guaranteed to be correct after one step as well.
453  * This transforms to the value (sqrt(2) - 1) / 2 for the
454  * ||X - X^2||_2 norm via theorem 1.
455  */
456  Interval<Treal> homoCopy = step[nSteps - 2].getHomo();
457  Interval<Treal> lumoCopy = step[nSteps - 2].getLumo();
458  homoCopy.puriStep(step[nSteps - 2].getPoly());
459  lumoCopy.puriStep(step[nSteps - 2].getPoly());
460  /* Note: we changed this from getActualThresh() to
461  getChosenThresh() to avoid ridiculously small values for
462  howAccurate, as happened earlier for cases when no truncation
463  occurred, i.e. when the matrices were small. */
464  howAccurate = step[nSteps - 1].getChosenThresh() / 100;
465  Treal highestPossibleAccuracy = 10.0 * step[nSteps - 1].getEigAccLoss();
466  ASSERTALWAYS(howAccurate >= 0);
467  ASSERTALWAYS(highestPossibleAccuracy >= 0);
468  howAccurate = howAccurate > highestPossibleAccuracy ?
469  howAccurate : highestPossibleAccuracy;
470 
471  if (homoCopy.length() > 0.2 || lumoCopy.length() > 0.2) {
472  /* Base decision on n0 and n1 and XmX2Norm from previous step */
473  if (step[nSteps - 2].getN0() / (n - nocc) > 0.5 &&
474  step[nSteps - 2].getN1() / nocc > 0.5 &&
475  step[nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
476  return true;
477  else
478  return false;
479  }
480  else {
481  /* Decision can probably be made from homo and lumo estimates */
482  bool computeHomo = true; /* Do we want to try to compute homo */
483  bool computeLumo = true; /* Do we want to try to compute lumo */
484  if (homoIsAccuratelyKnown() ||
485  homoCopy.upp() < 0.5 ||
486  template_blas_fabs(lumoCopy.low() - 0.5) < template_blas_fabs(homoCopy.low() - 0.5))
487  computeHomo = false;
488  if (lumoIsAccuratelyKnown() ||
489  lumoCopy.low() > 0.5 ||
490  template_blas_fabs(homoCopy.upp() - 0.5) < template_blas_fabs(lumoCopy.upp() - 0.5))
491  computeLumo = false;
492  return computeHomo || computeLumo;
493  }
494  }
495 
496 
497  template<typename Treal, typename Tvector, typename TdebugPolicy>
500  Treal lmax = eigFInterval.upp();
501  Treal lmin = eigFInterval.low();
502  Interval<Treal> altHomo(step[0].getHomo() * (lmin - lmax) + lmax);
503  Interval<Treal> altLumo(step[0].getLumo() * (lmin - lmax) + lmax);
504  homoF.intersect(altHomo);
505  lumoF.intersect(altLumo);
506  }
507 
508  template<typename Treal, typename Tvector, typename TdebugPolicy>
510  double accTime = 0;
511  for (int ind = 0; ind < nSteps; ++ind)
512  accTime += (double)step[ind].getTimeSquare();
513  return accTime;
514  }
515  template<typename Treal, typename Tvector, typename TdebugPolicy>
517  double accTime = 0;
518  for (int ind = 0; ind < nSteps; ++ind)
519  accTime += (double)step[ind].getTimeThresh();
520  return accTime;
521  }
522  template<typename Treal, typename Tvector, typename TdebugPolicy>
524  double accTime = 0;
525  for (int ind = 0; ind < nSteps; ++ind)
526  accTime += (double)step[ind].getTimeXmX2Norm();
527  return accTime;
528  }
529  template<typename Treal, typename Tvector, typename TdebugPolicy>
531  double accTime = 0;
532  for (int ind = 0; ind < nSteps; ++ind)
533  accTime += (double)step[ind].getTimeTotal();
534  return accTime;
535  }
536 
537 
538  template<typename Treal, typename Tvector, typename TdebugPolicy>
540  mTimings(std::ostream & s) const {
541  s<<"puriTime = [";
542  for (int ind = 0; ind < nSteps; ++ind) {
543  s<<
544  step[ind].getTimeSquare()<<" "<<
545  step[ind].getTimeThresh()<<" "<<
546  step[ind].getTimeXmX2Norm()<<" "<<
547  step[ind].getTimeTotal()<<" "<<
548  step[ind].getTimeXX2Write()<<" "<<
549  step[ind].getTimeXX2Read()<<" "<<
550  std::endl;
551  }
552  s<<"];\n"<<"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
553  "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
554  " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
555  s<<"figure; plot(puriTime(:,4),'-x')"<<std::endl;
556  }
557 
558  template<typename Treal, typename Tvector, typename TdebugPolicy>
560  mInfo(std::ostream & s) const {
561  s<<"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
562  s<<"% Norm for truncation: "<< getNormTypeString(normForTruncation)
563  << std::endl;
564  s<<"n = "<<n<<";"<<std::endl;
565  s<<"nocc = "<<nocc<<";"<<std::endl;
566  s<<"tolSubspaceError = "<<tolSubspaceError<<";"<<std::endl;
567  s<<"tolEigenvalError = "<<tolEigenvalError<<";"<<std::endl;
568  s<<"correct_occupation_was_forced_flag = "
569  <<correct_occupation_was_forced_flag<<";"<<std::endl;
570  s<<"% Each row of the following matrix contains purification info \n"
571  <<"% for one step.\n"
572  <<"% The columns are arranged as: \n"
573  <<"% ind, n0, n1, trace(X), trace(X*X), poly, actualThresh, delta, "
574  <<"correctOcc, XmX2low, XmX2upp, HOMOmid, LUMOmid, "
575  <<"nnz(X), nnz(X*X), chosenThresh, estimatedStepsLeft "
576  <<std::endl;
577  s<<"puriMat = [";
578  for (int ind = 0; ind < nSteps; ++ind) {
579  s<<
580  ind+1<<" "<<
581  step[ind].getN0()<<" "<<
582  step[ind].getN1()<<" "<<
583  step[ind].getTraceX()<<" "<<
584  step[ind].getTraceX2()<<" "<<
585  step[ind].getPoly()<<" "<<
586  step[ind].getActualThresh()<<" "<<
587  step[ind].getDelta()<<" "<<
588  step[ind].getCorrectOccupation()<<" "<<
589  step[ind].getXmX2EuclNorm().low()<<" "<<
590  step[ind].getXmX2EuclNorm().upp()<<" "<<
591  step[ind].getHomo().midPoint()<<" "<<
592  step[ind].getLumo().midPoint()<<" "<<
593  step[ind].getNnzX()<<" "<<
594  step[ind].getNnzX2()<<" "<<
595  step[ind].getChosenThresh()<<" "<<
596  step[ind].getEstimatedStepsLeft()<<" "<<
597  std::endl;
598  }
599  s<<"];\n";
600  s<<"if(1)\n";
601  s<<"ind = puriMat(:,1);\n";
602  s<<"figure; \n"
603  <<"plot(ind,puriMat(:,12),'x-b',ind,puriMat(:,13),'o-r')\n"
604  <<"legend('HOMO','LUMO'),xlabel('Iteration')\n"
605  <<"axis([0 max(ind) 0 1])\n";
606  s<<"figure; \n"
607  <<"plot(ind,100*puriMat(:,2)/(n-nocc),'o-r',...\n"
608  <<"ind, 100*puriMat(:,3)/nocc,'x-b')\n"
609  <<"legend('N0','N1'),xlabel('Iteration'), ylabel('Percentage')\n"
610  <<"axis([0 max(ind) 0 100])\n";
611  s<<"figure; \n"
612  <<"subplot(211)\n"
613  <<"plot(ind,100*puriMat(:,14)/(n*n),'o-r',...\n"
614  <<"ind, 100*puriMat(:,15)/(n*n),'x-b')\n"
615  <<"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
616  <<"ylabel('Percentage')\n"
617  <<"axis([0 max(ind) 0 100])\n";
618  s<<"subplot(212)\n"
619  <<"semilogy(ind,puriMat(:,16),'x-r',ind,puriMat(:,7),'o-b')\n"
620  <<"xlabel('Iteration'), ylabel('Threshold')\n"
621  <<"legend('Chosen threshold', 'Actual threshold')\n"
622  <<"axis([0 max(ind) min(puriMat(:,7))/10 max(puriMat(:,16))*10])\n";
623  s<<"figure; \n"
624  <<"plot(ind,ind(end:-1:1)-1,ind,puriMat(:,17),'x-r')\n"
625  <<"legend('Steps left', 'Estimated steps left')\n"
626  <<"xlabel('Iteration')\n";
627 
628  s<<"end\n";
629  }
630 
631  template<typename Treal, typename Tvector, typename TdebugPolicy>
633  mMemUsage(std::ostream & s) const {
634  s<<"puriMemUsage = [";
635  for (int ind = 0; ind < nSteps; ++ind) {
636  MemUsage::Values memUsageBeforeTrunc = step[ind].getMemUsageBeforeTrunc();
637  MemUsage::Values memUsageInXmX2Diff = step[ind].getMemUsageInXmX2Diff();
638  s<<
639  memUsageBeforeTrunc.res <<" "<<
640  memUsageBeforeTrunc.virt<<" "<<
641  memUsageBeforeTrunc.peak<<" "<<
642  memUsageInXmX2Diff.res <<" "<<
643  memUsageInXmX2Diff.virt<<" "<<
644  memUsageInXmX2Diff.peak<<" "<<
645  std::endl;
646  }
647  s<<"];\n";
648  s<<"figure; \n"
649  <<"plot(puriMemUsage(:,1),'x-b')\n"
650  << "hold on\n"
651  <<"plot(puriMemUsage(:,2),'o-r')\n"
652  <<"plot(puriMemUsage(:,3),'^-g')\n"
653  <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage before trunc [GB]')\n"
654  <<"%force y axis to start at 0\n"
655  <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
656  << std::endl;
657  s<<"figure; \n"
658  <<"plot(puriMemUsage(:,4),'x-b')\n"
659  << "hold on\n"
660  <<"plot(puriMemUsage(:,5),'o-r')\n"
661  <<"plot(puriMemUsage(:,6),'^-g')\n"
662  <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage in XmX2Diff [GB]')\n"
663  <<"%force y axis to start at 0\n"
664  <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
665  << std::endl;
666  }
667 
668  template<typename Treal, typename Tvector, typename TdebugPolicy>
670  getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
671  Tvector & eigVecHOMO) const {
672  bool haveHOMO = 0;
673  bool haveLUMO = 0;
674  for (int ind = 0; ind < nSteps; ++ind) {
675  if (!haveHOMO && step[ind].getHomoWasComputed()) {
676  eigVecHOMO = *step[ind].getEigVecPtr();
677  haveHOMO = 1;
678  }
679  if (!haveLUMO && step[ind].getLumoWasComputed()) {
680  eigVecLUMO = *step[ind].getEigVecPtr();
681  haveLUMO = 1;
682  }
683  }
684  }
685 
686 
687 } /* end namespace mat */
688 #undef __ID__
689 #endif
Treal upp() const
Definition: Interval.h:143
bool lumoIsAccuratelyKnown() const
LUMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:162
double getAccumulatedTimeXmX2Norm() const
Definition: PuriInfo.h:523
void getHOMOandLUMOeigVecs(Tvector &eigVecLUMO, Tvector &eigVecHOMO) const
Definition: PuriInfo.h:670
double getAccumulatedTimeSquare() const
Definition: PuriInfo.h:509
std::string getNormTypeString(normType nType)
Definition: matInclude.cc:55
bool homoIsAccuratelyKnown() const
HOMO estimation is considered to be accurate if the error is smaller than tolSubspaceError / 100 in s...
Definition: PuriInfo.h:152
Interval< Treal > homoF
Interval containing the HOMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:193
double getAccumulatedTimeTotal() const
Definition: PuriInfo.h:530
void decrease(Treal const value)
Definition: Interval.h:138
Treal low() const
Definition: Interval.h:142
float peak
Definition: matInclude.h:174
void improveCorrectOccupation()
Improves the correct occupation indicator Call AFTER convergence and ONLY if it is known that the occ...
Definition: PuriInfo.h:230
PuriStepInfo< Treal, Tvector, TdebugPolicy > * step
Definition: PuriInfo.h:180
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
static ergo_real distance(const ergo_real *a, const ergo_real *b)
Coomputes distance between two points, they do not need to be of the Vector3D type.
Definition: dft_common.cc:530
void getPolys(std::vector< int > &resultVector)
Definition: PuriInfo.h:120
PuriStepInfo< Treal, Tvector, TdebugPolicy > & operator()(int const ind)
Definition: PuriInfo.h:89
double getAccumulatedTimeThresh() const
Definition: PuriInfo.h:516
void estimateStepsLeft(int &stepsLeft, Treal &thresh) const
Estimates the number of steps (upper bound) needed for convergence based on homo/lumo information and...
Definition: PuriInfo.h:276
#define ASSERTALWAYS(x)
Definition: DebugPolicies.h:83
Definition: allocate.cc:30
int nocc
Number of occupied orbitals.
Definition: PuriInfo.h:178
Treal template_blas_fabs(Treal x)
Definition: Interval.h:44
void mTimings(std::ostream &file) const
Definition: PuriInfo.h:540
void puriStep(int poly)
Definition: Interval.h:228
Treal getOptimalThresh() const
Definition: PuriInfo.h:363
void forceCorrectOccupation()
Set the correctOccupation flag in the current step to 1.
Definition: PuriInfo.h:221
Contains information about the matrix before a purification step and about the step.
Definition: PuriStepInfo.h:55
float res
Definition: matInclude.h:172
Interval< Treal > getHomoF() const
Returns the best interval containing the homo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:112
void mMemUsage(std::ostream &file) const
Definition: PuriInfo.h:633
Contains information about a purification process.
Definition: PuriInfo.h:47
void mInfo(std::ostream &file) const
Definition: PuriInfo.h:560
float virt
Definition: matInclude.h:173
Treal midPoint() const
Definition: Interval.h:113
Treal getThreshIncreasingGap(Interval< Treal > const &middleGap) const
Definition: PuriInfo.h:399
Interval< Treal > const eigFInterval
Interval containing all eigenvalues before transformation to the [0, 1] interval. ...
Definition: PuriInfo.h:188
PuriStepInfo< Treal, Tvector, TdebugPolicy > & getNext()
Definition: PuriInfo.h:84
Treal subspaceError() const
Returns the subspace error introduced so far.
Definition: PuriInfo.h:97
bool correct_occupation_was_forced_flag
Correct occupation was assumed, not guaranteed.
Definition: PuriInfo.h:186
bool ShouldComputeXmX2EuclNormAccurately(Treal &howAccurate) const
Definition: PuriInfo.h:441
Interval< Treal > lumoF
Interval containing the LUMO eigenvalue before transformation to [0, 1].
Definition: PuriInfo.h:198
Interval< Treal > getLumoF() const
Returns the best interval containing the lumo eigenvalue (not transformed to [0, 1]) ...
Definition: PuriInfo.h:116
PuriStepInfo class.
bool cover(Treal const value) const
Definition: Interval.h:117
int nSteps
Number of taken steps.
Definition: PuriInfo.h:183
bool converged()
Definition: PuriInfo.h:138
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
int n
System size.
Definition: PuriInfo.h:177
bool empty() const
Definition: Interval.h:49
Treal const tolEigenvalError
Tolerated error in eigenvalues at convergence.
Definition: PuriInfo.h:207
Definition: matInclude.h:171
int maxSteps
Capacity of step array.
Definition: PuriInfo.h:181
bool correct_occupation_was_forced() const
Definition: PuriInfo.h:134
normType getNormForTruncation() const
Definition: PuriInfo.h:169
void improveHomoLumoF()
Definition: PuriInfo.h:499
Interval< Treal > getEigFInterval() const
Definition: PuriInfo.h:81
normType const normForTruncation
Norm used for truncation of small matrix elements.
Definition: PuriInfo.h:210
int getNSteps() const
Definition: PuriInfo.h:118
Treal const tolSubspaceError
Tolerated error in angle between correct and computed subspace.
Definition: PuriInfo.h:203
void improveInfo()
Improve homo / lumo values in each step.
Definition: PuriInfo.h:256
PuriInfo(int nn, int noc, Interval< Treal > eigFInt, Interval< Treal > hoF, Interval< Treal > luF, Treal toleratedEigenvalError, Treal toleratedSubspaceError, normType normForTruncation_, int maxS=100)
Definition: PuriInfo.h:49
void getThreshValues(std::vector< Treal > &resultVector)
Definition: PuriInfo.h:126
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:135
int getMaxSteps() const
Definition: PuriInfo.h:117
virtual ~PuriInfo()
Definition: PuriInfo.h:67