ergo
PuriInfo.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_PURIINFO
37 #define MAT_PURIINFO
38 #include <math.h>
39 #include <iomanip>
40 #include "PuriStepInfo.h"
41 #define __ID__ "$Id: PuriInfo.h 4437 2012-07-05 09:01:18Z elias $"
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  /* Tries to improve the homoF and lumoF eigenvalues.
120  * Called after the purification process.
121  */
124  void improveHomoLumoF();
125 
126  inline bool converged() {return step[nSteps - 1].converged();}
127 
128  double getAccumulatedTimeSquare() const;
129  double getAccumulatedTimeThresh() const;
130  double getAccumulatedTimeXmX2Norm() const;
131  double getAccumulatedTimeTotal() const;
132 
133  void mTimings(std::ostream & file) const;
134  void mInfo(std::ostream & file) const;
135  void mMemUsage(std::ostream & file) const;
140  bool homoIsAccuratelyKnown() const {
141  bool res = false;
142  for(int ind = 0; ind < nSteps; ++ind)
143  res = res || step[ind].homoIsAccuratelyKnown(tolSubspaceError / 100);
144  return res;
145  }
150  bool lumoIsAccuratelyKnown() const {
151  bool res = false;
152  for(int ind = 0; ind < nSteps; ++ind)
153  res = res || step[ind].lumoIsAccuratelyKnown(tolSubspaceError / 100);
154  return res;
155  }
156 
158  return normForTruncation;
159  }
160 
161  void getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
162  Tvector & eigVecHOMO) const;
163 
164  protected:
165  int n;
166  int nocc;
169  int maxSteps;
171  int nSteps;
191  Treal const tolSubspaceError;
195  Treal const tolEigenvalError;
202  Treal subspaceError(int end) const;
203 
204  private:
205  };
206 
207 
208  template<typename Treal, typename Tvector, typename TdebugPolicy>
210  if ( step[nSteps-1].getCorrectOccupation() )
211  return;
212  step[nSteps-1].setCorrectOccupation();
213  correct_occupation_was_forced_flag = true;
214  return;
215  }
216 
217  template<typename Treal, typename Tvector, typename TdebugPolicy>
219  if (step[nSteps-1].getCorrectOccupation()) {
220  Treal distance;
221  Interval<Treal> middleInt;
222  Interval<Treal> zeroOneInt(0.0,1.0);
223  for (int ind = 0; ind < nSteps; ind++) {
224  Treal XmX2Eucl = step[ind].getXmX2EuclNorm().upp();
225  if ( XmX2Eucl < 1 / (Treal)4) {
226  distance = (1 - template_blas_sqrt(1 - 4 * XmX2Eucl)) / 2;
227  middleInt = Interval<Treal>(distance, 1.0 - distance);
228  int i = ind;
229  while (!middleInt.empty() && i < nSteps - 1) {
230  middleInt.puriStep(step[i].getPoly());
231  middleInt.decrease(step[i+1].getActualThresh());
232  /* Make sure we stay in [0, 1] */
233  middleInt.intersect(zeroOneInt);
234  ++i;
235  } /* end while */
236  if (middleInt.cover(0.5))
237  step[ind].setCorrectOccupation();
238  } /* end if */
239  } /* end for */
240  }
241  }
242 
243  template<typename Treal, typename Tvector, typename TdebugPolicy>
245  for (int ind = nSteps - 2; ind >= 0; ind--) {
246  step[ind].exchangeInfoWithNext(step[ind + 1]);
247  }
248  for (int ind = 0; ind < nSteps - 1; ind++) {
249  step[ind].exchangeInfoWithNext(step[ind + 1]);
250  }
251  }
252 
253  template<typename Treal, typename Tvector, typename TdebugPolicy>
255  Treal error = 0;
256  for (int ind = 0; ind < end; ind++) {
257  error += step[ind].subspaceError();
258  }
259  return error;
260  }
261 
262  template<typename Treal, typename Tvector, typename TdebugPolicy>
264  estimateStepsLeft(int& stepsLeft, Treal& thresh) const {
265  stepsLeft = -1;
266  thresh = 0;
267  Interval<Treal> initialGap;
268  Interval<Treal> gap;
269  Treal tolError = tolSubspaceError - subspaceError(nSteps - 1);
270  Treal initialAltThresh = 1;
271  if (tolError <= 0.0)
272  return;
273 
274  /* Compute number of steps needed to converge. */
275  /* Compute initial interval */
276  /* nSteps == 0 means that a Purification object has not been
277  * associated to this instance yet.
278  * nSteps == 1 means that a Purification object has been associated
279  * to this instance but no steps have been taken yet.
280  */
281  Treal lastThresh = 0;
282  if (nSteps == 0 || nSteps == 1) {
283  Treal lmax = eigFInterval.upp();
284  Treal lmin = eigFInterval.low();
285  /* Compute transformed homo and lumo eigenvalues. */
286  initialGap = Interval<Treal>((lumoF.low() - lmax) / (lmin - lmax),
287  (homoF.upp() - lmax) / (lmin - lmax));
288  }
289  else {
290  initialGap = Interval<Treal>(step[nSteps - 2].getLumo().upp(),
291  step[nSteps - 2].getHomo().low());
292  initialAltThresh = getThreshIncreasingGap(initialGap);
293  lastThresh = step[nSteps - 2].getChosenThresh();
294  initialAltThresh = initialAltThresh > lastThresh / 10 ?
295  initialAltThresh : lastThresh / 10;
296  initialGap.puriStep(step[nSteps - 2].getPoly());
297  }
298  if (initialGap.empty())
299  return;
300 #if 0
301  /* Already converged? */
302  if (1.0 - initialGap.upp() < tolEigenvalError &&
303  initialGap.low() < tolEigenvalError) {
304  stepsLeft = 0;
305  thresh = 0;
306  return;
307  }
308 #endif
309  Treal thetaPerStep = 0.0; /* Tolerated subspace error per step. */
310  Treal altThresh = 0; /* Maximum threshold that guarantees increased
311  * gap.
312  */
313  Treal currThresh = 0; /* Chosen threshold. */
314  int stepsLeftPrev = -2;
315 
316  while (stepsLeft > stepsLeftPrev) {
317  gap = initialGap;
318  altThresh = initialAltThresh;
319  currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
320  currThresh = currThresh < altThresh ? currThresh : altThresh;
321  gap.decrease(currThresh);
322  lastThresh = currThresh;
323  stepsLeftPrev = stepsLeft;
324  stepsLeft = 0;
325  while (!gap.empty() &&
326  (1.0 - gap.upp() > tolEigenvalError ||
327  gap.low() > tolEigenvalError)) {
328  altThresh = getThreshIncreasingGap(gap);
329  altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
330 
331  gap.puriStep(gap.upp() + gap.low() < 1);
332 
333  currThresh = thetaPerStep * gap.length() / (1 + thetaPerStep);
334  currThresh = currThresh < altThresh ? currThresh : altThresh;
335  gap.decrease(currThresh);
336  lastThresh = currThresh;
337  ++stepsLeft;
338  }
339  thetaPerStep = tolError / (stepsLeft + 1);
340  }
341 
342  /* Compute optimal threshold for convergence of subspace. */
343  thetaPerStep = tolError / (stepsLeftPrev + 1);
344  thresh = thetaPerStep * initialGap.length() / (1 + thetaPerStep);
345  return;
346  }
347 
348 
349  /* FIXME !! */
350  template<typename Treal, typename Tvector, typename TdebugPolicy>
352  int stepsLeft = -1;
353  Treal thresh = 0.0;
354  estimateStepsLeft(stepsLeft, thresh);
355  if (nSteps > 0)
356  step[nSteps - 1].setEstimatedStepsLeft(stepsLeft);
357  if (stepsLeft < 0)
358  thresh = tolSubspaceError / 100; /* No reason */
359  if (nSteps > 1) {
360  Interval<Treal> middleGap =
361  Interval<Treal>(step[nSteps - 2].getLumo().upp(),
362  step[nSteps - 2].getHomo().low());
363  if (!middleGap.empty()) {
364  /* Get thresh that definitely gives increasing gap. */
365  Treal altThresh = getThreshIncreasingGap(middleGap);
366  /* Allow thresh to decrease only to a ten times smaller threshold */
367  Treal lastThresh = step[nSteps - 2].getChosenThresh();
368  altThresh = altThresh > lastThresh / 10 ? altThresh : lastThresh / 10;
369  thresh = thresh < altThresh ? thresh : altThresh;
370 #if 0
371  Interval<Treal> xmX2Eucl = step[nSteps - 2].getXmX2EuclNorm();
372  Treal tmp = 1 - xmX2Eucl.upp() * 4.0;
373  if (tmp > 0) {
374  Treal altThresh = 0.25 * (1 - template_blas_sqrt(tmp)) / 2;
375  thresh = thresh < altThresh ? thresh : altThresh;
376  }
377 #endif
378  }
379  }
380  //std::cout<<"nsteps left, thresh: "<<stepsLeft<<" , "<<thresh<<std::endl;
381  ASSERTALWAYS(thresh >= 0.0);
382  return thresh;
383  }
384 
385  template<typename Treal, typename Tvector, typename TdebugPolicy>
387  getThreshIncreasingGap(Interval<Treal> const & middleGap) const {
388  Treal x = middleGap.midPoint();
389  Treal delta = middleGap.upp() - x;
390  Treal thresh;
391 #if 0
392  thresh = delta * template_blas_fabs(2 * x - 1) / 10;
393 #else
394  /* After a purification step, truncation is done.
395  * The threshold t (measured by the Euclidean norm of the error matrix)
396  * has to satisfy t <= | 1 - 2x | * d / 2 where x and d are the midpoint
397  * and the length of the HOMO-LUMO gap respectively. In this way the
398  * gap is guaranteed to increase to the next step. Note that x = 0.5 and
399  * d = 0 gives zero threshold. x -> 0.5 as the process converges.
400  */
401  if (delta > 0.4)
402  /* This choice gives much better estimation of the remaining
403  * number of iterations.
404  */
405  /* Close to convergence we chose quadratical dependence on the
406  * midpoint since we expect the midpoint to go to 0.5 quadratically
407  * at convergence.
408  */
409  thresh = delta * template_blas_fabs(2 * x - 1) * template_blas_fabs(2 * x - 1);
410  else
411  thresh = delta * template_blas_fabs(2 * x - 1) / 2;
412  // thresh = thresh > 1e-7 ? thresh : 1e-7;
413 #endif
414  /**************************************************************
415  ELIAS NOTE 2010-11-18: I got assertion failure when using gcc
416  in Fedora 14 [g++ (GCC) 4.5.1 20100924 (Red Hat 4.5.1-4)] and
417  it turned out to be because the fabs call returned zero while
418  thresh was something like 1e-30 for single precision. Therefore
419  I added std::numeric_limits<Treal>::epsilon() in the assertion,
420  which seemed to solve the problem.
421  ***************************************************************/
422  ASSERTALWAYS(thresh <= delta * template_blas_fabs(2 * x - 1) + std::numeric_limits<Treal>::epsilon());
423  ASSERTALWAYS(thresh >= 0.0);
424  return thresh;
425  }
426 
427  template<typename Treal, typename Tvector, typename TdebugPolicy>
429  ShouldComputeXmX2EuclNormAccurately(Treal & howAccurate) const {
430 
431  if (nSteps == 0 || nSteps == 1) {
432  howAccurate = 0;
433  return false;
434  }
435  Treal ep = 0.207106781186547;
436  /* = (sqrt(2) - 1) / 2
437  * This value is obtained by noting that x = 1 / sqrt(2) -> x^2 = 1 / 2.
438  * If an interval ]1 - 1 / sqrt(2), 1 / sqrt(2)[ is empty from eigenvalues,
439  * and if the occ. count is correct, then the occupation
440  * count is guaranteed to be correct after one step as well.
441  * This transforms to the value (sqrt(2) - 1) / 2 for the
442  * ||X - X^2||_2 norm via theorem 1.
443  */
444  Interval<Treal> homoCopy = step[nSteps - 2].getHomo();
445  Interval<Treal> lumoCopy = step[nSteps - 2].getLumo();
446  homoCopy.puriStep(step[nSteps - 2].getPoly());
447  lumoCopy.puriStep(step[nSteps - 2].getPoly());
448  /* Note: we changed this from getActualThresh() to
449  getChosenThresh() to avoid ridiculously small values for
450  howAccurate, as happened earlier for cases when no truncation
451  occurred, i.e. when the matrices were small. */
452  howAccurate = step[nSteps - 1].getChosenThresh() / 100;
453  Treal highestPossibleAccuracy = 10.0 * step[nSteps - 1].getEigAccLoss();
454  ASSERTALWAYS(howAccurate >= 0);
455  ASSERTALWAYS(highestPossibleAccuracy >= 0);
456  howAccurate = howAccurate > highestPossibleAccuracy ?
457  howAccurate : highestPossibleAccuracy;
458 
459  if (homoCopy.length() > 0.2 || lumoCopy.length() > 0.2) {
460  /* Base decision on n0 and n1 and XmX2Norm from previous step */
461  if (step[nSteps - 2].getN0() / (n - nocc) > 0.5 &&
462  step[nSteps - 2].getN1() / nocc > 0.5 &&
463  step[nSteps - 2].getXmX2EuclNorm().midPoint() < ep)
464  return true;
465  else
466  return false;
467  }
468  else {
469  /* Decision can probably be made from homo and lumo estimates */
470  bool computeHomo = true; /* Do we want to try to compute homo */
471  bool computeLumo = true; /* Do we want to try to compute lumo */
472  if (homoIsAccuratelyKnown() ||
473  homoCopy.upp() < 0.5 ||
474  template_blas_fabs(lumoCopy.low() - 0.5) < template_blas_fabs(homoCopy.low() - 0.5))
475  computeHomo = false;
476  if (lumoIsAccuratelyKnown() ||
477  lumoCopy.low() > 0.5 ||
478  template_blas_fabs(homoCopy.upp() - 0.5) < template_blas_fabs(lumoCopy.upp() - 0.5))
479  computeLumo = false;
480  return computeHomo || computeLumo;
481  }
482  }
483 
484 
485  template<typename Treal, typename Tvector, typename TdebugPolicy>
488  Treal lmax = eigFInterval.upp();
489  Treal lmin = eigFInterval.low();
490  Interval<Treal> altHomo(step[0].getHomo() * (lmin - lmax) + lmax);
491  Interval<Treal> altLumo(step[0].getLumo() * (lmin - lmax) + lmax);
492  homoF.intersect(altHomo);
493  lumoF.intersect(altLumo);
494  }
495 
496  template<typename Treal, typename Tvector, typename TdebugPolicy>
498  double accTime = 0;
499  for (int ind = 0; ind < nSteps; ++ind)
500  accTime += (double)step[ind].getTimeSquare();
501  return accTime;
502  }
503  template<typename Treal, typename Tvector, typename TdebugPolicy>
505  double accTime = 0;
506  for (int ind = 0; ind < nSteps; ++ind)
507  accTime += (double)step[ind].getTimeThresh();
508  return accTime;
509  }
510  template<typename Treal, typename Tvector, typename TdebugPolicy>
512  double accTime = 0;
513  for (int ind = 0; ind < nSteps; ++ind)
514  accTime += (double)step[ind].getTimeXmX2Norm();
515  return accTime;
516  }
517  template<typename Treal, typename Tvector, typename TdebugPolicy>
519  double accTime = 0;
520  for (int ind = 0; ind < nSteps; ++ind)
521  accTime += (double)step[ind].getTimeTotal();
522  return accTime;
523  }
524 
525 
526  template<typename Treal, typename Tvector, typename TdebugPolicy>
528  mTimings(std::ostream & s) const {
529  s<<"puriTime = [";
530  for (int ind = 0; ind < nSteps; ++ind) {
531  s<<
532  step[ind].getTimeSquare()<<" "<<
533  step[ind].getTimeThresh()<<" "<<
534  step[ind].getTimeXmX2Norm()<<" "<<
535  step[ind].getTimeTotal()<<" "<<
536  step[ind].getTimeXX2Write()<<" "<<
537  step[ind].getTimeXX2Read()<<" "<<
538  std::endl;
539  }
540  s<<"];\n"<<"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
541  "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
542  " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
543  s<<"figure; plot(puriTime(:,4),'-x')"<<std::endl;
544  }
545 
546  template<typename Treal, typename Tvector, typename TdebugPolicy>
548  mInfo(std::ostream & s) const {
549  s<<"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
550  s<<"% Norm for truncation: "<< getNormTypeString(normForTruncation)
551  << std::endl;
552  s<<"n = "<<n<<";"<<std::endl;
553  s<<"nocc = "<<nocc<<";"<<std::endl;
554  s<<"tolSubspaceError = "<<tolSubspaceError<<";"<<std::endl;
555  s<<"tolEigenvalError = "<<tolEigenvalError<<";"<<std::endl;
556  s<<"correct_occupation_was_forced_flag = "
557  <<correct_occupation_was_forced_flag<<";"<<std::endl;
558  s<<"% Each row of the following matrix contains purification info \n"
559  <<"% for one step.\n"
560  <<"% The columns are arranged as: \n"
561  <<"% ind, n0, n1, trace(X), trace(X*X), poly, actualThresh, delta, "
562  <<"correctOcc, XmX2low, XmX2upp, HOMOmid, LUMOmid, "
563  <<"nnz(X), nnz(X*X), chosenThresh, estimatedStepsLeft "
564  <<std::endl;
565  s<<"puriMat = [";
566  for (int ind = 0; ind < nSteps; ++ind) {
567  s<<
568  ind+1<<" "<<
569  step[ind].getN0()<<" "<<
570  step[ind].getN1()<<" "<<
571  step[ind].getTraceX()<<" "<<
572  step[ind].getTraceX2()<<" "<<
573  step[ind].getPoly()<<" "<<
574  step[ind].getActualThresh()<<" "<<
575  step[ind].getDelta()<<" "<<
576  step[ind].getCorrectOccupation()<<" "<<
577  step[ind].getXmX2EuclNorm().low()<<" "<<
578  step[ind].getXmX2EuclNorm().upp()<<" "<<
579  step[ind].getHomo().midPoint()<<" "<<
580  step[ind].getLumo().midPoint()<<" "<<
581  step[ind].getNnzX()<<" "<<
582  step[ind].getNnzX2()<<" "<<
583  step[ind].getChosenThresh()<<" "<<
584  step[ind].getEstimatedStepsLeft()<<" "<<
585  std::endl;
586  }
587  s<<"];\n";
588  s<<"if(1)\n";
589  s<<"ind = puriMat(:,1);\n";
590  s<<"figure; \n"
591  <<"plot(ind,puriMat(:,12),'x-b',ind,puriMat(:,13),'o-r')\n"
592  <<"legend('HOMO','LUMO'),xlabel('Iteration')\n"
593  <<"axis([0 max(ind) 0 1])\n";
594  s<<"figure; \n"
595  <<"plot(ind,100*puriMat(:,2)/(n-nocc),'o-r',...\n"
596  <<"ind, 100*puriMat(:,3)/nocc,'x-b')\n"
597  <<"legend('N0','N1'),xlabel('Iteration'), ylabel('Percentage')\n"
598  <<"axis([0 max(ind) 0 100])\n";
599  s<<"figure; \n"
600  <<"subplot(211)\n"
601  <<"plot(ind,100*puriMat(:,14)/(n*n),'o-r',...\n"
602  <<"ind, 100*puriMat(:,15)/(n*n),'x-b')\n"
603  <<"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
604  <<"ylabel('Percentage')\n"
605  <<"axis([0 max(ind) 0 100])\n";
606  s<<"subplot(212)\n"
607  <<"semilogy(ind,puriMat(:,16),'x-r',ind,puriMat(:,7),'o-b')\n"
608  <<"xlabel('Iteration'), ylabel('Threshold')\n"
609  <<"legend('Chosen threshold', 'Actual threshold')\n"
610  <<"axis([0 max(ind) min(puriMat(:,7))/10 max(puriMat(:,16))*10])\n";
611  s<<"figure; \n"
612  <<"plot(ind,ind(end:-1:1)-1,ind,puriMat(:,17),'x-r')\n"
613  <<"legend('Steps left', 'Estimated steps left')\n"
614  <<"xlabel('Iteration')\n";
615 
616  s<<"end\n";
617  }
618 
619  template<typename Treal, typename Tvector, typename TdebugPolicy>
621  mMemUsage(std::ostream & s) const {
622  s<<"puriMemUsage = [";
623  for (int ind = 0; ind < nSteps; ++ind) {
624  MemUsage::Values memUsageBeforeTrunc = step[ind].getMemUsageBeforeTrunc();
625  MemUsage::Values memUsageInXmX2Diff = step[ind].getMemUsageInXmX2Diff();
626  s<<
627  memUsageBeforeTrunc.res <<" "<<
628  memUsageBeforeTrunc.virt<<" "<<
629  memUsageBeforeTrunc.peak<<" "<<
630  memUsageInXmX2Diff.res <<" "<<
631  memUsageInXmX2Diff.virt<<" "<<
632  memUsageInXmX2Diff.peak<<" "<<
633  std::endl;
634  }
635  s<<"];\n";
636  s<<"figure; \n"
637  <<"plot(puriMemUsage(:,1),'x-b')\n"
638  << "hold on\n"
639  <<"plot(puriMemUsage(:,2),'o-r')\n"
640  <<"plot(puriMemUsage(:,3),'^-g')\n"
641  <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage before trunc [GB]')\n"
642  <<"%force y axis to start at 0\n"
643  <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
644  << std::endl;
645  s<<"figure; \n"
646  <<"plot(puriMemUsage(:,4),'x-b')\n"
647  << "hold on\n"
648  <<"plot(puriMemUsage(:,5),'o-r')\n"
649  <<"plot(puriMemUsage(:,6),'^-g')\n"
650  <<"legend('Resident','Virtual','Peak'),xlabel('Iteration'),ylabel('Mem usage in XmX2Diff [GB]')\n"
651  <<"%force y axis to start at 0\n"
652  <<"axissaved = axis; axissaved(3) = 0; axis(axissaved);\n"
653  << std::endl;
654  }
655 
656  template<typename Treal, typename Tvector, typename TdebugPolicy>
658  getHOMOandLUMOeigVecs(Tvector & eigVecLUMO,
659  Tvector & eigVecHOMO) const {
660  bool haveHOMO = 0;
661  bool haveLUMO = 0;
662  for (int ind = 0; ind < nSteps; ++ind) {
663  if (!haveHOMO && step[ind].getHomoWasComputed()) {
664  eigVecHOMO = *step[ind].getEigVecPtr();
665  haveHOMO = 1;
666  }
667  if (!haveLUMO && step[ind].getLumoWasComputed()) {
668  eigVecLUMO = *step[ind].getEigVecPtr();
669  haveLUMO = 1;
670  }
671  }
672  }
673 
674 
675 } /* end namespace mat */
676 #undef __ID__
677 #endif