ergo
Purification_scaled.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 PUR_PURIFICATION_SCALED
37 #define PUR_PURIFICATION_SCALED
38 #include <iomanip>
39 #include <sys/time.h>
40 #include <stdexcept>
41 #include "matrix_proxy.h"
42 #include "matInclude.h"
43 #include "Step.h"
44 namespace pur {
45  template<typename Tmatrix>
47  public:
48  typedef typename Tmatrix::real real;
49 
50  Purification_scaled(Tmatrix & F_and_D,
55  mat::Interval<real> const & hoF,
57  mat::Interval<real> const & luF,
59  real const toleratedEigenvalError,
62  real const toleratedSubspaceError,
65  int const max_steps,
66  mat::normType normForTruncation,
68  bool const use_scaling
70  );
71 
72 
73  void purify();
74 
76  mat::Interval<real> & luF );
81  void mInfo(std::ostream & file) const;
82 
83  void mTime(std::ostream & file) const;
84 
85  class Time;
86  private:
87 
88  static int get_poly( real const homo, real const lumo ) {
89  // return fabs(homo*homo + lumo*lumo - 1) >
90  // fabs(2*homo - homo*homo + 2*lumo - lumo*lumo - 1);
91  return homo + lumo < 1;
92  }
93 
94  static int estimated_steps_left( mat::Interval<real> eig_homo,
95  mat::Interval<real> eig_lumo,
96  real const toleratedEigenvalError,
97  int const max_steps,
98  bool const use_scaling );
106  static real get_threshold( real const gap,
107  int const steps_left,
108  real const acc_error,
109  real const toleratedSubspaceError );
119  mat::Interval<real>& lumo,
120  bool & homo_was_computed,
121  bool & lumo_was_computed );
129  real accumulated_error();
145  bool converged();
146 
147  Tmatrix & X;
148  Tmatrix X2;
149  int n;
150  std::vector<Step<real> > puri_steps;
152  real const tolEigenvalError;
154  real const tolSubspaceError;
156  unsigned int current_step;
161  };
162 
163  template<typename Tmatrix>
165  Purification_scaled( Tmatrix & F_and_D,
167  mat::Interval<real> const & hoF,
168  mat::Interval<real> const & luF,
169  real const toleratedEigenvalError,
170  real const toleratedSubspaceError,
171  int const max_steps,
172  mat::normType normForTruncation,
173  bool const use_scaling)
174  : X(F_and_D),
175  puri_steps(max_steps+1),
176  eigFInt(eigFInt),
177  tolEigenvalError(toleratedEigenvalError),
178  tolSubspaceError(toleratedSubspaceError),
179  current_step(0),
180  normTruncation(normForTruncation),
181  use_scaling(use_scaling),
182  homo_is_computed(0),
183  lumo_is_computed(0) {
184  Time t_tot; t_tot.tic();
185  n = X.get_nrows();
186  using namespace mat;
187  real lmin = eigFInt.low();
188  real lmax = eigFInt.upp();
189  X.add_identity(-lmax); /* Scale to [0, 1] interval and negate */
190  X *= ((real)1.0 / (lmin - lmax));
191  /* Compute transformed homo and lumo eigenvalues. */
192  mat::Interval<real> homo = hoF;
193  mat::Interval<real> lumo = luF;
194  // homo and lumo must be in the [lmin, lmax] interval
195  homo.intersect( mat::Interval<real>(lmin,lmax) );
196  lumo.intersect( mat::Interval<real>(lmin,lmax) );
197  if ( homo.empty() )
198  throw "homo interval empty in Purification_scaled constructor due to incorrect input.";
199  if ( lumo.empty() )
200  throw "lumo interval empty in Purification_scaled constructor due to incorrect input.";
201  homo = (homo - lmax) / (lmin - lmax);
202  lumo = (lumo - lmax) / (lmin - lmax);
203 
204  int steps_left = estimated_steps_left( homo, lumo,
206  2*puri_steps.size(),
207  use_scaling );
208  real acc_error = accumulated_error();
209  real chosen_thresh = get_threshold( homo.low() - lumo.upp(),
210  steps_left, acc_error,
212  Time t_thr; t_thr.tic();
213  real actual_thresh = X.thresh( chosen_thresh, normTruncation );
214  double wall_sec_thresh = t_thr.toc();
215  // Eigenvalues may have moved due to truncation. Therefore, we
216  // increase the intervals containing the HOMO and LUMO
217  // eigenvalues according to Weyl's theorem.
218  homo.increase( actual_thresh );
219  lumo.increase( actual_thresh );
220  // We do not believe that homo or lumo eigenvalues are outside
221  // [0, 1] here
222  mat::Interval<real> zero_one_int(0.0,1.0);
223  homo.intersect(zero_one_int);
224  lumo.intersect(zero_one_int);
225  real const ONE = 1.0;
226  Time t_sqr; t_sqr.tic();
227  X2 = ONE * X * X;
228  double wall_sec_square = t_sqr.toc();
229  double wall_sec_total = t_tot.toc();
231  homo, lumo,
232  X.trace(), X2.trace(),
233  chosen_thresh,
234  actual_thresh,
235  X.nnz(), X2.nnz(),
236  wall_sec_thresh,
237  wall_sec_square,
238  0,
239  wall_sec_total);
240  }
241 
242  template<typename Tmatrix>
244  using namespace mat;
245  while ( !converged() && current_step < puri_steps.size()-1 ) {
246  Time t_tot; t_tot.tic();
249 
250  int poly = get_poly( homo.low(), lumo.upp() );
251  // Scaling
252  real alpha = 1; // alpha == 1 means no scaling
253  if (use_scaling) {
254  if (poly) {
255  alpha = 2 / ( 1+homo.upp() );
256  }
257  else {
258  alpha = 2 / ( 2-lumo.low() );
259  }
260  }
261  if (poly) {
262  /* Polynomial 2 * alpha * x - alpha^2 * x^2 */
263  X2 *= (real)-(alpha*alpha);
264  X2 += ((real)2.0 * alpha) * X;
265  }
266  else {
267  /* Polynomial ( alpha * x + (1-alpha) )^2 */
268  X2 *= (real)(alpha*alpha);
269  X2 += ((real)2.0 * alpha * (1-alpha)) * X;
270  X2.add_identity( (1-alpha)*(1-alpha) );
271  }
272  // Apply polynomial to homo and lumo intervals
273  homo.puriStep(poly, alpha);
274  lumo.puriStep(poly, alpha);
275  // Transfer content of X2 to X clearing previous content of X if any
276  // In current implementation this is needed regardless of which
277  // polynomial is used
278  X2.transfer(X);
279 
280  int steps_left = estimated_steps_left( homo, lumo,
282  2*puri_steps.size(),
283  use_scaling );
284  real acc_error = accumulated_error();
285  assert(acc_error >= 0);
286  real chosen_thresh = get_threshold( homo.low() - lumo.upp(),
287  steps_left, acc_error,
289  Time t_thr; t_thr.tic();
290  real actual_thresh = X.thresh( chosen_thresh, normTruncation );
291  double wall_sec_thresh = t_thr.toc();
292  // Eigenvalues may have moved due to truncation. Therefore, we
293  // increase the intervals containing the HOMO and LUMO
294  // eigenvalues according to Weyl's theorem.
295  homo.increase( actual_thresh );
296  lumo.increase( actual_thresh );
297 
298  real const ONE = 1.0;
299  Time t_sqr; t_sqr.tic();
300  X2 = ONE * X * X;
301  double wall_sec_square = t_sqr.toc();
302  bool homo_was_computed = false;
303  bool lumo_was_computed = false;
304  Time t_norm; t_norm.tic();
306  homo_was_computed,
307  lumo_was_computed );
308  double wall_sec_XmX2norm = t_norm.toc();
309  ++current_step;
310  double wall_sec_total = t_tot.toc();
311  puri_steps[current_step] = Step<real>( poly, alpha,
312  homo, lumo,
313  X.trace(), X2.trace(),
314  chosen_thresh,
315  actual_thresh,
316  X.nnz(), X2.nnz(),
317  wall_sec_thresh,
318  wall_sec_square,
319  wall_sec_XmX2norm,
320  wall_sec_total);
321  if ( homo_was_computed ) {
323  homo_is_computed = true;
324  }
325  if ( lumo_was_computed ) {
327  lumo_is_computed = true;
328  }
329  } // end while
330  } // end purify()
331 
332  template<typename Tmatrix>
335  mat::Interval<real> lumo,
336  real const toleratedEigenvalError,
337  int const max_steps,
338  bool const use_scaling) {
339 
340  int steps = 1;
341  // There is always at least one truncation left to do and to be on
342  // the safe side we expect that truncation will result in an extra
343  // step which means that the estimated number of steps is always
344  // larger than or equal to 2.
345  while ((1 - homo.low() > toleratedEigenvalError &&
346  lumo.upp() > toleratedEigenvalError) && steps < max_steps) {
347  steps++;
348  int poly = get_poly( homo.low(), lumo.upp() );
349  real alpha = 1;
350  if (use_scaling) {
351  if (poly) {
352  alpha = 2 / ( 1+homo.upp() );
353  }
354  else {
355  alpha = 2 / ( 2-lumo.low() );
356  }
357  }
358  homo.puriStep(poly, alpha);
359  lumo.puriStep(poly, alpha);
360  } // end while
361  return steps;
362  }
363 
364  template<typename Tmatrix>
366  get_threshold( real const gap,
367  int const steps_left,
368  real const acc_error,
369  real const toleratedSubspaceError ) {
370  real acceptable_error = toleratedSubspaceError - acc_error;
371  /*
372  Set number of steps left to be at least smallest_factor since
373  then there will always be some room left for another
374  truncation. That is, we can make some error in the last
375  iterations in case the number of steps was underestimated. The
376  number smallest_factor can be replaced but should always be
377  larger than 1. This number can be used to control the
378  convergence of homo and lumo eigenvalues in the final
379  iterations.
380  */
381  real smallest_factor = 1.01;
382  real steps_left_pessimistic = steps_left > smallest_factor ? steps_left : smallest_factor;
383  assert(acceptable_error >= 0);
384  return ( gap * acceptable_error/steps_left_pessimistic ) /
385  ( 1 + acceptable_error/steps_left_pessimistic );
386  }
387 
388 #if 1
389  // The following new version differs from the previous one in that
390  // the norm ||X-X^2|| is computed accurately not only if we are sure to
391  // compute a wanted eigenvalue but also if there is a chance that we
392  // do so. Sometimes, if the computed eigenvalue is only in one of
393  // the previously known homo and lumo eigenvalue intervals, one can
394  // after the computation know which eigenvalue was
395  // computed. However, sometimes it is not possible to know even after
396  // computation which eigenvalue was computed, and then the
397  // computation was done in vain. In the previous purification
398  // version of 2008 (or earlier?) the information (the norm value)
399  // was stored since further information may become available in
400  // later iterations. That approach is not used here.
401  template<typename Tmatrix>
404  mat::Interval<real>& lumo,
405  bool & homo_was_computed,
406  bool & lumo_was_computed ) {
407  homo_was_computed = false;
408  lumo_was_computed = false;
409  bool computeAccurately = false;
410  bool homo_may_be_computed = false;
411  bool lumo_may_be_computed = false;
412  if ( homo.low() > 0.5 && homo.low() + lumo.low() < 1 ) {
413  homo_may_be_computed = true;
414  // Necessary but not sufficient conditions to compute homo
415  // eigenvalue satisfied
416  if ( !homo_is_computed )
417  computeAccurately = true;
418  }
419  if ( lumo.upp() < 0.5 && homo.upp() + lumo.upp() > 1 ) {
420  lumo_may_be_computed = true;
421  // Necessary but not sufficient conditions to compute lumo
422  // eigenvalue satisfied
423  if ( !lumo_is_computed )
424  computeAccurately = true;
425  }
426  // real const XmX2ENIsSmallValue = 0.207106781186547;
427  // A value here of 0.2475 means that the interval [0.45, 0.55]
428  // must be empty from eigenvalues.
429  // (0.45-0.45*0.45 == 0.55-0.55*0.55 == 0.2475)
430  real const XmX2ENIsSmallValue = 0.2475;
431  real diffAcc = puri_steps[current_step].chosen_thresh / 100;
432  mat::Interval<real> normXmX2;
433  if ( computeAccurately )
434  normXmX2 = Tmatrix::diffIfSmall(X, X2,
435  mat::euclNorm, diffAcc,
436  XmX2ENIsSmallValue);
437  else
438  normXmX2 = Tmatrix::diffIfSmall(X, X2,
439  mat::frobNorm, diffAcc,
440  XmX2ENIsSmallValue);
441 #if 0
442  // Code for computing upper bound using mixed norm. Not clear if
443  // this would have any effect.
444  real norm_upper_bound = Tmatrix::mixed_diff(X, X2, diffAcc);
445  norm_upper_bound = norm_upper_bound < normXmX2.upp() ? norm_upper_bound : normXmX2.upp();
446  normXmX2 = mat::Interval<real>(normXmX2.low(), norm_upper_bound);
447 #endif
448  // The norm normXmX2 must be smaller than or equal to 0.25
449  {
450  real lowerBound = normXmX2.low();
451  real upperBound = normXmX2.upp() < (real)0.25 ? normXmX2.upp() : (real)0.25;
452  normXmX2 = mat::Interval<real>( lowerBound , upperBound );
453  }
454 
455  // Try to improve homo or lumo
456  mat::Interval<real> tmp = mat::sqrtInt((normXmX2 * (real)(-4.0)) + (real)1.0);
457  // Let's see if we can see if we have computed the homo or lumo
458  // eigenvalue (or perhaps none of them)
459  //
460  // if we have computed homo and homo > 1/2 =>
461  mat::Interval<real> homoTmp = (tmp + (real)1.0) / (real)2.0;
462  // if we have computed lumo and lumo < 1/2 =>
463  mat::Interval<real> lumoTmp =
464  ((tmp * (real)(-1.0)) + (real)1.0) / (real)2.0;
465  if (homo_may_be_computed && !lumo.overlap(lumoTmp) ) {
466  homo.intersect_always_non_empty(homoTmp);
467  if ( computeAccurately && normXmX2.length() <= 2.1*diffAcc )
468  homo_was_computed = true;
469  }
470  if (lumo_may_be_computed && !homo.overlap(homoTmp)) {
471  lumo.intersect_always_non_empty(lumoTmp);
472  if ( computeAccurately && normXmX2.length() <= 2.1*diffAcc )
473  lumo_was_computed = true;
474  }
475  } // end new version of improve_homo_lumo_based_on_normXmX2
476 
477 #else
478  template<typename Tmatrix>
481  mat::Interval<real>& lumo,
482  bool & homo_was_computed,
483  bool & lumo_was_computed ) {
484  homo_was_computed = false;
485  lumo_was_computed = false;
486  bool computeAccurately = false;
487  bool homoIsComputable = false;
488  bool lumoIsComputable = false;
489  if ( homo.upp() + lumo.upp() < 1 ) {
490  // Possibly homo can be computed
491  // Check if homo can be and needs to be computed
492  if (homo.low() > 0.5) {
493  homoIsComputable = true;
494  if ( !homo_is_computed )
495  computeAccurately = true;
496  }
497  }
498  if ( homo.low() + lumo.low() > 1 ){
499  // Possibly lumo can be computed
500  // Check if lumo can be and needs to be computed
501  if (lumo.upp() < 0.5) {
502  lumoIsComputable = true;
503  if ( !lumo_is_computed )
504  computeAccurately = true;
505  }
506  }
507 
508  // real const XmX2ENIsSmallValue = 0.207106781186547;
509  // A value here of 0.2475 means that the interval [0.45, 0.55]
510  // must be empty from eigenvalues.
511  // (0.45-0.45*0.45 == 0.55-0.55*0.55 == 0.2475)
512  real const XmX2ENIsSmallValue = 0.2475;
513  real diffAcc = puri_steps[current_step].chosen_thresh / 100;
514  mat::Interval<real> normXmX2;
515  if ( computeAccurately )
516  normXmX2 = Tmatrix::diffIfSmall(X, X2,
517  mat::euclNorm, diffAcc,
518  XmX2ENIsSmallValue);
519  else
520  normXmX2 = Tmatrix::diffIfSmall(X, X2,
521  mat::frobNorm, diffAcc,
522  XmX2ENIsSmallValue);
523  // The norm normXmX2 must be smaller than or equal to 0.25
524  {
525  real lowerBound = normXmX2.low();
526  real upperBound = normXmX2.upp() < (real)0.25 ? normXmX2.upp() : (real)0.25;
527  normXmX2 = mat::Interval<real>( lowerBound , upperBound );
528  }
529  // Try to improve homo or lumo
530  mat::Interval<real> tmp = mat::sqrtInt((normXmX2 * (real)(-4.0)) + (real)1.0);
531  if (homoIsComputable) {
532  // homo > 1/2 =>
533  mat::Interval<real> homoTmp = (tmp + (real)1.0) / (real)2.0;
534  homo.intersect_always_non_empty(homoTmp);
535  if ( computeAccurately && normXmX2.length() <= 2.1*diffAcc )
536  homo_was_computed = true;
537  }
538  if (lumoIsComputable) {
539  // lumo < 1/2 =>
540  mat::Interval<real> lumoTmp =
541  ((tmp * (real)(-1.0)) + (real)1.0) / (real)2.0;
542  lumo.intersect_always_non_empty(lumoTmp);
543  if ( computeAccurately && normXmX2.length() <= 2.1*diffAcc )
544  lumo_was_computed = true;
545  }
546  } // end get_normXmX2
547 #endif
548 
549  template<typename Tmatrix>
551  real total_error = 0;
552  for (unsigned int step = 0; step <= current_step; step++) {
553  real normE = puri_steps[step].actual_thresh;
554  real gap = puri_steps[step].eig_homo.low() - puri_steps[step].eig_lumo.upp();
555  total_error += normE / (gap - normE);
556  }
557  return total_error;
558  }
559 
560  template<typename Tmatrix>
562  for (unsigned int step = current_step; step >= 1; step-- )
563  puri_steps[step].propagate_homo_to_previous( puri_steps[step-1] );
564  }
565  template<typename Tmatrix>
567  for (unsigned int step = current_step; step >= 1; step-- )
568  puri_steps[step].propagate_lumo_to_previous( puri_steps[step-1] );
569  }
570 
571  template<typename Tmatrix>
575  return ( 1 - homo.low() <= tolEigenvalError &&
576  lumo.upp() <= tolEigenvalError );
577  }
578 
579  template<typename Tmatrix>
582  mat::Interval<real> & luF) {
583  if ( !converged() )
584  throw "Attempt to get homo and lumo eigenvalue intervals when purification has not converged.";
585  mat::Interval<real> homo = puri_steps[0].eig_homo;
586  mat::Interval<real> lumo = puri_steps[0].eig_lumo;
587  homo.increase( puri_steps[0].actual_thresh );
588  lumo.increase( puri_steps[0].actual_thresh );
589  real lmin = eigFInt.low();
590  real lmax = eigFInt.upp();
591  hoF = homo * (lmin - lmax) + lmax;
592  luF = lumo * (lmin - lmax) + lmax;
593  }
594 
595  template<typename Tmatrix>
597  mInfo(std::ostream & s) const {
598  s<<"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
599  s<<"% Norm for truncation: "<< mat::getNormTypeString(normTruncation)
600  << std::endl;
601  s<<"n = "<<n<<";"<<std::endl;
602  s<<"tolSubspaceError = "<<tolSubspaceError<<";"<<std::endl;
603  s<<"tolEigenvalError = "<<tolEigenvalError<<";"<<std::endl;
604  s<<"nIter = "<<current_step+1<<";"<<std::endl;
605 
606  s<<std::setprecision(15);
607  s<<"% traceX and traceX2 \n";
608  s<<"traces = [";
609  for (unsigned int ind = 0; ind <= current_step; ++ind) {
610  s<< puri_steps[ind].traceX << " "
611  << puri_steps[ind].traceX2 << std::endl;
612  }
613  s<<"];\n";
614  s<<"% poly \n";
615  s<<"poly = [";
616  for (unsigned int ind = 0; ind <= current_step; ++ind) {
617  s<< puri_steps[ind].poly << std::endl;
618  }
619  s<<"];\n";
620  s<<"% chosenThresh actualThresh \n";
621  s<<"thresh = [";
622  for (unsigned int ind = 0; ind <= current_step; ++ind) {
623  s<< puri_steps[ind].chosen_thresh << " "
624  << puri_steps[ind].actual_thresh << std::endl;
625  }
626  s<<"];\n";
627  s<<"% nnzX nnzX2 \n";
628  s<<"nnzVals = [";
629  for (unsigned int ind = 0; ind <= current_step; ++ind) {
630  s<< puri_steps[ind].nnzX << " "
631  << puri_steps[ind].nnzX2 << std::endl;
632  }
633  s<<"];\n";
634  s<<"% lumo_low lumo_upp homo_low homo_upp \n";
635  s<<"homo_lumo_eigs = [";
636  for (unsigned int ind = 0; ind <= current_step; ++ind) {
637  s<< puri_steps[ind].eig_lumo_orig.low() << " "
638  << puri_steps[ind].eig_lumo_orig.upp() << " "
639  << puri_steps[ind].eig_homo_orig.low() << " "
640  << puri_steps[ind].eig_homo_orig.upp() << std::endl;
641  }
642  s<<"];\n";
643 
644  s<<"figure; \n"
645  <<"plot(1:nIter,(homo_lumo_eigs(:,3)+homo_lumo_eigs(:,4))/2,'-b',"
646  << "1:nIter,(homo_lumo_eigs(:,1)+homo_lumo_eigs(:,2))/2,'-r')\n"
647  <<"legend('HOMO','LUMO'),xlabel('Iteration')\n"
648  <<"hold on\n"
649  <<"for ind = 1:nIter \n"
650  <<" plot([ind ind],[homo_lumo_eigs(ind,3) homo_lumo_eigs(ind,4)], '-b')\n"
651  <<" plot([ind ind],[homo_lumo_eigs(ind,1) homo_lumo_eigs(ind,2)], '-r')\n"
652  <<"end\n"
653  <<"axis([0 nIter 0 1])\n";
654 
655  s<<"figure; \n"
656  <<"subplot(211)\n"
657  <<"plot(1:nIter,100*nnzVals(:,1)/(n*n),'o-r',...\n"
658  <<"1:nIter, 100*nnzVals(:,2)/(n*n),'x-b')\n"
659  <<"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
660  <<"ylabel('Percentage')\n"
661  <<"axis([0 nIter 0 100])\n";
662  s<<"subplot(212)\n"
663  <<"semilogy(1:nIter,thresh(:,1),'x-r',1:nIter,thresh(:,2),'o-b')\n"
664  <<"xlabel('Iteration'), ylabel('Threshold')\n"
665  <<"legend('Chosen threshold', 'Actual threshold')\n"
666  <<"axis([0 nIter min(thresh(:,2))/10 max(thresh(:,1))*10])\n";
667 
668  }
669 
670  template<typename Tmatrix>
672  mTime(std::ostream & s) const {
673  s<<"% PURIFICATION TIMINGS IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
674  s<<"% Norm for truncation: "<< mat::getNormTypeString(normTruncation)
675  << std::endl;
676  s<<"n = "<<n<<";"<<std::endl;
677  s<<"tolSubspaceError = "<<tolSubspaceError<<";"<<std::endl;
678  s<<"tolEigenvalError = "<<tolEigenvalError<<";"<<std::endl;
679  s<<"wall_sec_thresh = [";
680  for (unsigned int ind = 0; ind <= current_step; ++ind) {
681  s << puri_steps[ind].wall_sec_thresh
682  << " " << std::endl;
683  }
684  s<<"];\n";
685  s<<"wall_sec_square = [";
686  for (unsigned int ind = 0; ind <= current_step; ++ind) {
687  s << puri_steps[ind].wall_sec_square
688  << " " << std::endl;
689  }
690  s<<"];\n";
691  s<<"wall_sec_XmX2norm = [";
692  for (unsigned int ind = 0; ind <= current_step; ++ind) {
693  s << puri_steps[ind].wall_sec_XmX2norm
694  << " " << std::endl;
695  }
696  s<<"];\n";
697  s<<"wall_sec_total = [";
698  for (unsigned int ind = 0; ind <= current_step; ++ind) {
699  s << puri_steps[ind].wall_sec_total
700  << " " << std::endl;
701  }
702  s<<"];\n";
703  // Plotting
704  s<<"puriTime = [wall_sec_square wall_sec_thresh wall_sec_XmX2norm];\n";
705  s<<"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
706  "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
707  " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
708  s<<"figure; plot(wall_sec_total,'-x')"<<std::endl;
709  }
710 
711  template<typename Tmatrix>
712  class Purification_scaled<Tmatrix>::Time {
713  public:
714  Time() {}
715  void tic() {ticTime = get_wall_seconds();}
716  double toc() {return get_wall_seconds() - ticTime;}
717  private:
718  double ticTime;
719  static double get_wall_seconds() {
720  struct timeval tv;
721  if(gettimeofday(&tv, NULL) != 0)
722  throw std::runtime_error("Error in get_wall_seconds(), in gettimeofday().");
723  double seconds = tv.tv_sec + (double)tv.tv_usec / 1000000;
724  return seconds;
725  }
726  };
727 
728 
729 } /* end namespace pur */
730 #endif
std::string getNormTypeString(normType nType)
Definition: matInclude.cc:55
bool overlap(Interval const &other) const
Definition: Interval.h:123
Treal upp() const
Definition: Interval.h:143
Definition: Purification_scaled.h:46
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
Step class.
bool converged()
Definition: Purification_scaled.h:572
real const tolEigenvalError
Tolerated error in eigenvalues at convergence.
Definition: Purification_scaled.h:152
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
static double get_wall_seconds()
Definition: bench_gemm_only.cc:47
double ticTime
Definition: Purification_scaled.h:718
Proxy structs used by the matrix API.
void mInfo(std::ostream &file) const
Definition: Purification_scaled.h:597
std::vector< Step< real > > puri_steps
Definition: Purification_scaled.h:150
double toc()
Definition: Purification_scaled.h:716
Tmatrix X2
Definition: Purification_scaled.h:148
Time()
Definition: Purification_scaled.h:714
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:131
bool empty() const
Definition: Interval.h:49
Definition: allocate.cc:30
bool use_scaling
Definition: Purification_scaled.h:158
static int get_poly(real const homo, real const lumo)
Definition: Purification_scaled.h:88
mat::Interval< real > const & eigFInt
Definition: Purification_scaled.h:151
void purify()
Definition: Purification_scaled.h:243
bool homo_is_computed
Definition: Purification_scaled.h:159
void puriStep(int poly)
Definition: Interval.h:228
Definition: matInclude.h:135
static real get_threshold(real const gap, int const steps_left, real const acc_error, real const toleratedSubspaceError)
Compute acceptable error due to truncation in an iteration based on the current band gap...
Definition: Purification_scaled.h:366
mat::normType normTruncation
Definition: Purification_scaled.h:157
unsigned int current_step
Definition: Purification_scaled.h:156
void propagate_homo_information()
Go through step vector and improve homo information in each step possible.
Definition: Purification_scaled.h:561
Treal low() const
Definition: Interval.h:142
void propagate_lumo_information()
Go through step vector and improve lumo information in each step possible.
Definition: Purification_scaled.h:566
static int estimated_steps_left(mat::Interval< real > eig_homo, mat::Interval< real > eig_lumo, real const toleratedEigenvalError, int const max_steps, bool const use_scaling)
Estimates number of remaining iterations.
Definition: Purification_scaled.h:334
Interval< Treal > sqrtInt(Interval< Treal > const &other)
Definition: Interval.h:217
Purification_scaled(Tmatrix &F_and_D, mat::Interval< real > const &eigFInt, mat::Interval< real > const &hoF, mat::Interval< real > const &luF, real const toleratedEigenvalError, real const toleratedSubspaceError, int const max_steps, mat::normType normForTruncation, bool const use_scaling)
Definition: Purification_scaled.h:165
void get_homo_lumo_intervals(mat::Interval< real > &hoF, mat::Interval< real > &luF)
Computed eigenvalues of F.
Definition: Purification_scaled.h:581
void mTime(std::ostream &file) const
Definition: Purification_scaled.h:672
Definition: Purification_scaled.h:712
ergo_real real
Definition: cubature_rules.h:33
Definition: Step.h:42
real accumulated_error()
Total accumulated error due to removal of small matrix elements up to current iteration.
Definition: Purification_scaled.h:550
Tmatrix & X
Definition: Purification_scaled.h:147
static double get_wall_seconds()
Definition: Purification_scaled.h:719
void tic()
Definition: Purification_scaled.h:715
bool lumo_is_computed
Definition: Purification_scaled.h:160
Copyright(c) Emanuel Rubensson 2006.
Definition: matInclude.h:135
real const tolSubspaceError
Tolerated subspace error.
Definition: Purification_scaled.h:154
void intersect_always_non_empty(Interval const &other)
Definition: Interval.h:78
Tmatrix::real real
Definition: Purification_scaled.h:48
int n
Definition: Purification_scaled.h:149
Definition: Purification_scaled.h:44
normType
Definition: matInclude.h:135
void improve_homo_lumo_based_on_normXmX2(mat::Interval< real > &homo, mat::Interval< real > &lumo, bool &homo_was_computed, bool &lumo_was_computed)
Computes interval containing spectral norm of X-X^2 matrix.
Definition: Purification_scaled.h:403