ergo
purification_general.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4  * and Anastasia Kruchinina.
5  *
6  * This program is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program. If not, see <http://www.gnu.org/licenses/>.
18  *
19  * Primary academic reference:
20  * Ergo: An open-source program for linear-scaling electronic structure
21  * calculations,
22  * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23  * Kruchinina,
24  * SoftwareX 7, 107 (2018),
25  * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26  *
27  * For further information about Ergo, see <http://www.ergoscf.org>.
28  */
29 
39 #ifndef HEADER_PURIFICATION_GENERAL
40 #define HEADER_PURIFICATION_GENERAL
41 
42 #include <iostream>
43 #include <fstream>
44 #include <sstream>
45 
46 #include "matrix_typedefs.h" // definitions of matrix types and interval type
47 #include "realtype.h" // definitions of types
48 #include "matrix_utilities.h"
50 #include "output.h"
51 #include "matrix_proxy.h"
52 
53 #include "puri_info.h"
54 #include "constants.h"
55 #include "utilities.h"
56 #include "units.h"
57 
58 #include "files_dense.h"
59 #include "files_sparse.h"
60 #include "files_sparse_bin.h"
61 
62 #include "get_eigenvectors.h"
63 
64 typedef ergo_real real;
65 
66 /**************** DEFINE ******************/
67 
68 // number of additional iterations
69 
70 /* After that stopping criterion says to stop, we save matrix X into the file,
71  * perform additional iterations and read matrix X back.
72  * It is done just for testing purposes in case we wish to run a few SCF iterations
73  * and use results in the iterations given by the stopping criterion,
74  * but at the same time we want to see what happen after stoping criterion says
75  * to stop the iterations. */
76 #define NUM_ADDITIONAL_ITERATIONS 0
77 
78 // enable more output
79 //#define DEBUG_PURI_OUTPUT
80 #define PURI_OUTPUT_NNZ
81 
82 
83 // enable printf output instead of output to ergoscf.out
84 // do it if you need run just recursive expansion, not the whole SCF calculations
85 //#define ENABLE_PRINTF_OUTPUT
86 
87 
88 //#define SAVE_MATRIX_IN_EACH_ITERATION // save F and X_i for every i
89 
90 
91 
92 /**********************************/
93 
94 
95 extern real eucl_acc;
96 extern real mixed_acc;
99 extern int EIG_EMPTY;
100 extern int EIG_SQUARE_INT;
101 extern int EIG_PROJECTION_INT;
102 extern int EIG_POWER_INT;
103 extern int EIG_LANCZOS_INT;
104 
105 
110 template<typename MatrixType>
112 {
113 public:
114  typedef ergo_real real;
117  typedef std::vector<int> VectorTypeInt;
118  typedef std::vector<real> VectorTypeReal;
121 
123  {
124  initialized_flag = false;
125  puri_is_prepared_flag = false;
126  computed_spectrum_bounds = false;
127 
128  // eigenvectors staff
131  eigVecHOMO = NULL;
132  eigVecLUMO = NULL;
133  }
134 
137  virtual void initialize(const MatrixType& F_,
138  const IntervalType& lumo_bounds_,
139  const IntervalType& homo_bounds_,
140  int maxit_,
141  real error_sub_,
142  real error_eig_,
143  int use_new_stopping_criterion_,
144  NormType norm_truncation,
145  NormType norm_stop_crit,
146  int nocc_
147  );
148 
149 
152  virtual bool is_initialized() const { return initialized_flag; }
153 
156  virtual bool puri_is_prepared() const { return puri_is_prepared_flag; }
157 
160  virtual void PurificationStart();
161 
163  virtual void prepare_to_purification();
164 
166  virtual void purification_process();
167 
169  virtual void eigenvalue_bounds_estimation();
170 
171  virtual ~PurificationGeneral() {}
172 
177  virtual void clear() { X.clear(); Xsq.clear(); }
178 
179 
184  void set_spectrum_bounds(real eigmin, real eigmax);
185 
189  void get_spectrum_bounds(real& eigmin, real& eigmax);
190 
193 
194 
195  virtual real total_subspace_error(int it);
196 
198  static real get_epsilon()
199  { return mat::getMachineEpsilon<real>(); }
200 
203  { return std::numeric_limits<real>::max(); }
204 
207  { return std::numeric_limits<real>::min(); }
208 
210  void gen_matlab_file_norm_diff(const char *filename) const;
212  void gen_matlab_file_threshold(const char *filename) const;
214  void gen_matlab_file_nnz(const char *filename) const;
216  void gen_matlab_file_eigs(const char *filename) const;
218  void gen_matlab_file_time(const char *filename) const;
220  void gen_matlab_file_cond_num(const char *filename) const;
221 
223  void gen_python_file_nnz(const char *filename) const;
224 
225 
226 
228  void set_eigenvectors_params(string eigenvectors_method_,
229  string eigenvectors_iterative_method_,
230  real eigensolver_accuracy_,
231  int eigensolver_maxiter_,
232  int scf_step_,
233  int use_prev_vector_as_initial_guess_,
234  int try_eigv_on_next_iteration_if_fail_,
235  VectorType *eigVecLUMO_,
236  VectorType *eigVecHOMO_
237  );
238 
241 
242 
243  void compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F); // for testing
244 
245 
253 protected:
254 
259  void save_matrix_now(string str);
260 
263 
266  void compute_X();
267 
269  void map_bounds_to_0_1();
270 
275  virtual void check_standard_stopping_criterion(const real XmX2_norm, int& stop);
276 
283  virtual void check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace,
284  int& stop, real& estim_order);
285 
287  virtual void stopping_criterion(IterationInfo& iter_info, int& stop, real& estim_order);
288 
291 
301 
302  void compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo);
303 
304 
306  double get_nnz_X(size_t& nnzX )
307  { nnzX = X.nnz(); return (double)(((double)nnzX) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
308 
310  double get_nnz_X()
311  { return (double)(((double)X.nnz()) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
312 
314  double get_nnz_Xsq(size_t& nnzXsq )
315  { nnzXsq = Xsq.nnz(); return (double)(((double)nnzXsq) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
316 
318  double get_nnz_Xsq()
319  { return (double)(((double)Xsq.nnz()) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
320 
325  void estimate_homo_lumo(const VectorTypeReal& XmX2_norm_mixed,
326  const VectorTypeReal& XmX2_norm_frob,
327  const VectorTypeReal& XmX2_trace);
328 
329 
330  void get_frob_norm_est(const VectorTypeReal& XmX2_norm_frob,
331  const std::vector<real>& h_in,
332  const std::vector<real>& l_in,
333  VectorTypeReal& YmY2_norm_frob_est);
334 
335 
336 
337  void get_eigenvalue_estimates(const VectorTypeReal& XmX2_norm_mixed,
338  const VectorTypeReal& XmX2_norm_frob,
339  const VectorTypeReal& XmX2_trace);
340 
341 
342  void propagate_values_in_each_iter(real value_unocc, real value_occ,
343  VectorTypeReal& out_unocc,
344  VectorTypeReal& out_occ,
345  int nmax);
346 
347 
349 
350  void get_iterations_for_lumo_and_homo(int& chosen_iter_lumo,
351  int& chosen_iter_homo);
352 
356 
357  void output_norms_and_traces(IterationInfo& iter_info) const;
358 
359  void get_interval_with_lambda(real& eigVal, VectorType& eigVec, bool& is_homo, bool& is_lumo, const int iter);
360  void get_eigenvalue_of_F_from_eigv_of_Xi(real& eigVal, const VectorType& eigVec);
361 
362  void save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it);
363 
365 
367  void find_shifts_every_iter();
369 
370 
371  void writeToTmpFile(MatrixType& A) const;
372  void readFromTmpFile(MatrixType& A) const;
373 
374  /*
375  * virtual void update_bounds(const real value);
376  * real compute_chemical_potential(const int it);
377  * real get_lower_bound_in_estim_homo_lumo(const int it);
378  */
379 
380  virtual void set_init_params() = 0;
381  virtual void truncate_matrix(real& thresh, int it) = 0;
382  virtual void estimate_number_of_iterations(int& estim_num_iter) = 0;
383  virtual void purify_X(const int it) = 0;
384  virtual void purify_bounds(const int it) = 0;
385  virtual void save_other_iter_info(IterationInfo& iter_info, int it) = 0;
386  virtual void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it) = 0;
387  virtual void return_constant_C(const int it, real& Cval) = 0;
388 
389  /*virtual real apply_inverse_poly(const int it, real x) = 0;*/
390  virtual real apply_poly(const int it, real x) = 0;
391  virtual void apply_poly_to_eigs(const int it, real& homo, real& lumo) = 0;
392  virtual real compute_derivative(const int it, real x, real& DDf) = 0;
393 
394 
395 
396  /* PARAMETERS */
397 
400 
404  int maxit;
407  int nocc;
445  /*EIGENVECTORS STAFF*/
446 
448 
460 
461 
464 
465 
466  /* Eigenvectors */
470  // accuracy of the eigenvalue problem solver
471  // when residual is less then this value,stop
473  // maximum number of iterations for eigensolver
475 
478 
480 
483 
486 
489 
490 
493 
496 
503  int scf_step;
504 
512 };
513 
514 /**************************************/
515 
516 
517 /******************* INIT *******************/
518 
519 
520 template<typename MatrixType>
522  const IntervalType& lumo_bounds_,
523  const IntervalType& homo_bounds_,
524  int maxit_,
525  real error_sub_,
526  real error_eig_,
527  int use_new_stopping_criterion_,
528  NormType norm_truncation_,
529  NormType norm_stop_crit_,
530  int nocc_
531  )
532 {
533  X = F_;
534  maxit = maxit_;
535  assert(maxit >= 1);
536  error_sub = error_sub_;
537  error_eig = error_eig_;
538  use_new_stopping_criterion = use_new_stopping_criterion_;
539  normPuriTrunc = norm_truncation_;
540  normPuriStopCrit = norm_stop_crit_;
541  nocc = nocc_;
542 
543  initialized_flag = true;
544 
545  // save bounds for the matrix F
546  lumo_bounds_F = lumo_bounds_;
547  homo_bounds_F = homo_bounds_;
548 
549  /* Use this function to enable printf output of the purification
550  * work if you want to run just the purification, not the whole scf calculations */
551 #ifdef ENABLE_PRINTF_OUTPUT
553 #endif
554 
555  size_t nnzX;
556  double nnzXp = get_nnz_X(nnzX);
557  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Creating purification object: N = %d"
558  " , nocc = %d , NNZ = %lu <-> %.5lf %%",
559  X.get_nrows(), nocc, nnzX, nnzXp);
560 
561 
562  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen norm for the truncation: ");
563  switch (normPuriTrunc)
564  {
565  case mat::mixedNorm:
567  break;
568 
569  case mat::euclNorm:
571  break;
572 
573  case mat::frobNorm:
575  break;
576 
577  default:
578  throw std::runtime_error("Unknown norm in PurificationGeneral");
579  }
580 
581 
582 #ifdef USE_CHUNKS_AND_TASKS
583  if ((normPuriTrunc == mat::mixedNorm) || (normPuriTrunc == mat::euclNorm))
584  {
585  normPuriTrunc = mat::frobNorm;
586  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change norm for the truncation to Frobenius.");
587  }
588 #endif
589 
590  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen norm for the stopping criterion: ");
591  switch (normPuriStopCrit)
592  {
593  case mat::mixedNorm:
595  break;
596 
597  case mat::euclNorm:
599  break;
600 
601  case mat::frobNorm:
603  break;
604 
605  default:
606  throw std::runtime_error("Unknown norm in PurificationGeneral");
607  }
608 
609 
610 #ifdef USE_CHUNKS_AND_TASKS
611  if ((normPuriStopCrit == mat::mixedNorm) || (normPuriStopCrit == mat::euclNorm))
612  {
613  normPuriStopCrit = mat::frobNorm;
614  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change norm the stopping criterion to Frobenius.");
615  }
616 #endif
617 
618  if (this->use_new_stopping_criterion == 1)
619  {
620  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen the NEW stopping criterion.");
621  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Allowed error in subspace %e", (double)error_sub);
622  }
623  else
624  {
625  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen the OLD stopping criterion.");
626  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Allowed error in subspace %e , in eigenvalues %e", (double)error_sub, (double)error_eig);
627  }
628 
629  check_stopping_criterion_iter = -1; // will be set in function estimate_number_of_iterations()
630  compute_eigenvectors_in_this_SCF_cycle = false; // can be set to true in prepare_to_purification()
631 
632  set_init_params();
633 
634  additional_iterations = NUM_ADDITIONAL_ITERATIONS;
635  info.stopping_criterion = use_new_stopping_criterion; // 1 if new, 0 if not
636  info.error_subspace = error_sub;
637 
638  info.debug_output = 0;
639 #ifdef DEBUG_PURI_OUTPUT
640  info.debug_output = 1;
641 #endif
642 }
643 
644 
645 /*******************************************************/
646 
647 
648 template<typename MatrixType>
650  string eigenvectors_iterative_method_,
651  real eigensolver_accuracy_,
652  int eigensolver_maxiter_,
653  int scf_step_,
654  int use_prev_vector_as_initial_guess_,
655  int try_eigv_on_next_iteration_if_fail_,
656  VectorType *eigVecLUMO_,
657  VectorType *eigVecHOMO_
658 )
659 {
660  // before this was an input parameter
661  number_of_eigenvalues = 1;
662  if (number_of_eigenvalues > 1)
663  {
664  throw std::runtime_error("Error in set_eigenvectors_params() : cannot compute more than 1 eigenpair.");
665  }
666 
667  // can be NULL, empty or containing eigenvector from the previous cycle
668  eigVecLUMO = eigVecLUMO_;
669  eigVecHOMO = eigVecHOMO_;
670 
671  eigensolver_accuracy = eigensolver_accuracy_;
672  eigensolver_maxiter = eigensolver_maxiter_;
673  assert(eigensolver_maxiter >= 1);
674  scf_step = scf_step_;
675  eigenvectors_method_str = eigenvectors_method_;
676  eigenvectors_iterative_method_str = eigenvectors_iterative_method_;
677  eigenvectors_method = get_int_eig_method(eigenvectors_method_);
678  eigenvectors_iterative_method = get_int_eig_iter_method(eigenvectors_iterative_method_);
679  use_prev_vector_as_initial_guess = use_prev_vector_as_initial_guess_;
680  try_eigv_on_next_iteration_if_fail = try_eigv_on_next_iteration_if_fail_;
681 
682  info.lumo_eigenvector_is_computed = false;
683  info.homo_eigenvector_is_computed = false;
684 
685  iter_for_homo = -1;
686  iter_for_lumo = -1;
687 
688  // no given method for computing eigenvectors
689  if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && (eigenvectors_method == EIG_EMPTY))
690  {
691 
692  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "No given method for computing eigenvectors."
693  "Eigenvectors will not be computed in this SCF cycle. Set eigenvectors to NULL.");
694 
695  delete eigVecLUMO; // not NULL here
696  delete eigVecHOMO; // not NULL here
697  eigVecLUMO = NULL;
698  eigVecHOMO = NULL;
699  }
700  else
701  {
702  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen method to compute eigenvectors: %s", eigenvectors_method_str.c_str());
703  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen iterative method to compute eigenvectors: %s", eigenvectors_iterative_method_str.c_str());
704  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen eigensolver accuracy: %g", (double)eigensolver_accuracy);
705  }
706 
707  // reuse eigenvector computed in the previous SCF cycle as an initial guess
708  if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && use_prev_vector_as_initial_guess)
709  {
710  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use eigenvectors from the previous SCF cycle as an initial guess for this SCF cycle");
711  }
712 
713 
714 #ifndef USE_CHUNKS_AND_TASKS
715  if (compute_eigenvectors_in_each_iteration)
716  {
717  /* if we compute eigenvectors in every iteration, we save initial guess in the separate vector eigVecLUMORef*/
718  if ((eigVecLUMO != NULL) && use_prev_vector_as_initial_guess)
719  {
721  if (X.is_empty())
722  {
723  throw std::runtime_error("Error in set_eigenvectors_params() : cannot save initial guess for LUMO!");
724  }
725  X.getCols(cols);
726  eigVecLUMORef.resetSizesAndBlocks(cols);
727  eigVecLUMORef = *eigVecLUMO;
728  }
729 
730  if ((eigVecHOMO != NULL) && use_prev_vector_as_initial_guess)
731  {
733  if (X.is_empty())
734  {
735  throw std::runtime_error("Error in set_eigenvectors_params() : cannot save initial guess for HOMO!");
736  }
737  X.getCols(cols);
738  eigVecHOMORef.resetSizesAndBlocks(cols);
739  eigVecHOMORef = *eigVecHOMO;
740  }
741  }
742 #endif
743 }
744 
745 
746 template<typename MatrixType>
748 {
749  if (eigenvectors_method == "square")
750  {
751  return EIG_SQUARE_INT;
752  }
753  if (eigenvectors_method == "projection")
754  {
755  return EIG_PROJECTION_INT;
756  }
757  if (eigenvectors_method == "")
758  {
759  return EIG_EMPTY;
760  }
761  throw std::runtime_error("Error in get_int_eig_method(): unknown method to compute eigenvectors");
762 }
763 
764 
765 template<typename MatrixType>
766 int PurificationGeneral<MatrixType>::get_int_eig_iter_method(string eigenvectors_iterative_method)
767 {
768  if (eigenvectors_iterative_method == "power")
769  {
770  return EIG_POWER_INT;
771  }
772  if (eigenvectors_iterative_method == "lanczos")
773  {
774  return EIG_LANCZOS_INT;
775  }
776  if (eigenvectors_iterative_method == "")
777  {
778  return EIG_EMPTY;
779  }
780  throw std::runtime_error("Error in get_int_eig_iter_method(): unknown iterative method to compute eigenvectors");
781 }
782 
783 
784 template<typename MatrixType>
786 {
787  real XmX2_fro_norm = iter_info.XmX2_fro_norm;
788  real XmX2_eucl = iter_info.XmX2_eucl;
789  real XmX2_mixed_norm = iter_info.XmX2_mixed_norm;
790  real XmX2_trace = iter_info.XmX2_trace;
791  int it = iter_info.it;
792 
793  assert(it >= 0);
794 
795 #ifndef USE_CHUNKS_AND_TASKS
796  if (normPuriStopCrit == mat::euclNorm)
797  {
798  assert(XmX2_eucl >= 0);
799  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "||X-X^2||_2 = %e", (double)XmX2_eucl);
800  }
801 
802  if (normPuriStopCrit == mat::mixedNorm)
803  {
804  assert(XmX2_fro_norm >= 0);
805  assert(XmX2_mixed_norm >= 0);
806  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "||X-X^2||_F = %e , ||X-X^2||_mixed = %e", (double)XmX2_fro_norm, (double)XmX2_mixed_norm);
807  }
808 #endif
809 
810  if (normPuriStopCrit == mat::frobNorm)
811  {
812  assert(XmX2_fro_norm >= 0);
813  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "||X-X^2||_F = %e", (double)XmX2_fro_norm);
814  }
815 
816  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "trace(X-X^2) = %e", (double)XmX2_trace);
817  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "trace(X) = %e", (double)X.trace());
818 }
819 
820 
821 /****************** PURIFICATION_START ********************/
822 
823 
824 template<typename MatrixType>
826 {
827  Util::TimeMeter total_time; // measure total time of the purification process
828 
829  prepare_to_purification();
830 
831 #ifdef DEBUG_PURI_OUTPUT
832  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Starting recursive expansion");
833 #endif
834  purification_process();
835 
836  eigenvalue_bounds_estimation();
837 
838 
839 
840  /* COMPUTE EIGENVECTORS WITH PROJECTION METHOD */
841  if (info.converged == 1)
842  {
843  if (compute_eigenvectors_in_this_SCF_cycle && (eigenvectors_method == EIG_PROJECTION_INT))
844  {
845  compute_eigenvectors_without_diagonalization_last_iter_proj();
846  }
847  }
848  else
849  {
850  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot compute eigenvectors using projection method since the purification did not converge");
851  }
852 
853  /* CHECK IF EIGENVECTORS ARE CORRECT - EITHER COMPUTED WITH SQUARE METHOD OR PROJECTION */
854  check_eigenvectors_at_the_end();
855 
856  total_time.print(LOG_AREA_DENSFROMF, "Recursive expansion");
857  double total_time_stop = total_time.get_elapsed_wall_seconds();
858  info.total_time = total_time_stop;
859 }
860 
861 
862 template<typename MatrixType>
864 {
865  if (!is_initialized())
866  {
867  throw std::runtime_error("Error in prepare_to_purification() : function is called for an uninitialized class.");
868  }
869 
870  if (!computed_spectrum_bounds)
871  {
872  Util::TimeMeter total_time_spectrum_bounds;
873  compute_spectrum_bounds();
874  total_time_spectrum_bounds.print(LOG_AREA_DENSFROMF, "compute_spectrum_bounds");
875  double total_time_spectrum_bounds_stop = total_time_spectrum_bounds.get_elapsed_wall_seconds();
876  info.time_spectrum_bounds = total_time_spectrum_bounds_stop;
877  }
878 
879  info.set_spectrum_bounds(spectrum_bounds.low(), spectrum_bounds.upp());
880  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Spectrum of F: \t [ %.12lf , %.12lf ]", (double)spectrum_bounds.low(), (double)spectrum_bounds.upp());
881 
882  map_bounds_to_0_1();
883  set_truncation_parameters();
884 
885  if ((eigVecLUMO != NULL) || (eigVecHOMO != NULL))
886  {
887  // check if we have non-overlapping homo and lumo bounds
888  if (1 - homo_bounds.upp() - lumo_bounds.upp() >= TOL_OVERLAPPING_BOUNDS)
889  {
890  compute_eigenvectors_in_this_SCF_cycle = true;
891  info.compute_eigenvectors_in_this_SCF_cycle = compute_eigenvectors_in_this_SCF_cycle;
892  F = X; // Save original matrix F, needed for computation of the Rayleigh quotient.
893  // Quotient is needed to compute eigenvalue corresponding to the eigenvector
894  // and see which eigenvalue it is, homo or lumo.
895  writeToTmpFile(F);
896  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "calling determine_iteration_for_eigenvectors()");
897  determine_iteration_for_eigenvectors(); // on which iterations we should compute eigenvectors
898  }
899  else
900  {
901  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Homo and lumo inner bounds are too close (< %g), "
902  "homo and lumo eigenvectors will not be computed", (double)TOL_OVERLAPPING_BOUNDS);
903  }
904  }
905 
906  compute_X(); // F->X, put eigenvalues to the [0,1]
907 
908  puri_is_prepared_flag = true;
909 }
910 
911 
912 template<typename MatrixType>
914 {
915  if (!puri_is_prepared())
916  {
917  throw std::runtime_error("Error in purification_process() : "
918  "first expect a call for function prepare_to_purification().");
919  }
920 
921  int it;
922  int stop = 0;
923  real thresh;
924  real Xsquare_time_stop = -1, total_time_stop = -1, trunc_time_stop = -1, purify_time_stop = -1;
925  real frob_diff_time_stop = -1, eucl_diff_time_stop = -1, trace_diff_time_stop = -1, mixed_diff_time_stop = -1, stopping_criterion_time_stop = -1;
926  int maxit_tmp = maxit;
927  real estim_order = -1;
928  real XmX2_trace = -1;
929  real XmX2_fro_norm = -1;
930  real XmX2_mixed_norm = -1;
931  real XmX2_eucl = -1;
932 
933  int already_stop_before = 0; // flag for checking stopping criterion, needed in case additional_iterations != 0
934 
935  info.Iterations.clear();
936  info.Iterations.reserve(100);
937  //std::cout << "Cap: " << this->info.Iterations.capacity() << std::endl;
938 
939  IterationInfo iter_info; // 0-th iterations
940 
941  // 0 iteration
942  it = 0;
943 
944 
945  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "BEFORE ITERATIONS:");
946 
947 
948 
949 #ifdef SAVE_MATRIX_IN_EACH_ITERATION
950  {
951  ostringstream str;
952  str << it;
953  save_matrix_now(str.str());
954  }
955 #endif
956 
957  Util::TimeMeter total_time; // for this iteration
958 
959 #ifdef PURI_OUTPUT_NNZ
960  double nnzX = get_nnz_X();
961 #endif
962 
963  // truncate matrix
964  Util::TimeMeter trunc_time;
965  truncate_matrix(thresh, it);
966  trunc_time.print(LOG_AREA_DENSFROMF, "truncate_matrix");
967  trunc_time_stop = trunc_time.get_elapsed_wall_seconds();
968 
969 #ifdef PURI_OUTPUT_NNZ
970  if((double)thresh >= 0) // value is available
971  {
972  double nnzXafter = get_nnz_X();
973  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e , nnz before %.2lf %% , nnz after %.2lf %%", (double)thresh, nnzX, nnzXafter);
974  }
975  else
976  {
977  double nnzXafter = get_nnz_X();
978  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "nnz before %.2lf %% , nnz after %.2lf %%", nnzX, nnzXafter);
979  }
980 #else
981  if((double)thresh >= 0) // value is available
982  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e", (double)thresh);
983 #endif
984 
986 
987  Util::TimeMeter Xsquare_time;
988  Xsq = (real)1.0 * X * X;
989  Xsquare_time.print(LOG_AREA_DENSFROMF, "square");
990  Xsquare_time_stop = Xsquare_time.get_elapsed_wall_seconds();
991 
992 #ifdef PURI_OUTPUT_NNZ
993  size_t nnzXsq;
994  double nnzXsqp = get_nnz_Xsq(nnzXsq);
995  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NNZ Xsq = %lu <-> %.2lf %%", nnzXsq, nnzXsqp);
996 #endif
997 
998 
999  // compute frob norm of X-X2
1000  Util::TimeMeter frob_diff_time;
1001  XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1002  frob_diff_time.print(LOG_AREA_DENSFROMF, "Frobenius norm of X-X^2");
1003  frob_diff_time_stop = frob_diff_time.get_elapsed_wall_seconds();
1004 
1005  if (normPuriStopCrit == mat::mixedNorm)
1006  {
1007 #ifndef USE_CHUNKS_AND_TASKS
1008  Util::TimeMeter mixed_diff_time;
1009  XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq, mixed_acc);
1010  mixed_diff_time.print(LOG_AREA_DENSFROMF, "Mixed norm of X-X^2");
1011  mixed_diff_time_stop = mixed_diff_time.get_elapsed_wall_seconds();
1012 #endif
1013  }
1014 
1015  // compute trace of X-X2
1016  Util::TimeMeter trace_diff_time;
1017  XmX2_trace = X.trace() - Xsq.trace();
1018  trace_diff_time.print(LOG_AREA_DENSFROMF, "Trace of X-X^2");
1019  trace_diff_time_stop = trace_diff_time.get_elapsed_wall_seconds();
1020 
1021  if (normPuriStopCrit == mat::euclNorm)
1022  {
1023 #ifndef USE_CHUNKS_AND_TASKS
1024  Util::TimeMeter eucl_diff_time;
1025  XmX2_eucl = MatrixType::eucl_diff(X, Xsq, eucl_acc);
1026  eucl_diff_time.print(LOG_AREA_DENSFROMF, "Spectral norm of X-X^2");
1027  eucl_diff_time_stop = eucl_diff_time.get_elapsed_wall_seconds();
1028 #endif
1029  }
1030 
1031 
1032  iter_info.it = it;
1033  iter_info.XmX2_fro_norm = XmX2_fro_norm;
1034  iter_info.XmX2_eucl = XmX2_eucl;
1035  iter_info.XmX2_mixed_norm = XmX2_mixed_norm;
1036  iter_info.XmX2_trace = XmX2_trace;
1037 #ifdef DEBUG_PURI_OUTPUT
1038  output_norms_and_traces(iter_info);
1039 #endif
1040 
1041 
1042  if (compute_eigenvectors_in_this_SCF_cycle)
1043  {
1044  compute_eigenvectors_without_diagonalization(it, iter_info);
1045  }
1046 
1047  ostringstream str_out;
1048  str_out << "Iteration " << it;
1049  total_time.print(LOG_AREA_DENSFROMF, str_out.str().c_str());
1050  total_time_stop = total_time.get_elapsed_wall_seconds();
1051 
1052 // SAVE INFO ABOUT ITERATION
1053  {
1054  iter_info.gap = 1 - homo_bounds.upp() - lumo_bounds.upp(); // = VecGap[0]
1055  iter_info.threshold_X = thresh;
1056  iter_info.Xsquare_time = Xsquare_time_stop;
1057  iter_info.trunc_time = trunc_time_stop;
1058  iter_info.purify_time = 0;
1059  iter_info.NNZ_X = get_nnz_X();
1060  iter_info.NNZ_X2 = get_nnz_Xsq();
1061  iter_info.total_time = total_time_stop;
1062  iter_info.homo_bound_low = homo_bounds.low();
1063  iter_info.lumo_bound_low = lumo_bounds.low();
1064  iter_info.homo_bound_upp = homo_bounds.upp();
1065  iter_info.lumo_bound_upp = lumo_bounds.upp();
1066  stopping_criterion_time_stop = 0; // we are not checking stopping criterion on the 1 iteration
1067  iter_info.stopping_criterion_time = stopping_criterion_time_stop;
1068  iter_info.eucl_diff_time = eucl_diff_time_stop;
1069  iter_info.frob_diff_time = frob_diff_time_stop;
1070  iter_info.mixed_diff_time = mixed_diff_time_stop;
1071  iter_info.trace_diff_time = trace_diff_time_stop;
1072 
1073  save_other_iter_info(iter_info, it);
1074  }
1075  /**************/
1076 
1077  info.Iterations.push_back(iter_info); // add info about 0 iteration
1078 
1079 
1080  output_current_memory_usage(LOG_AREA_DENSFROMF, "Before iteration 1");
1081  Util::TimeMeter timeMeterWriteAndReadAll;
1082  std::string sizesStr = mat::FileWritable::writeAndReadAll();
1083  timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
1084  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
1085 
1086 
1087  it = 1;
1088  while (it <= maxit_tmp)
1089  {
1090  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "ITERATION %d :", it);
1091 
1092  IterationInfo iter_info; // i-th iteration
1093 
1094  Util::TimeMeter total_time; // for this iteration
1095 
1096  Util::TimeMeter purify_time;
1097  purify_X(it);
1098  purify_time.print(LOG_AREA_DENSFROMF, "purify_X");
1099  purify_time_stop = purify_time.get_elapsed_wall_seconds();
1100 
1101 
1102 #ifdef SAVE_MATRIX_IN_EACH_ITERATION
1103  {
1104  ostringstream str;
1105  str << it;
1106  save_matrix_now(str.str());
1107  }
1108 #endif
1109 
1110 #ifdef PURI_OUTPUT_NNZ
1111  double nnzX = get_nnz_X();
1112 #endif
1113 
1114  // truncate matrix
1115  Util::TimeMeter trunc_time;
1116  truncate_matrix(thresh, it);
1117  trunc_time.print(LOG_AREA_DENSFROMF, "truncate_matrix");
1118  trunc_time_stop = trunc_time.get_elapsed_wall_seconds();
1119 
1120 #ifdef PURI_OUTPUT_NNZ
1121  double nnzXafter = get_nnz_X();
1122  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e , nnz before %.2lf %% , nnz after %.2lf %%", thresh, nnzX, nnzXafter);
1123 #else
1124  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e", thresh);
1125 #endif
1126 
1127 
1128  output_current_memory_usage(LOG_AREA_DENSFROMF, "Before X square");
1129 
1130  Util::TimeMeter Xsquare_time;
1131  Xsq = (real)1.0 * X * X;
1132  Xsquare_time.print(LOG_AREA_DENSFROMF, "square");
1133  Xsquare_time_stop = Xsquare_time.get_elapsed_wall_seconds();
1134 #ifdef PURI_OUTPUT_NNZ
1135  size_t nnzXsq;
1136  double nnzXsqp = get_nnz_Xsq(nnzXsq);
1137  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NNZ Xsq = %lu <-> %.2lf %%", nnzXsq, nnzXsqp);
1138 #endif
1139 
1140 
1141  // update bounds for homo and lumo using corresponding polynomial
1142  purify_bounds(it);
1143 
1144  // compute frob norm of X-X2
1145  Util::TimeMeter frob_diff_time;
1146  XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1147  frob_diff_time.print(LOG_AREA_DENSFROMF, "Frobenius norm of X-X^2");
1148  frob_diff_time_stop = frob_diff_time.get_elapsed_wall_seconds();
1149 
1150  if (normPuriStopCrit == mat::mixedNorm)
1151  {
1152 #ifndef USE_CHUNKS_AND_TASKS
1153  Util::TimeMeter mixed_diff_time;
1154  XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq, mixed_acc);
1155  //XmX2_mixed_norm = Xsq.mixed(mixed_acc);
1156  mixed_diff_time.print(LOG_AREA_DENSFROMF, "Mixed norm of X-X^2");
1157  mixed_diff_time_stop = mixed_diff_time.get_elapsed_wall_seconds();
1158 #endif
1159  }
1160 
1161  if (normPuriStopCrit == mat::euclNorm)
1162  {
1163 #ifndef USE_CHUNKS_AND_TASKS
1164  Util::TimeMeter eucl_diff_time;
1165  XmX2_eucl = MatrixType::eucl_diff(X, Xsq, eucl_acc);
1166  eucl_diff_time.print(LOG_AREA_DENSFROMF, "Spectral norm of X-X^2");
1167  eucl_diff_time_stop = eucl_diff_time.get_elapsed_wall_seconds();
1168 #endif
1169  }
1170 
1171  // compute trace of X-X2
1172  Util::TimeMeter trace_diff_time;
1173  XmX2_trace = X.trace() - Xsq.trace();
1174  trace_diff_time.print(LOG_AREA_DENSFROMF, "Trace of X-X^2");
1175  trace_diff_time_stop = trace_diff_time.get_elapsed_wall_seconds();
1176 
1177 
1178  // CHECK FOR A NEGATIVE TRACE
1179  if (XmX2_trace < -1e10) // here is definitively some misconvergence
1180  {
1181  throw std::runtime_error("Error in purification_process() : trace of X-X^2 is negative, seems as a"
1182  " misconvergence of the recursive expansion.");
1183  }
1184 
1185  iter_info.it = it;
1186  iter_info.XmX2_fro_norm = XmX2_fro_norm;
1187  iter_info.XmX2_eucl = XmX2_eucl;
1188  iter_info.XmX2_mixed_norm = XmX2_mixed_norm;
1189  iter_info.XmX2_trace = XmX2_trace;
1190 #ifdef DEBUG_PURI_OUTPUT
1191  output_norms_and_traces(iter_info);
1192 #endif
1193 
1194 
1195 // SAVE INFO ABOUT ITERATION
1196  {
1197  iter_info.threshold_X = thresh;
1198  iter_info.Xsquare_time = Xsquare_time_stop;
1199  iter_info.trunc_time = trunc_time_stop;
1200  iter_info.purify_time = purify_time_stop;
1201 
1202  iter_info.eucl_diff_time = eucl_diff_time_stop;
1203  iter_info.frob_diff_time = frob_diff_time_stop;
1204  iter_info.mixed_diff_time = mixed_diff_time_stop;
1205  iter_info.trace_diff_time = trace_diff_time_stop;
1206  }
1207 
1208 
1209  /* COMPUTE EIGENVECTORS WITH FOLDED SPECTRUM METHOD (SHIFT_AND_SQUARE) */
1210  if (compute_eigenvectors_in_this_SCF_cycle)
1211  {
1212  compute_eigenvectors_without_diagonalization(it, iter_info);
1213  }
1214 
1215 
1216  // check stopping criterion (call function on every iteration
1217  // larger or equal to check_stopping_criterion_iter)
1218  if (it >= check_stopping_criterion_iter)
1219  {
1220  Util::TimeMeter stopping_criterion_time;
1221  stopping_criterion(iter_info, stop, estim_order);
1222  stopping_criterion_time.print(LOG_AREA_DENSFROMF, "stopping_criterion");
1223  stopping_criterion_time_stop = stopping_criterion_time.get_elapsed_wall_seconds();
1224  iter_info.stopping_criterion_time = stopping_criterion_time_stop;
1225  }
1226  else
1227  {
1228  stop = 0;
1229  estim_order = -1;
1230  }
1231 
1232  // if we reach stopping iteration or maximum allowed number or iterations
1233  // and we are not already stop (in case we have additional_iterations != 0)
1234  if ((already_stop_before == 0) && ((stop == 1) || (it == maxit)))
1235  {
1236  if (stop == 1)
1237  {
1238  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "PURIFICATION CONVERGED after %d iterations", it);
1239  info.converged = 1;
1240  }
1241  else // if it == maxit
1242  {
1243  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NOT CONVERGED. Reached the maximum number of iterations %d", maxit);
1244  info.converged = 0;
1245  }
1246 
1247  assert(maxit_tmp <= (int)VecPoly.size());
1248  maxit_tmp = it + additional_iterations;
1249  already_stop_before = 1;
1250  }
1251 
1252  ostringstream str_out;
1253  str_out << "Iteration " << it;
1254  total_time.print(LOG_AREA_DENSFROMF, str_out.str().c_str());
1255  total_time_stop = total_time.get_elapsed_wall_seconds();
1256 
1257 
1258 
1259  /******************/
1260 
1261 
1262  iter_info.NNZ_X = get_nnz_X();
1263  iter_info.NNZ_X2 = get_nnz_Xsq();
1264 
1265  iter_info.homo_bound_low = homo_bounds.low();
1266  iter_info.homo_bound_upp = homo_bounds.upp();
1267  iter_info.lumo_bound_low = lumo_bounds.low();
1268  iter_info.lumo_bound_upp = lumo_bounds.upp();
1269 
1270  iter_info.total_time = total_time_stop;
1271  iter_info.constantC = constant_C;
1272  if (use_new_stopping_criterion)
1273  {
1274  iter_info.order = estim_order;
1275  }
1276 
1277  save_other_iter_info(iter_info, it);
1278 
1279  /*******************/
1280 
1281  info.Iterations.push_back(iter_info); // add info about i-th iteration to the info
1282 
1283  it++;
1284  } //while
1285 
1286 
1287  // output number of non-zeros in the obtained density matrix
1288  double nnzD = get_nnz_X();
1289  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Number of non-zeros in D is %.2lf %%", nnzD);
1290 
1291 
1292  output_current_memory_usage(LOG_AREA_DENSFROMF, "After the last iteration");
1293  Util::TimeMeter timeMeterWriteAndReadAll_end;
1295  timeMeterWriteAndReadAll_end.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
1296  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
1297 
1298 
1299  info.total_it = maxit_tmp;
1300  info.additional_iterations = additional_iterations;
1301 
1302  real acc_err_sub = total_subspace_error(maxit_tmp - additional_iterations);
1303 #ifdef DEBUG_PURI_OUTPUT
1304  if (acc_err_sub != -1)
1305  {
1306  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "TOTAL accumulated subspace error is %e", acc_err_sub);
1307  }
1308 #endif
1309  info.accumulated_error_subspace = acc_err_sub;
1310 }
1311 
1312 
1313 /* EIGENVALUE BOUND ESTIMATION */
1314 
1315 template<typename MatrixType>
1317 {
1318  if (info.converged == 1)
1319  {
1320  // estimate eigenvalues of the matrix F
1321  VectorTypeReal norms_mixed, norms_frob, traces;
1322  // use mixed norm instead of the Frobenius if it is possible
1323  if (normPuriStopCrit == mat::mixedNorm)
1324  {
1325  info.get_vec_mixed_norms(norms_mixed);
1326  }
1327  info.get_vec_frob_norms(norms_frob);
1328  info.get_vec_traces(traces);
1329  get_eigenvalue_estimates(norms_mixed, norms_frob, traces);
1330 
1331 
1332  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Estimated bounds for the eigenvalues for the Fock matrix:");
1333  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO: [ %.12lf , %.12lf ]", (double)lumo_bounds_F_new.low(), (double)lumo_bounds_F_new.upp());
1334  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO: [ %.12lf , %.12lf ]", (double)homo_bounds_F_new.low(), (double)homo_bounds_F_new.upp());
1335  }
1336  else
1337  {
1338  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot estimate eigenvalues since the purification did not converge");
1339  }
1340 
1341 
1342  info.homo_estim_low_F = homo_bounds_F_new.low();
1343  info.homo_estim_upp_F = homo_bounds_F_new.upp();
1344  info.lumo_estim_low_F = lumo_bounds_F_new.low();
1345  info.lumo_estim_upp_F = lumo_bounds_F_new.upp();
1346 }
1347 
1348 
1349 /******************************************************************************************************************************/
1350 
1351 
1352 template<typename MatrixType>
1354 {
1355  if ((eigVecHOMO == NULL) || eigVecHOMO->is_empty())
1356  {
1357  return;
1358  }
1359  info.eigValHOMO = 0;
1360  info.homo_eigenvector_is_computed = false;
1361 
1362  /* NOTE: clear() gives a zero vector, which may not be empty (i.e. still save some data about the structure)!
1363  * Thus we use clear_structure() to get an empty vector. */
1364  eigVecHOMO->clear_structure();
1365 }
1366 
1367 
1368 template<typename MatrixType>
1370 {
1371  if ((eigVecLUMO == NULL) || eigVecLUMO->is_empty())
1372  {
1373  return;
1374  }
1375  info.eigValLUMO = 0;
1376  info.lumo_eigenvector_is_computed = false;
1377  eigVecLUMO->clear_structure();
1378 }
1379 
1380 
1381 template<typename MatrixType>
1383 {
1384  if (info.converged != 1)
1385  {
1386  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Discard all computed eigenvectors since the purification did not converge");
1387  discard_homo_eigenvector();
1388  discard_lumo_eigenvector();
1389  return;
1390  }
1391 
1392 
1393 
1394  if (compute_eigenvectors_in_this_SCF_cycle)
1395  {
1396  int ITER_DIFF = 2;
1397  // check if we passed iter_for_homo iteration
1398  if (!eigVecHOMO->is_empty() && (info.total_it < iter_for_homo))
1399  {
1400  do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "HOMO eigenvector was not computed. Iteration for homo: %d, total number of iterations: %d",
1401  iter_for_homo, info.total_it);
1402  discard_homo_eigenvector();
1403  }
1404  else
1405  {
1406  // if yes, was it one of the last iterations?
1407  if (!eigVecHOMO->is_empty() && (info.total_it - iter_for_homo < ITER_DIFF) && info.homo_eigenvector_is_computed)
1408  {
1409  do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "HOMO eigenvector was computed in one of the last recursive expansion iterations (%d of total %d). "
1410  "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
1411  "thus result may not be accurate!", iter_for_homo, info.total_it);
1412  }
1413  }
1414 
1415  // check if we passed iter_for_lumo iteration
1416  if (!eigVecLUMO->is_empty() && (info.total_it < iter_for_lumo))
1417  {
1418  do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "LUMO eigenvector was not computed. Iteration for lumo: %d, total number of iterations: %d",
1419  iter_for_lumo, info.total_it);
1420  discard_lumo_eigenvector();
1421  }
1422  else
1423  {
1424  // if yes, was it one of the last iterations?
1425  if (!eigVecLUMO->is_empty() && (info.total_it - iter_for_lumo < ITER_DIFF) && info.lumo_eigenvector_is_computed)
1426  {
1427  do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "LUMO eigenvector was computed in one of the last recursive expansion iterations (%d of total %d). "
1428  "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
1429  "thus result may not be accurate!", iter_for_lumo, info.total_it);
1430  }
1431  }
1432 
1433 
1434 
1435  real low_lumo_F_bound = info.lumo_estim_low_F;
1436  real upp_lumo_F_bound = info.lumo_estim_upp_F;
1437  real low_homo_F_bound = info.homo_estim_low_F;
1438  real upp_homo_F_bound = info.homo_estim_upp_F;
1439 
1440  // For small cases bounds can be too good or even slightly wrong
1441  // due to numerical errors; thus we allow a small flexibility
1442  real flex_tolerance = THRESHOLD_EIG_TOLERANCE;
1443 
1444  if (info.homo_eigenvector_is_computed) // check if eigenvector was computed
1445  {
1446  if ((low_homo_F_bound - flex_tolerance <= eigValHOMO) && (eigValHOMO <= upp_homo_F_bound + flex_tolerance))
1447  {
1448  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO eigenvalue is %lf , HOMO bounds are [ %lf , %lf ]",
1449  (double)eigValHOMO, (double)low_homo_F_bound, (double)upp_homo_F_bound);
1450  info.eigValHOMO = eigValHOMO;
1451  }
1452  else
1453  {
1454  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO eigenvalue is outside HOMO bounds [ %lf , %lf ],"
1455  " discard computed HOMO eigenvector.",
1456  (double)low_homo_F_bound, (double)upp_homo_F_bound);
1457  discard_homo_eigenvector();
1458  }
1459  }
1460  else
1461  {
1462  discard_homo_eigenvector();
1463  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO eigenvector is not computed.");
1464  }
1465 
1466  if (info.lumo_eigenvector_is_computed) // check if eigenvector was computed
1467  {
1468  if ((low_lumo_F_bound - flex_tolerance <= eigValLUMO) && (eigValLUMO <= upp_lumo_F_bound + flex_tolerance))
1469  {
1470  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO eigenvalue is %lf , LUMO bounds are [ %lf , %lf ]",
1471  (double)eigValLUMO, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
1472  info.eigValLUMO = eigValLUMO;
1473  }
1474  else
1475  {
1476  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed LUMO eigenvalue is outside LUMO bounds [ %lf , %lf ],"
1477  " discard computed LUMO eigenvector.",
1478  (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
1479  discard_lumo_eigenvector();
1480  }
1481  }
1482  else
1483  {
1484  discard_lumo_eigenvector();
1485  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO eigenvector is not computed.");
1486  }
1487 
1488  if (info.homo_eigenvector_is_computed && info.lumo_eigenvector_is_computed)
1489  {
1490  real gap = eigValLUMO - eigValHOMO;
1491  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO-LUMO gap is %lf = %lf eV", (double)gap, (double)(gap / UNIT_one_eV));
1492  }
1493  }
1494 }
1495 
1496 
1497 /***************** COMPUTE_X *****************/
1498 
1499 template<typename MatrixType>
1501 {
1502  if (spectrum_bounds.empty())
1503  {
1504  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval spectrum_bounds is empty in compute_X().");
1505  }
1506 
1507 #ifdef DEBUG_PURI_OUTPUT
1508  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Put eigenvalues of F to the interval [0,1] in reverse order.");
1509 #endif
1510 
1511  real eigmin = spectrum_bounds.low();
1512  real eigmax = spectrum_bounds.upp();
1513  X.add_identity(-eigmax);
1514  X *= ((real)1.0 / (eigmin - eigmax));
1515 }
1516 
1517 
1518 template<typename MatrixType>
1520 {
1521  if (spectrum_bounds.empty())
1522  {
1523  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval spectrum_bounds is empty in map_bounds_to_0_1().");
1524  }
1525 
1526 #ifdef DEBUG_PURI_OUTPUT
1527  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Transform homo and lumo bounds...");
1528 #endif
1529 
1530  real eigmin = spectrum_bounds.low();
1531  real eigmax = spectrum_bounds.upp();
1532 
1533  // Compute transformed homo and lumo eigenvalues.
1534 
1535  homo_bounds = homo_bounds_F;
1536  lumo_bounds = lumo_bounds_F;
1537  // homo and lumo must be in the [lmin, lmax] interval
1538  homo_bounds.intersect(spectrum_bounds);
1539  lumo_bounds.intersect(spectrum_bounds);
1540 
1541 #ifdef DEBUG_PURI_OUTPUT
1542  if (homo_bounds.empty())
1543  {
1544  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
1545  }
1546  if (lumo_bounds.empty())
1547  {
1548  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
1549  }
1550 #endif
1551 
1552  if (!mat::Interval<real>::intersect(homo_bounds, lumo_bounds).empty())
1553  {
1554  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Bounds for homo and lumo of F are overlapping.");
1555  }
1556 
1557  // homo_bounds : (1-homo) bounds for matrix X0, later for matrix Xi
1558  // homo_bounds_X0 : homo bounds for matrix X0
1559  homo_bounds = (homo_bounds - eigmax) / (eigmin - eigmax);
1560  homo_bounds_X0 = homo_bounds;
1561  homo_bounds = IntervalType(1 - homo_bounds.upp(), 1 - homo_bounds.low());
1562 
1563  // lumo_bounds : lumo bounds for matrix X0, later for matrix Xi
1564  // lumo_bounds_X0 : lumo bounds for matrix X0
1565  lumo_bounds = (lumo_bounds - eigmax) / (eigmin - eigmax);
1566  lumo_bounds_X0 = lumo_bounds;
1567 
1568  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO bounds of X: \t [ %.12lf , %.12lf ]", (double)(1 - homo_bounds.upp()), (double)(1 - homo_bounds.low()));
1569  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO bounds of X: \t [ %.12lf , %.12lf ]", (double)lumo_bounds.low(), (double)lumo_bounds.upp());
1570 
1571 }
1572 
1573 
1574 /**************************************************************************************/
1575 
1576 
1577 template<typename MatrixType>
1579 {
1580  int estim_num_iter = 0;
1581 
1582  estimate_number_of_iterations(estim_num_iter);
1583  info.estim_total_it = estim_num_iter;
1584 
1585  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "ESTIMATED NUMBER OF ITERATIONS IS %d", estim_num_iter);
1586 
1587  if (estim_num_iter < maxit)
1588  {
1589  // Estimated number of iterations estim_num_iter is less than maxit, set maxit to estim_num_iter
1590  maxit = estim_num_iter;
1591  }
1592 
1593  // error due to truncation allowed in each iteration
1594  error_per_it = error_sub / estim_num_iter;
1595 
1596 #ifdef DEBUG_OUTPUT
1597  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "The total allowed subspace error is %e", error_sub);
1598  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Then error due to truncation allowed in each iteration is %e", error_per_it);
1599 #endif
1600 }
1601 
1602 
1603 /****************** SET_SPECTRUM_BOUNDS *****************/
1604 
1605 template<typename MatrixType>
1607 {
1608  spectrum_bounds = IntervalType(eigmin, eigmax);
1609  computed_spectrum_bounds = true;
1610 }
1611 
1612 
1613 /****************** GET_SPECTRUM_BOUNDS *****************/
1614 
1615 template<typename MatrixType>
1617 {
1618  if (!computed_spectrum_bounds)
1619  {
1620  compute_spectrum_bounds();
1621  }
1622  eigmin = spectrum_bounds.low();
1623  eigmax = spectrum_bounds.upp();
1624 }
1625 
1626 
1627 /****************** COMPUTE_SPECTRUM_BOUNDS *****************/
1628 
1629 template<typename MatrixType>
1631 {
1632  // find approximations using Gershgorin bounds
1633  real eigminG, eigmaxG, eigmin, eigmax;
1634 
1635  Util::TimeMeter total_time_spectrum_bounds;
1636  X.gershgorin(eigminG, eigmaxG);
1637  total_time_spectrum_bounds.print(LOG_AREA_DENSFROMF, "gershgorin");
1638 
1639 
1640  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Gershgorin bounds: [ %.12lf , %.12lf ]", (double)eigminG, (double)eigmaxG);
1641 
1642 
1643  /* ELIAS NOTE 2016-02-08: Expand Gershgorin bounds by a very small
1644  * amount to avoid problems of misconvergence in case one of the
1645  * bounds is exact and the gap is zero (in such cases we want
1646  * convergence failure, not convergence to a solution with wrong
1647  * occupation). */
1648  real smallNumberToExpandWith = template_blas_sqrt(mat::getMachineEpsilon<real>());
1649  eigminG -= smallNumberToExpandWith;
1650  eigmaxG += smallNumberToExpandWith;
1651 #ifdef DEBUG_PURI_OUTPUT
1652  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "EXPANDED Gershgorin bounds: [ %.12lf , %.12lf ]", (double)eigminG, (double)eigmaxG);
1653 #endif
1654 
1655  eigmin = eigminG;
1656  eigmax = eigmaxG;
1657 
1658  // Lanczos helps us to improve Gershgorin bounds
1659 #if 1 // 0 - without Lanczos, 1 - with Lanczos
1660 #ifndef USE_CHUNKS_AND_TASKS
1661 
1662  // try to impove with Lanczos algorithm
1663  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Trying to impove bounds using Lanczos algorithm...");
1664 
1665  real acc = 1e3 * template_blas_sqrt(get_epsilon()); // this accuracy may be too high for single precision
1666  MatrixType Xshifted(X);
1667  Xshifted.add_identity((real)(-1.0) * eigminG); // Xsh = X - eigmin*I
1668 
1669  int maxIter = 100;
1670  try
1671  {
1672  eigmax = Xshifted.eucl(acc, maxIter) + eigminG + acc;
1673  }
1674  catch (mat::AcceptableMaxIter& e)
1675  {
1676 
1677  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Lanczos failed to find extreme upper eigenvalue within maxiter... using Gershgorin bound");
1678 
1679  eigmax = eigmaxG;
1680  }
1681 
1682  // Now we want to create Fshifted = ( F - lambdaMaxGers*id ) but we
1683  // do this starting from the existing Fshifted, correcting it back
1684  // to F and then subtracting lambdaMaxGers*id.
1685  Xshifted.add_identity((real)(1.0) * eigminG); // Now Fshifted = F.
1686  Xshifted.add_identity((real)(-1.0) * eigmaxG);
1687 
1688  try
1689  {
1690  eigmin = -Xshifted.eucl(acc, maxIter) + eigmaxG - acc;
1691  }
1692  catch (mat::AcceptableMaxIter& e)
1693  {
1694 
1695  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Lanczos failed to find extreme lower eigenvalue within maxiter... using Gershgorin bound");
1696 
1697  eigmin = eigminG;
1698  }
1699 #endif // USE_CHUNKS_AND_TASKS
1700 #endif
1701 
1702  // extreme case of matrix with 1 element
1703  if (eigmin == eigmax)
1704  {
1705  real tol = 1e-2;
1706  eigmin -= tol;
1707  eigmax += tol;
1708  }
1709 
1710  spectrum_bounds = IntervalType(eigmin, eigmax);
1711 
1712  computed_spectrum_bounds = true;
1713 }
1714 
1715 
1716 /****************** CHECK_STOPPING_CRITERION **********************/
1717 
1718 template<typename MatrixType>
1720 {
1721  int it = iter_info.it;
1722  real XmX2_norm_it = -1, XmX2_norm_itm2 = -1;
1723 
1724  if (it < check_stopping_criterion_iter)
1725  {
1726  return; // do not check the stopping criterion
1727  }
1728  /* ELIAS NOTE 2018-06-25: it sometimes happens that
1729  check_stopping_criterion_iter is -1 which gives a problem here
1730  when it=1 since in that case we go out-of-bounds when trying to
1731  access Iterations[it - 2] below. To get around this problem I
1732  did an ugly fix by adding the check below, simply doing return
1733  in that case. This should probably be solved in a better way, so
1734  that check_stopping_criterion_iter is always set, hopefully
1735  Anastasia can look at that later. */
1736  /* FIXME ELIAS: ask Anastasia about this! */
1737  if(use_new_stopping_criterion && check_stopping_criterion_iter == -1 && it < 2)
1738  return;
1739  if (use_new_stopping_criterion)
1740  {
1741  assert(it >= 2);
1742  // if spectral norm is used for the etimation of the order
1743  if (normPuriStopCrit == mat::euclNorm)
1744  {
1745  XmX2_norm_it = iter_info.XmX2_eucl;
1746  XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_eucl;
1747  }
1748  else
1749  if (normPuriStopCrit == mat::frobNorm)
1750  {
1751  XmX2_norm_it = iter_info.XmX2_fro_norm;
1752  XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_fro_norm;
1753  }
1754  else
1755  if (normPuriStopCrit == mat::mixedNorm)
1756  {
1757  XmX2_norm_it = iter_info.XmX2_mixed_norm;
1758  XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_mixed_norm;
1759  }
1760  else
1761  {
1762  throw std::runtime_error("Error in stopping_criterion() : unknown matrix norm.");
1763  }
1764 
1765  real XmX2_trace = iter_info.XmX2_trace;
1766  check_new_stopping_criterion(it, XmX2_norm_it, XmX2_norm_itm2, XmX2_trace, stop, estim_order);
1767  }
1768  else // use standard stopping criterion
1769  {
1770  if (normPuriStopCrit == mat::euclNorm)
1771  {
1772  XmX2_norm_it = iter_info.XmX2_eucl;
1773  }
1774  if (normPuriStopCrit == mat::frobNorm)
1775  {
1776  XmX2_norm_it = iter_info.XmX2_fro_norm;
1777  }
1778  if (normPuriStopCrit == mat::mixedNorm)
1779  {
1780  XmX2_norm_it = iter_info.XmX2_mixed_norm;
1781  }
1782 
1783  estim_order = -1;
1784  check_standard_stopping_criterion(XmX2_norm_it, stop);
1785  }
1786 }
1787 
1788 
1789 template<typename MatrixType>
1791 {
1792 
1793  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Checking standard stopping criterion... ");
1794 
1795  bool homoLumo_converged = (homo_bounds.upp() < error_eig &&
1796  lumo_bounds.upp() < error_eig);
1797  bool XmX2norm_converged = XmX2_norm < error_eig;
1798  if (homoLumo_converged || XmX2norm_converged)
1799  {
1800  stop = 1;
1801  }
1802 
1803 #ifdef DEBUG_PURI_OUTPUT
1804  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "stop = %d", stop);
1805 #endif
1806 }
1807 
1808 
1809 template<typename MatrixType>
1810 void PurificationGeneral<MatrixType>::check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace,
1811  int& stop, real& estim_order)
1812 {
1813  // must do at least 2 iterations
1814  if (it < 2)
1815  {
1816  estim_order = -1;
1817  return;
1818  }
1819 
1820 
1821  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Checking stopping criterion... ");
1822 
1823 
1824  real C; // constant may depend on the purification method
1825  return_constant_C(it, C);
1826  this->constant_C = C;
1827  if (C == -1)
1828  {
1829  estim_order = -1;
1830  return;
1831  }
1832 
1833  estim_order = template_blas_log(XmX2_norm_it / C) / template_blas_log(XmX2_norm_itm2); // rate of convergence - due to overflow may be Inf
1834 
1835  if ((VecPoly[it - 1] != VecPoly[it]) && // alternating polynomials
1836  (XmX2_norm_itm2 < 1) && // assumption for frob and mixed norms, always true for the spectral norm
1837  (XmX2_norm_it >= C * template_blas_pow(XmX2_norm_itm2, (real)ORDER))) // r <= 2 (or smaller value)
1838  {
1839  stop = 1;
1840  }
1841 
1842 
1843  if ((stop != 1) && (XmX2_norm_it < get_epsilon() * template_blas_sqrt(template_blas_sqrt(get_epsilon()))))
1844  {
1845  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "************************************************************************************************************");
1846  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "The norm value went much below machine precision, therefore we stop here since n_max can be underestimated.");
1847  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "************************************************************************************************************");
1848  stop = 1;
1849  }
1850 
1851 
1852 
1853 #ifdef DEBUG_PURI_OUTPUT
1854  if ((VecPoly[it - 1] != VecPoly[it]) && (XmX2_norm_itm2 < 1))
1855  {
1856  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "e_i =%e, C*e_{i-1}^q = %e", XmX2_norm_it, C * template_blas_pow(XmX2_norm_itm2, (real)ORDER));
1857  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Order of convergence = %lf, stop = %d", (double)estim_order, stop);
1858  }
1859  else
1860  {
1861  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Order of convergence cannot be computed");
1862  }
1863 #endif
1864 }
1865 
1866 
1867 /********************* TOTAL_SUBSPACE_ERROR ****************/
1868 
1869 template<typename MatrixType>
1872 {
1873  assert(it <= (int)VecGap.size());
1874 
1875  real error = 0;
1876  real normE;
1877  for (int i = 0; i <= it; ++i)
1878  {
1879  if (VecGap[i] == -1)
1880  {
1881  return -1; // gap is not known
1882  }
1883  normE = info.Iterations[i].threshold_X;
1884  error += normE / (VecGap[i] - normE);
1885  }
1886 
1887  return error;
1888 }
1889 
1890 
1891 template<typename MatrixType>
1893 {
1894  return info.estim_total_it;
1895 }
1896 
1897 
1898 template<typename MatrixType>
1900 {
1901  if (info.converged == 1)
1902  {
1903  return info.total_it;
1904  }
1905  else
1906  {
1907  return -1;
1908  }
1909 }
1910 
1911 
1912 /************ GET ESTIMATE OF EIGENVALUES OF F FROM PURIFICATION **************/
1913 template<typename MatrixType>
1915  VectorTypeReal& out_unocc,
1916  VectorTypeReal& out_occ,
1917  int nmax)
1918 {
1919  out_occ.clear();
1920  out_occ.resize(nmax);
1921  out_unocc.clear();
1922  out_unocc.resize(nmax);
1923 
1924  out_unocc[0] = value_unocc;
1925  out_occ[0] = value_occ;
1926  real occ, unocc;
1927 
1928  for (int i = 1; i < nmax; ++i) // note: i starts from 1
1929  {
1930  occ = out_occ[i - 1];
1931  unocc = out_unocc[i - 1];
1932  apply_poly_to_eigs(i, occ, unocc);
1933  out_occ[i] = occ;
1934  out_unocc[i] = unocc;
1935  }
1936 }
1937 
1938 
1939 template<typename MatrixType>
1941  const VectorTypeReal& XmX2_norm_frob,
1942  const VectorTypeReal& XmX2_trace)
1943 {
1944  estimate_homo_lumo(XmX2_norm_mixed, XmX2_norm_frob, XmX2_trace);
1945 
1946  VectorTypeReal h_in, l_in;
1947 
1948  real eigmax = spectrum_bounds.upp();
1949  real eigmin = spectrum_bounds.low();
1950  real homo_in_0 = 1 - (homo_bounds_F_new.low() - eigmax) / (eigmin - eigmax);
1951  real lumo_in_0 = (lumo_bounds_F_new.upp() - eigmax) / (eigmin - eigmax);
1952 
1953  int n = get_exact_number_of_puri_iterations();
1954  assert(n > 0);
1955 
1956  propagate_values_in_each_iter(lumo_in_0, homo_in_0, l_in, h_in, n);
1957 
1958  VectorTypeReal YmY2_norm_frob_est;
1959  get_frob_norm_est(XmX2_norm_frob, h_in, l_in, YmY2_norm_frob_est);
1960 }
1961 
1962 
1963 template<typename MatrixType>
1965  const VectorTypeReal& h_in,
1966  const VectorTypeReal& l_in,
1967  VectorTypeReal& YmY2_norm_frob_est)
1968 {
1969  int n = get_exact_number_of_puri_iterations();
1970 
1971  assert(n > 0);
1972 
1973  YmY2_norm_frob_est.clear();
1974  YmY2_norm_frob_est.resize(n);
1975 
1976  real th, tl;
1977  for (int i = 0; i < n; ++i)
1978  {
1979  th = h_in[i] - h_in[i] * h_in[i];
1980  tl = l_in[i] - l_in[i] * l_in[i];
1981  YmY2_norm_frob_est[i] = template_blas_sqrt(XmX2_norm_frob[i] * XmX2_norm_frob[i] - th * th - tl * tl);
1982  }
1983 }
1984 
1985 
1986 
1987 /*
1988  ||X-X^2||^2_F / trace(X-X^2) <= ||X-X^2||_2 <= ||X-X^2||_m ( <= ||X-X^2||_F),
1989  ||X-X^2||_m - mixed norm of X-X^2
1990  */
1991 template<typename MatrixType>
1993  const VectorTypeReal& XmX2_norm_frob,
1994  const VectorTypeReal& XmX2_trace)
1995 {
1996  homo_bounds_F_new = intervalType(-1e22, 1e22);
1997  lumo_bounds_F_new = intervalType(-1e22, 1e22);
1998 
1999  // do not use additional iterations in the estimation
2000  int total_it = info.total_it - info.additional_iterations;
2001 
2002  // lumo_out, lumo_in, 1-homo_out, 1-homo_in
2003  VectorTypeReal bounds_from_it(4);
2004  VectorTypeReal final_bounds(4, 1); // set all to one
2005 
2006  // criterion for the eligible iterations for the estimation of the bounds
2007  real STOP_NORM = gammaStopEstim - gammaStopEstim * gammaStopEstim;
2008  real vi, wi, mi;
2009  real temp_value;
2010 
2011  VectorTypeReal XmX2_norm_out;
2012  if (XmX2_norm_mixed.size() == XmX2_norm_frob.size())
2013  {
2014  XmX2_norm_out = XmX2_norm_mixed;
2015  }
2016  else
2017  {
2018  XmX2_norm_out = XmX2_norm_frob;
2019  }
2020 
2021  for (int it = total_it; it >= 0; it--)
2022  {
2023  vi = XmX2_norm_frob[it];
2024  wi = XmX2_trace[it];
2025  mi = XmX2_norm_out[it];
2026 
2027  if (vi >= STOP_NORM)
2028  {
2029  break;
2030  }
2031 
2032  if (wi <= template_blas_sqrt(get_epsilon()))
2033  {
2034  continue;
2035  }
2036 
2037  // lumo bounds
2038  temp_value = vi * vi / wi;
2039  bounds_from_it[0] = 0.5 * (1 - template_blas_sqrt(1 - 4 * temp_value));
2040  bounds_from_it[1] = 0.5 * (1 - template_blas_sqrt(1 - 4 * mi));
2041 
2042  // bounds for 1-homo
2043  bounds_from_it[2] = bounds_from_it[0];
2044  bounds_from_it[3] = bounds_from_it[1];
2045 
2046 
2047  apply_inverse_poly_vector(it, bounds_from_it);
2048 
2049 
2050  final_bounds[0] = std::min(final_bounds[0], bounds_from_it[0]); // outer
2051  final_bounds[1] = std::min(final_bounds[1], bounds_from_it[1]); // inner
2052 
2053  final_bounds[2] = std::min(final_bounds[2], bounds_from_it[2]); // outer
2054  final_bounds[3] = std::min(final_bounds[3], bounds_from_it[3]); // inner
2055  }
2056 
2057  // get bounds for F
2058  real maxeig = spectrum_bounds.upp();
2059  real mineig = spectrum_bounds.low();
2060  lumo_bounds_F_new = IntervalType(maxeig * (1 - final_bounds[1]) + mineig * final_bounds[1],
2061  maxeig * (1 - final_bounds[0]) + mineig * final_bounds[0]);
2062  homo_bounds_F_new = IntervalType(mineig * (1 - final_bounds[2]) + maxeig * final_bounds[2],
2063  mineig * (1 - final_bounds[3]) + maxeig * final_bounds[3]);
2064 }
2065 
2066 
2067 /*************************** SAVE MATRIX **************************/
2068 
2069 template<typename MatrixType>
2071 {
2072 #ifdef USE_CHUNKS_AND_TASKS
2073  vector<int> Itmp, I, Jtmp, J;
2074  vector<real> Vtmp, V;
2075  X.get_all_values(Itmp, Jtmp, Vtmp);
2076 
2077  size_t nnz = 0;
2078  // Count nonzeros
2079  for (size_t i = 0; i < Itmp.size(); i++)
2080  {
2081  nnz += (Vtmp[i] != 0);
2082  }
2083 
2084  I.reserve(nnz);
2085  J.reserve(nnz);
2086  V.reserve(nnz);
2087  // Extract nonzeros
2088  for (size_t i = 0; i < Itmp.size(); i++)
2089  {
2090  if (Vtmp[i] != 0)
2091  {
2092  I.push_back(Itmp[i]);
2093  J.push_back(Jtmp[i]);
2094  V.push_back(Vtmp[i]);
2095  }
2096  }
2097 
2098  string name = "X_" + str + ".mtx";
2099  if (write_matrix_to_mtx(name.c_str(), I, J, V, X.get_nrows()) == -1)
2100  {
2101  throw std::runtime_error("Error in save_matrix_now : error in write_matrix_to_mtx.\n");
2102  }
2103 #endif
2104 }
2105 
2106 
2107 /********************************************************
2108  *** COMPUTATION OF EIGENVECTORS ***
2109  *********************************************************/
2110 
2111 
2112 // FUNCTION FOR COMPARISON OF THE RESULTS
2113 #ifdef USE_CHUNKS_AND_TASKS
2114 
2115 template<typename MatrixType>
2117 {
2118  throw std::runtime_error("compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
2119 }
2120 
2121 
2122 template<typename MatrixType>
2124 {
2125  throw std::runtime_error("compute_eigenvectors_without_diagonalization_last_iter_proj() is not implemented for CHTMatrix.");
2126 }
2127 
2128 
2129 template<typename MatrixType>
2131 {
2132  throw std::runtime_error("compute_eigenvectors_without_diagonalization() is not implemented for CHTMatrix.");
2133 }
2134 
2135 
2136 #else
2137 // Assumption: F is not in the file
2138 template<typename MatrixType>
2140 {
2141  real start_shift = homo_bounds_F.upp();
2142  real end_shift = lumo_bounds_F.low();
2143  real dist = (end_shift - start_shift) / 15;
2144  real acc_eigv = eigensolver_accuracy;
2145  real eigval;
2146  int eig_num = 1;
2147 
2148  MatrixType Fsq;
2149 
2150  Fsq = (real)1.0 * F * F; // F^2
2151 
2152  real sigma;
2153 
2154  // choose different shifts
2155  for (sigma = start_shift; sigma < end_shift; sigma += dist)
2156  {
2157  MatrixType M(Fsq);
2158  M.add_identity(sigma * sigma); // F^2 + sigma^2I
2159  M += ((real) - 2 * sigma) * F; // F^2 + sigma^2I - 2sigma*F
2160  M = ((real) - 1.0) * M; // -(F-sigma*I)^2
2161 
2162  vector<real> eigValTmp(1); // here will be computed eigenvalues of M
2163  vector<int> num_iter_solver(1, -1); // number of iterations
2164 
2166  F.getRows(rows);
2167  vector<VectorType> eigVec(1, VectorType(rows)); // here will be computed eigenvectors
2168  eigVec[0].rand(); // initial guess
2169  eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
2170  eigenvectors_iterative_method_str, num_iter_solver,
2171  eigensolver_maxiter_for_F);
2172 
2173  eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2174 
2175  printf("sigma = %lf , eigval = %lf , iters = %d\n", (double)sigma, (double)eigval, num_iter_solver[0]);
2176  }
2177 
2178  sigma = end_shift;
2179  {
2180  MatrixType M(Fsq);
2181  M.add_identity(sigma * sigma); // F^2 + sigma^2I
2182  M += ((real) - 2 * sigma) * F; // F^2 + sigma^2I - 2sigma*F
2183  M = ((real) - 1.0) * M; // -(F-sigma*I)^2
2184 
2185  vector<real> eigValTmp(1); // here will be computed eigenvalues of M
2186  vector<int> num_iter_solver(1, -1); // number of iterations
2187 
2189  F.getRows(rows);
2190  vector<VectorType> eigVec(1, VectorType(rows)); // here will be computed eigenvectors
2191  eigVec[0].rand(); // initial guess
2192  eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
2193  eigenvectors_iterative_method_str, num_iter_solver,
2194  eigensolver_maxiter_for_F);
2195 
2196  eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2197 
2198  printf("sigma = %lf , eigval = %lf , iters = %d\n", (double)sigma, (double)eigval, num_iter_solver[0]);
2199  }
2200 }
2201 
2202 
2203 template<typename MatrixType>
2205 {
2206  real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2207 
2208  /* flags deciding on which iteration to compute eigenvectors */
2209  bool compute_homo_now = false;
2210  bool compute_lumo_now = false;
2211 
2212 
2213  if (compute_eigenvectors_in_each_iteration)
2214  {
2215  if (eigenvectors_method == EIG_SQUARE_INT)
2216  {
2217  // can we compute homo in this iteration?
2218  if (eigVecHOMO != NULL)
2219  {
2220  for (size_t i = 0; i < good_iterations_homo.size(); ++i)
2221  {
2222  if (good_iterations_homo[i] == it)
2223  {
2224  compute_homo_now = true;
2225  break;
2226  }
2227  }
2228  }
2229 
2230  // can we compute lumo in this iteration?
2231  if (eigVecLUMO != NULL)
2232  {
2233  for (size_t i = 0; i < good_iterations_lumo.size(); ++i)
2234  {
2235  if (good_iterations_lumo[i] == it)
2236  {
2237  compute_lumo_now = true;
2238  break;
2239  }
2240  }
2241  }
2242  }
2243  else // eigenvectors_method == EIG_PROJECTION_INT
2244  {
2245  // for projection method we can compute for each X_i, including X_0
2246  if (eigVecHOMO != NULL)
2247  {
2248  compute_homo_now = true;
2249  }
2250  if (eigVecLUMO != NULL)
2251  {
2252  compute_lumo_now = true;
2253  }
2254  }
2255  }
2256  else
2257  {
2258  // compute just in chosen iterations
2259  if (eigVecHOMO != NULL)
2260  {
2261  if (it == iter_for_homo)
2262  {
2263  compute_homo_now = true;
2264  }
2265  }
2266  if (eigVecLUMO != NULL)
2267  {
2268  if (it == iter_for_lumo)
2269  {
2270  compute_lumo_now = true;
2271  }
2272  }
2273  }
2274 
2275  if (compute_homo_now && !info.homo_eigenvector_is_computed)
2276  {
2277  if (eigenvectors_method == EIG_SQUARE_INT)
2278  {
2279  Util::TimeMeter homo_total_time;
2280 
2281  MatrixType M(Xsq); // M = Xsq
2282  writeToTmpFile(Xsq);
2283 
2284  real sigma = SIGMA_HOMO_VEC[it]; // take precomputed shift
2285  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "choose shift %lf", (double)sigma);
2286 
2287  /* GET MATRIX */
2288 
2289  M.add_identity(sigma * sigma); // X^2 + sigma^2I
2290  M += ((real) - 2 * sigma) * X;
2291  M = ((real) - 1.0) * M; // -(X-sigma*I)^2
2292 
2293  Util::TimeMeter homo_solver_time;
2294  compute_eigenvector(M, eigVecHOMO, it, true);
2295  homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2296  homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
2297 
2298  homo_total_time.print(LOG_AREA_DENSFROMF, "compute homo eigenvector");
2299  homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
2300 
2301  iter_info.homo_eig_solver_time = homo_solver_time_stop;
2302  iter_info.orbital_homo_time = homo_total_time_stop;
2303 
2304 
2305  M.clear();
2306  readFromTmpFile(Xsq);
2307  }
2308  else // Projection method
2309  {
2310  if (compute_eigenvectors_in_each_iteration)
2311  {
2312  // Save matrix X_i into the mtx file. We will use it later to compute homo eigenvectors.
2313  ostringstream name;
2314  name << "homo_" << it;
2315  Util::TimeMeter homo_time_save_matrix;
2316  save_matrix_now(name.str());
2317  homo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving homo matrix into mtx");
2318  }
2319  else
2320  {
2321  // Save matrix X_i. We will use it later to compute homo eigenvector.
2322  Util::TimeMeter homo_time_save_matrix;
2323  writeToTmpFile(X);
2324  X_homo = X;
2325  readFromTmpFile(X);
2326  homo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving homo matrix using writeToFile");
2327  }
2328  }
2329  }
2330 
2331 
2332 
2333  if (compute_lumo_now && !info.lumo_eigenvector_is_computed)
2334  {
2335  if (eigenvectors_method == EIG_SQUARE_INT)
2336  {
2337  Util::TimeMeter lumo_total_time;
2338 
2339  MatrixType M(Xsq); // M = Xsq
2340  writeToTmpFile(Xsq);
2341 
2342  real sigma = SIGMA_LUMO_VEC[it]; // take precomputed shift
2343  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "choose shift %lf", (double)sigma);
2344 
2345  /* GET MATRIX */
2346 
2347  M.add_identity(sigma * sigma); // X^2 + sigma^2I
2348  M += ((real) - 2 * sigma) * X;
2349  M = ((real) - 1.0) * M; // -(X-sigma*I)^2
2350 
2351  Util::TimeMeter lumo_solver_time;
2352  compute_eigenvector(M, eigVecLUMO, it, false);
2353  lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2354  lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
2355 
2356  lumo_total_time.print(LOG_AREA_DENSFROMF, "compute lumo eigenvector");
2357  lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
2358 
2359  iter_info.lumo_eig_solver_time = lumo_solver_time_stop;
2360  iter_info.orbital_lumo_time = lumo_total_time_stop;
2361 
2362  M.clear();
2363  readFromTmpFile(Xsq);
2364  }
2365  else // Projection method
2366  {
2367  if (compute_eigenvectors_in_each_iteration)
2368  {
2369  // Save matrix X_i into the mtx file. We will use it later to compute lumo eigenvectors.
2370  ostringstream name;
2371  name << "lumo_" << it;
2372  Util::TimeMeter lumo_time_save_matrix;
2373  save_matrix_now(name.str());
2374  lumo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving lumo matrix into mtx");
2375  }
2376  else
2377  {
2378  // Save matrix X_i. We will use it later to compute lumo eigenvector.
2379  Util::TimeMeter lumo_time_save_matrix;
2380  writeToTmpFile(X);
2381  X_lumo = X;
2382  readFromTmpFile(X);
2383  lumo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving lumo matrix using writeToFile");
2384  }
2385  }
2386  }
2387 
2388 
2389  if (compute_eigenvectors_in_each_iteration)
2390  {
2391  // compute again in the next iteration
2392  info.homo_eigenvector_is_computed = false;
2393  info.lumo_eigenvector_is_computed = false;
2394  }
2395 }
2396 
2397 
2398 template<typename MatrixType>
2400 {
2401  real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2402  real DX_mult_time_homo_stop, DX_mult_time_lumo_stop;
2403 
2404  output_current_memory_usage(LOG_AREA_DENSFROMF, "Before computing eigenvectors:");
2405  Util::TimeMeter timeMeterWriteAndReadAll;
2406  std::string sizesStr = mat::FileWritable::writeAndReadAll();
2407  timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
2408  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
2409 
2410  // we can clear here X^2 !
2411  Xsq.clear();
2412  assert(eigVecHOMO != NULL);
2413  assert(eigVecLUMO != NULL);
2414 
2415  if (compute_eigenvectors_in_each_iteration)
2416  {
2417  int total_it = info.total_it;
2418 
2419  /*
2420  * Since we are computing eigenvectors in every iteration,
2421  * use the same matrix for homo and for lumo
2422  */
2423  for (int it = 0; it <= total_it; ++it)
2424  {
2425  // read matrix from mtx file
2426  MatrixType Xi;
2427  ostringstream name;
2428  name << "X_lumo_" << it << ".mtx";
2431  X.getRows(rows);
2432  X.getCols(cols);
2433  // TODO
2434  // Xi.read_from_mtx(name.str(), rows, cols);
2435 
2436  Util::TimeMeter homo_total_time;
2437 
2438  MatrixType DXi(X);
2439 
2440  Util::TimeMeter DX_mult_time_homo; // separate timing for multiplication D*X
2441  MatrixType TMP(DXi);
2442  // DXi = 1 * TMP * Xi + 0 * DXi
2443  MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, Xi, 0, DXi);
2444  DX_mult_time_homo.print(LOG_AREA_DENSFROMF, "computing D*X");
2445  DX_mult_time_homo_stop = DX_mult_time_homo.get_elapsed_wall_seconds();
2446 
2447  // get HOMO matrix
2448  MatrixType Zh(X);
2449  Zh -= DXi; // D-DXi
2450 
2451  Util::TimeMeter homo_solver_time;
2452  compute_eigenvector(Zh, eigVecHOMO, it, true);
2453  homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2454  homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
2455 
2456  homo_total_time.print(LOG_AREA_DENSFROMF, "computing homo eigenvector");
2457  homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
2458 
2459  Util::TimeMeter lumo_total_time;
2460 
2461  // get LUMO matrix
2462  DXi -= Xi;
2463  DXi = (real)(-1) * DXi; // Xi-DXi
2464 
2465  Util::TimeMeter lumo_solver_time;
2466  compute_eigenvector(DXi, eigVecLUMO, it, false);
2467  lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2468  lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
2469 
2470  lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
2471  lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
2472 
2473  info.Iterations[it].DX_mult_homo_time = DX_mult_time_homo_stop;
2474  info.Iterations[it].DX_mult_lumo_time = 0; // dummy in this case since we reuse matrix
2475 
2476  info.Iterations[it].homo_eig_solver_time = homo_solver_time_stop;
2477  info.Iterations[it].lumo_eig_solver_time = lumo_solver_time_stop;
2478 
2479  info.Iterations[it].orbital_homo_time = homo_total_time_stop;
2480  info.Iterations[it].orbital_lumo_time = lumo_total_time_stop;
2481  }
2482 
2483  output_current_memory_usage(LOG_AREA_DENSFROMF, "After computing eigenvectors:");
2484  return;
2485  }
2486 
2487 
2488  // check if we passed iter_for_homo iteration and saved the matrix
2489  if (info.total_it >= iter_for_homo)
2490  {
2491  // reading X_homo matrix from the bin file
2492  Util::TimeMeter X_homo_read;
2493  readFromTmpFile(X_homo);
2494  X_homo_read.print(LOG_AREA_DENSFROMF, "reading X matrix (for homo) using readFromFile");
2495 
2496 
2497  Util::TimeMeter homo_total_time; // total time for homo
2498 
2499  MatrixType DXi(X);
2500 
2501  // multiplying D*Xi
2502  Util::TimeMeter DX_mult_time_homo;
2503  MatrixType TMP(DXi);
2504  // DXi = 1 * TMP * X_homo + 0 * DXi
2505  MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, X_homo, 0, DXi);
2506  DX_mult_time_homo.print(LOG_AREA_DENSFROMF, "computing D*X (for homo)");
2507  DX_mult_time_homo_stop = DX_mult_time_homo.get_elapsed_wall_seconds();
2508 
2509  // get HOMO matrix
2510  // note: we may need DXi for computing lumo, do not overwrite it
2511  MatrixType Zh(X);
2512  Zh -= DXi; // D-DXi
2513  Util::TimeMeter homo_solver_time;
2514  compute_eigenvector(Zh, eigVecHOMO, iter_for_homo, true);
2515  homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2516  homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
2517 
2518  info.Iterations[iter_for_homo].homo_eig_solver_time = homo_solver_time_stop; // note: here is included just time for compute_eigenvector()
2519  info.Iterations[iter_for_homo].DX_mult_homo_time = DX_mult_time_homo_stop;
2520  Zh.clear();
2521 
2522  homo_total_time.print(LOG_AREA_DENSFROMF, "computing homo eigenvector");
2523  homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
2524  info.Iterations[iter_for_homo].orbital_homo_time = homo_total_time_stop;
2525 
2526 
2527  // if we are computing both homo and lumo in the same iteration
2528  if (iter_for_homo == iter_for_lumo)
2529  {
2530  Util::TimeMeter lumo_total_time;
2531 
2532  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Reuse matrix D*X_i for lumo computations"); // no matrix read
2533 
2534  // get LUMO matrix
2535  DXi -= X_homo;
2536  DXi = (real)(-1) * DXi; // Xi-DXi
2537 
2538  Util::TimeMeter lumo_solver_time;
2539  compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo, false);
2540  lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2541  lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
2542 
2543 
2544  lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
2545  lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
2546 
2547  info.Iterations[iter_for_lumo].DX_mult_lumo_time = 0; // we reuse DX matrix
2548  info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop; // note: here is included just time for compute_eigenvector()
2549  info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
2550  } // LUMO
2551  } // HOMO
2552 
2553  X_homo.clear();
2554 
2555  // check if we passed iter_for_lumo iteration and saved the matrix, and that it was not the same iteration
2556  // as iter_for_homo
2557  if ((info.total_it >= iter_for_lumo) && (iter_for_homo != iter_for_lumo))
2558  {
2559  Util::TimeMeter X_lumo_read;
2560  readFromTmpFile(X_lumo);
2561  X_lumo_read.print(LOG_AREA_DENSFROMF, "reading X matrix (for lumo) using readFromFile");
2562 
2563  Util::TimeMeter lumo_total_time;
2564 
2565  MatrixType DXi(X); // D
2566 
2567  Util::TimeMeter DX_mult_time_lumo;
2568  MatrixType TMP(DXi);
2569  // DXi = 1 * TMP * X_lumo + 0 * DXi
2570  MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, X_lumo, 0, DXi);
2571  DX_mult_time_lumo.print(LOG_AREA_DENSFROMF, "computing D*X (for lumo)");
2572  DX_mult_time_lumo_stop = DX_mult_time_lumo.get_elapsed_wall_seconds();
2573 
2574  // get LUMO matrix
2575  DXi -= X_lumo;
2576  DXi = (real)(-1) * DXi; // Xi-DXi
2577 
2578  Util::TimeMeter lumo_solver_time;
2579  compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo, false);
2580  lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
2581  lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
2582  info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop;
2583  info.Iterations[iter_for_lumo].DX_mult_lumo_time = DX_mult_time_lumo_stop;
2584 
2585  lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
2586  lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
2587  info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
2588 
2589  X_lumo.clear();
2590  } // LUMO
2591 }
2592 
2593 
2594 #endif // USE_CHUNKS_AND_TASKS
2595 
2596 
2597 
2598 template<typename MatrixType>
2600 {
2601  if (eigenvectors_method == EIG_SQUARE_INT)
2602  {
2603  get_iterations_for_lumo_and_homo(iter_for_lumo, iter_for_homo);
2604  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvector for HOMO will be computed on the iteration %d. ", iter_for_homo);
2605  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvector for LUMO will be computed on the iteration %d. ", iter_for_lumo);
2606  }
2607  else if (eigenvectors_method == EIG_PROJECTION_INT)
2608  {
2609  // Compute eigenvectors for X_0.
2610  // Projection method is used just for comparison, thus either eigenvectors
2611  // are computed only on the first iteration or in each iteration.
2612  iter_for_lumo = 0;
2613  iter_for_homo = 0;
2614  }
2615  else
2616  {
2617  throw std::runtime_error("Error in determine_iteration_for_eigenvectors: unknown method for computing eigenvectors.");
2618  }
2619 }
2620 
2621 
2627 template<typename MatrixType>
2629  int& chosen_iter_homo)
2630 {
2631  int maxiter = maxit;
2632  // Inner bounds (from the previous SCF cycle {i-1}, expanded
2633  // with norm of F_i - F_{i-1}).
2634  real homo0 = 1 - homo_bounds.upp(); // inner bounds
2635  real lumo0 = lumo_bounds.upp();
2636  real homoi = homo0, lumoi = lumo0;
2637  real dummy = 0;
2638  real Dh_homo, Dh_lumo, Dgh_homo, Dgh_lumo,
2639  Dgh_homo_max = get_min_double(), Dgh_lumo_max = get_min_double();
2640 
2641  chosen_iter_lumo = -1;
2642  chosen_iter_homo = -1;
2643 
2644  good_iterations_homo.clear();
2645  good_iterations_lumo.clear();
2646 
2647  find_shifts_every_iter();
2648 
2649  for (int i = 1; i <= maxiter; ++i)
2650  {
2651  homoi = apply_poly(i, homoi); // apply POLY[i] on homo
2652  lumoi = apply_poly(i, lumoi); // apply POLY[i] on lumo
2653 
2654  Dh_homo = compute_derivative(i, homo0, dummy);
2655  Dh_lumo = compute_derivative(i, lumo0, dummy);
2656 
2657  // derivative in every iteration
2658  Dgh_homo = 2 * (homoi - SIGMA_HOMO_VEC[i]) * Dh_homo;
2659  Dgh_lumo = 2 * (lumoi - SIGMA_LUMO_VEC[i]) * Dh_lumo;
2660 
2661  if (homoi >= SIGMA_HOMO_VEC[i])
2662  {
2663  good_iterations_homo.push_back(i);
2664  if (Dgh_homo >= Dgh_homo_max)
2665  {
2666  Dgh_homo_max = Dgh_homo;
2667  chosen_iter_homo = i;
2668  }
2669  } // else we cannot be sure which eigenvalue we are computing
2670 
2671  if (lumoi <= SIGMA_LUMO_VEC[i])
2672  {
2673  good_iterations_lumo.push_back(i);
2674  if (template_blas_fabs(Dgh_lumo) >= Dgh_lumo_max) // derivative for lumo is negative
2675  {
2676  Dgh_lumo_max = template_blas_fabs(Dgh_lumo);
2677  chosen_iter_lumo = i;
2678  }
2679  } // else we cannot be sure which eigenvalue we are computing
2680  }
2681 
2682  if ((chosen_iter_homo == -1) || (chosen_iter_lumo == -1))
2683  {
2684  throw "Error in get_iterations_for_lumo_and_homo() : something went wrong, cannot choose iteration to compute HOMO or LUMO eigenvector.";
2685  }
2686 }
2687 
2688 
2689 template<typename MatrixType>
2691 {
2692  int maxiter = this->maxit;
2693 
2694  this->ITER_ERROR_VEC.clear();
2695  this->ITER_ERROR_VEC.resize(maxiter + 1);
2696 
2697  real error = error_per_it;
2698  if (error_per_it == 0)
2699  {
2700  error = this->get_epsilon();
2701  }
2702 
2703  for (int i = 1; i <= maxiter; ++i)
2704  {
2705  this->ITER_ERROR_VEC[i] = (error * this->VecGap[i]) / (1 + error);
2706  }
2707 }
2708 
2709 
2713 template<typename MatrixType>
2715 {
2716  int maxiter = maxit;
2717 
2718  SIGMA_HOMO_VEC.resize(maxiter + 1);
2719  SIGMA_LUMO_VEC.resize(maxiter + 1);
2720 
2721  // Inner bounds (from the previous SCF cycle {i-1}, expanded
2722  // with norm of F_i - F_{i-1}).
2723  real homo = 1 - homo_bounds.upp(); // inner bounds
2724  real lumo = lumo_bounds.upp();
2725  real homo_out = 1 - homo_bounds.low(); // outer bounds
2726  real lumo_out = lumo_bounds.low();
2727 
2728  for (int i = 1; i <= maxiter; ++i)
2729  {
2730  homo = apply_poly(i, homo); // apply POLY[i] on homo
2731  lumo = apply_poly(i, lumo); // apply POLY[i] on lumo
2732 
2733  homo_out = apply_poly(i, homo_out); // apply POLY[i] on homo_out
2734  lumo_out = apply_poly(i, lumo_out); // apply POLY[i] on lumo_out
2735 
2736  SIGMA_HOMO_VEC[i] = (homo_out + lumo) / 2;
2737  SIGMA_LUMO_VEC[i] = (lumo_out + homo) / 2;
2738  }
2739 }
2740 
2741 
2742 template<typename MatrixType>
2744 {
2745  int maxiter = maxit;
2746 
2747  EIG_REL_GAP_HOMO_VEC.resize(maxiter + 1);
2748  EIG_REL_GAP_LUMO_VEC.clear();
2749  EIG_REL_GAP_LUMO_VEC.resize(maxiter + 1);
2750  EIG_ABS_GAP_HOMO_VEC.clear();
2751  EIG_ABS_GAP_HOMO_VEC.resize(maxiter + 1);
2752  EIG_ABS_GAP_LUMO_VEC.clear();
2753  EIG_ABS_GAP_LUMO_VEC.resize(maxiter + 1);
2754 
2755  // Inner bounds (from the previous SCF cycle {i-1}, expanded
2756  // with norm of F_i - F_{i-1}).
2757  real homo0 = homo_bounds_X0.low(); // inner bounds
2758  real lumo0 = lumo_bounds_X0.upp();
2759  real one = 1.0;
2760  real zero = 0.0;
2761 
2762  real homo_map, lumo_map, one_map, zero_map, sigma;
2763 
2764  real homo = homo0, lumo = lumo0;
2765 
2766  for (int i = 1; i <= maxiter; ++i)
2767  {
2768  homo = apply_poly(i, homo); // apply POLY[i] on homo
2769  lumo = apply_poly(i, lumo); // apply POLY[i] on lumo
2770 
2771  sigma = SIGMA_HOMO_VEC[i];
2772 
2773  homo_map = (homo - sigma) * (homo - sigma);
2774  lumo_map = (lumo - sigma) * (lumo - sigma);
2775  one_map = (one - sigma) * (one - sigma);
2776  zero_map = (zero - sigma) * (zero - sigma);
2777 
2778  EIG_ABS_GAP_HOMO_VEC[i] = min(lumo_map - homo_map, one_map - homo_map);
2779  EIG_REL_GAP_HOMO_VEC[i] = EIG_ABS_GAP_HOMO_VEC[i] /
2780  max(zero_map - homo_map, one_map - homo_map);
2781  }
2782 
2783  homo = homo0, lumo = lumo0;
2784 
2785  for (int i = 1; i <= maxiter; ++i)
2786  {
2787  homo = apply_poly(i, homo); // apply POLY[i] on homo
2788  lumo = apply_poly(i, lumo); // apply POLY[i] on lumo
2789 
2790  sigma = SIGMA_LUMO_VEC[i];
2791 
2792  homo_map = (homo - sigma) * (homo - sigma);
2793  lumo_map = (lumo - sigma) * (lumo - sigma);
2794  zero_map = (zero - sigma) * (zero - sigma);
2795  one_map = (one - sigma) * (one - sigma);
2796 
2797  EIG_ABS_GAP_LUMO_VEC[i] = min(homo_map - lumo_map, zero_map - lumo_map);
2798  EIG_REL_GAP_LUMO_VEC[i] = EIG_ABS_GAP_LUMO_VEC[i] /
2799  max(zero_map - lumo_map, one_map - lumo_map);
2800  }
2801 
2802 }
2803 
2804 
2805 #ifdef USE_CHUNKS_AND_TASKS
2806 template<typename MatrixType>
2807 void PurificationGeneral<MatrixType>::compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
2808 {
2809  throw "compute_eigenvector() is not implemented for CHTMatrix.";
2810 }
2811 
2812 
2813 #else
2814 template<typename MatrixType>
2815 void PurificationGeneral<MatrixType>::compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
2816 {
2817  assert(eigVecHOMOorLUMO != NULL);
2818  real acc_eigv = eigensolver_accuracy;
2819 
2820  #ifdef DEBUG_PURI_OUTPUT
2821  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Starting compute_eigenvector()");
2822  #endif
2823 
2824 
2825  if (compute_eigenvectors_in_each_iteration && use_prev_vector_as_initial_guess)
2826  {
2827  // copy vector from previous SCF cycle to be an initial guess
2828  if (is_homo)
2829  {
2830  *eigVecHOMO = eigVecHOMORef;
2831  }
2832  else
2833  {
2834  *eigVecLUMO = eigVecLUMORef;
2835  }
2836  }
2837 
2838 
2839  vector<real> eigValTmp(number_of_eigenvalues); // here will be computed eigenvalues of M
2841 
2842 
2843  X.getRows(rows);
2844 
2845  /*
2846  * Apparently the std::vector constructor calls the copy constructor which is not allowed if the data structure of VectorType is not set.
2847  * In VectorGeneral class were added a new constructor receiving data structure.
2848  */
2849  vector<VectorType> eigVec(number_of_eigenvalues, VectorType(rows)); // here will be computed eigenvectors
2850  if (use_prev_vector_as_initial_guess) // copy vector from previous SCF cycle to be an initial guess
2851  {
2852  use_prev_vector_as_initial_guess = 0; // we firstly check if we have an eigenvector from the previous SCF cycle
2853  if (is_homo)
2854  {
2855  if (!eigVecHOMO->is_empty())
2856  {
2857  eigVec[0] = *eigVecHOMO;
2858  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use HOMO eigenvector computed in the previous SCF cycle as initial guess");
2859  use_prev_vector_as_initial_guess = 1;
2860  }
2861  }
2862  else
2863  {
2864  if (!eigVecLUMO->is_empty())
2865  {
2866  eigVec[0] = *eigVecLUMO;
2867  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use LUMO eigenvector computed in the previous SCF cycle as initial guess");
2868  use_prev_vector_as_initial_guess = 1;
2869  }
2870  }
2871  }
2872  else // use random initial guess
2873  {
2874  eigVec[0].rand();
2875  }
2876 
2877 
2878  vector<int> num_iter_solver(number_of_eigenvalues, -1); // number of iterations
2879 
2880  Util::TimeMeter computeEigenvectors_time;
2881  // run non-deflated version
2882  int eig_num = 1;
2883 
2884  eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
2885  eigenvectors_iterative_method_str, num_iter_solver,
2886  eigensolver_maxiter);
2887  double eigv_elapsed_wall_sec = computeEigenvectors_time.get_elapsed_wall_seconds();
2888  computeEigenvectors_time.print(LOG_AREA_DENSFROMF, "eigensolver");
2889 
2890  if (num_iter_solver.empty())
2891  {
2892  throw std::runtime_error("Error in compute_eigenvector() : (num_iter_solver.empty())");
2893  }
2894 
2895 
2896  // initialize
2897  bool is_homo_tmp = false, is_lumo_tmp = false;
2898 
2899  if (num_iter_solver[0] == eigensolver_maxiter)
2900  {
2901  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigensolver did not converge within maxiter = %d iterations.", eigensolver_maxiter);
2902  }
2903  else
2904  {
2905  is_homo_tmp = is_homo, is_lumo_tmp = !is_homo;
2906  real eigVal;
2907  get_interval_with_lambda(eigVal, eigVec[0], is_homo_tmp, is_lumo_tmp, it); // compute also the corresponding eigenvalue of F
2908  if (is_homo_tmp)
2909  {
2910  really_good_iterations_homo.push_back(it);
2911  *eigVecHOMO = eigVec[0];
2912  info.homo_eigenvector_is_computed = true;
2913  info.homo_eigenvector_is_computed_in_iter = it;
2914  info.homo_eigensolver_iter = num_iter_solver[0];
2915  info.homo_eigensolver_time = eigv_elapsed_wall_sec;
2916  eigValHOMO = eigVal;
2917  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() for HOMO in iteration %d "
2918  ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
2919  }
2920  else if (is_lumo_tmp)
2921  {
2922  really_good_iterations_lumo.push_back(it); // iterations where we computed lumo (for any reason)
2923  *eigVecLUMO = eigVec[0];
2924  info.lumo_eigenvector_is_computed = true;
2925  info.lumo_eigenvector_is_computed_in_iter = it;
2926  info.lumo_eigensolver_iter = num_iter_solver[0];
2927  info.lumo_eigensolver_time = eigv_elapsed_wall_sec;
2928  eigValLUMO = eigVal;
2929  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() for LUMO in iteration %d "
2930  ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
2931  }
2932  else
2933  {
2934  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() in iteration %d "
2935  ": number of eigensolver iterations is %d", it, num_iter_solver[0]);
2936  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error in compute_eigenvector: wrong eigenvalue ");
2937  }
2938  }
2939 
2940  if (compute_eigenvectors_in_each_iteration)
2941  {
2942  save_eigenvectors_to_file(is_homo_tmp, is_lumo_tmp, it);
2943  }
2944  }
2945 
2946 
2947 #endif // USE_CHUNKS_AND_TASKS
2948 
2949 
2950 template<typename MatrixType>
2953 {
2954 #ifndef USE_CHUNKS_AND_TASKS
2955  A.writeToFile();
2956 #endif
2957 }
2958 
2959 
2960 template<typename MatrixType>
2963 {
2964 #ifndef USE_CHUNKS_AND_TASKS
2965  A.readFromFile();
2966 #endif
2967 }
2968 
2969 
2970 template<typename MatrixType>
2972  save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it)
2973 {
2974  if (is_lumo || is_homo)
2975  {
2976  std::vector<real> fullVector;
2977  string eig_name;
2978  if (is_homo)
2979  {
2980  eigVecHOMO->fullvector(fullVector);
2981  eig_name = "homo";
2982  }
2983  else
2984  {
2985  eigVecLUMO->fullvector(fullVector);
2986  eig_name = "lumo";
2987  }
2988 
2989  // save vector to file
2990  ostringstream name;
2991  if (scf_step != -1)
2992  {
2993  name << eig_name << "_" << it << "_scf_step_" << scf_step << ".txt";
2994  }
2995  else
2996  {
2997  name << eig_name << "_" << it << ".txt";
2998  }
2999  if (write_vector(name.str().c_str(), fullVector) == -1)
3000  {
3001  throw std::runtime_error("Error in save_eigenvectors_to_file() : error in write_vector.");
3002  }
3003  name.str("");
3004  }
3005 }
3006 
3007 
3008 // FIXME: add const to vector
3009 template<typename MatrixType>
3011 {
3012  readFromTmpFile(F);
3013  real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3014  writeToTmpFile(F);
3015  eigVal = approx_eig;
3016 }
3017 
3018 
3019 template<typename MatrixType>
3020 void PurificationGeneral<MatrixType>::get_interval_with_lambda(real& eigVal, VectorType& eigVec, bool& is_homo, bool& is_lumo, const int iter)
3021 {
3022  assert(is_homo || is_lumo);
3023 
3024  bool is_homo_init = is_homo;
3025  bool is_lumo_init = is_lumo;
3026 
3027  /* COMPUTE EIGENVALUE OF F CORRESPONDING TO THE COMPUTED EIGENVECTOR USING RAYLEIGH QUOTIENT. */
3028 
3029  /*
3030  * Note: For the square method we compute eigenvalues in the current
3031  * iteration during the purification. The bounds lumo_bounds and
3032  * homo_bounds are changing in every iteration according to the
3033  * polynomial (without expansion by tau). Thus we should use this
3034  * bounds if square method is used. However, for the projection
3035  * method we should used bounds which are saved into the info object.
3036  */
3037  real low_lumo_F_bound, low_homo_F_bound;
3038  real upp_lumo_F_bound, upp_homo_F_bound;
3039 
3040  if (eigenvectors_method == EIG_SQUARE_INT)
3041  // for square method use bounds from the previous SCF cycle (i.e. bounds expanded with norm ||F_i-F_{i-1}||)
3042  {
3043  low_lumo_F_bound = lumo_bounds_F.low();
3044  upp_lumo_F_bound = lumo_bounds_F.upp();
3045  low_homo_F_bound = homo_bounds_F.low();
3046  upp_homo_F_bound = homo_bounds_F.upp();
3047  }
3048  else if (eigenvectors_method == EIG_PROJECTION_INT)
3049  // for projection method we can use new bounds computed in this SCF cycle
3050  {
3051  low_lumo_F_bound = info.lumo_estim_low_F;
3052  upp_lumo_F_bound = info.lumo_estim_upp_F;
3053  low_homo_F_bound = info.homo_estim_low_F;
3054  upp_homo_F_bound = info.homo_estim_upp_F;
3055  }
3056  else
3057  {
3058  throw std::runtime_error("Error in get_interval_with_lambda() : unexpected eigenvectors_method value.");
3059  }
3060 
3061 #ifdef DEBUG_PURI_OUTPUT
3062  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Check Rayleigh quotient...");
3063 #endif
3064 
3065  readFromTmpFile(F);
3066  real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3067  writeToTmpFile(F);
3068 
3069  eigVal = approx_eig;
3070 
3071  real flex_tolerance = THRESHOLD_EIG_TOLERANCE;
3072 
3073  // it is HOMO
3074  if ((approx_eig <= upp_homo_F_bound + flex_tolerance) && (approx_eig >= low_homo_F_bound - flex_tolerance))
3075  {
3076  is_homo = true;
3077  is_lumo = false;
3078 
3079  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO eigenvalue of F is %lf, "
3080  "HOMO bounds are [ %lf , %lf ]", (double)approx_eig, (double)low_homo_F_bound, (double)upp_homo_F_bound);
3081 
3082  iter_for_homo = iter; // We already computed homo in this iteration (in case we thought that it was lumo)
3083  // Do we want to recompute lumo in the next iteration?
3084  if (is_lumo_init && (eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3085  {
3086  iter_for_lumo = iter_for_lumo + 1;
3087  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute LUMO in the next iteration.");
3088  }
3089  }
3090  // it is LUMO
3091  else if ((approx_eig <= upp_lumo_F_bound + flex_tolerance) && (approx_eig >= low_lumo_F_bound - flex_tolerance))
3092  {
3093  is_homo = false;
3094  is_lumo = true;
3095 
3096  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed LUMO eigenvalue of F is %lf, "
3097  "LUMO interval [ %lf , %lf ]", (double)approx_eig, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
3098 
3099  iter_for_lumo = iter; // We already computed lumo in this iteration (in case we thought that it was homo)
3100  // Do we want to recompute homo in the next iteration?
3101  if (is_homo_init && (eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3102  {
3103  iter_for_homo = iter_for_homo + 1;
3104  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute HOMO in the next iteration.");
3105  }
3106  }
3107  else
3108  {
3109  is_homo = false;
3110  is_lumo = false;
3111 
3112  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvalue is outside of both intervals for homo and lumo.");
3113  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvalue is %lf, HOMO interval [ %lf , %lf ], LUMO interval [ %lf , %lf ]",
3114  (double)approx_eig, (double)low_homo_F_bound, (double)upp_homo_F_bound, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
3115 
3116  // Do we want to recompute homo (or lumo) in the next iteration?
3117  // We will try to compute HOMO, however, it can be LUMO.
3118  // If it will be LUMO, we save it as computed LUMO and continue with computing HOMO in the next iteration.
3119  if ((eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3120  {
3121  iter_for_homo = iter_for_homo + 1;
3122  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute HOMO (or LUMO) in the next iteration.");
3123  }
3124  }
3125 }
3126 
3127 
3128 
3129 /******************************************************************/
3130 /*********************** MATLAB FUNCTIONS *************************/
3131 /******************************************************************/
3132 
3133 template<typename MatrixType>
3135 {
3136  std::ofstream f;
3137  f.open(filename, std::ios::out);
3138  if (!f.is_open())
3139  {
3140  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3141  return;
3142  }
3143 
3144  int it = info.total_it;
3145  // save POLY
3146  f << "POLY = [";
3147  for (int i = 0; i <= it; ++i)
3148  {
3149  f << VecPoly[i] << " ";
3150  }
3151  f << "];" << std::endl;
3152 
3153  // choose norm
3154  if (normPuriStopCrit == mat::euclNorm)
3155  {
3156  f << "X_norm = [";
3157  for (int i = 0; i <= it; ++i) // including 0 iteration = original matrix X
3158  {
3159  f << (double)info.Iterations[i].XmX2_eucl << " ";
3160  }
3161  f << "];" << std::endl;
3162  f << " norm_letter = '2';" << std::endl;
3163  }
3164  else if (normPuriStopCrit == mat::frobNorm)
3165  {
3166  f << "X_norm = [";
3167  for (int i = 0; i <= it; ++i) // including 0 iteration = original matrix X
3168  {
3169  f << (double)info.Iterations[i].XmX2_fro_norm << " ";
3170  }
3171  f << "];" << std::endl;
3172  f << " norm_letter = 'F';" << std::endl;
3173  }
3174  else if (normPuriStopCrit == mat::mixedNorm)
3175  {
3176  f << "X_norm = [";
3177  for (int i = 0; i <= it; ++i) // including 0 iteration = original matrix X
3178  {
3179  f << (double)info.Iterations[i].XmX2_mixed_norm << " ";
3180  }
3181  f << "];" << std::endl;
3182  f << " norm_letter = 'M';" << std::endl;
3183  }
3184  else
3185  {
3186  throw "Wrong norm in PurificationGeneral::gen_matlab_file_norm_diff()";
3187  }
3188 
3189  f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
3190  f << "it = " << it << ";" << std::endl;
3191  f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3192  f << "fighandle = figure; clf;" << std::endl;
3193  f << "MARKER = ['o', '+'];" << std::endl;
3194  f << "semilogy(0:stop_iteration, X_norm(1:stop_iteration+1), '-b', plot_props{:});" << std::endl;
3195  f << "hold on" << std::endl;
3196  f << "for i = 1:stop_iteration+1" << std::endl;
3197  f << " if POLY(i) == 1" << std::endl;
3198  f << " h1 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3199  f << " else" << std::endl;
3200  f << " h2 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3201  f << " end" << std::endl;
3202  f << "end" << std::endl;
3203  f << "if stop_iteration ~= it" << std::endl;
3204  f << "h3 = semilogy(stop_iteration+1:it, X_norm(stop_iteration+2:it+1), '-.vr', plot_props{:});" << std::endl;
3205  f << "semilogy(stop_iteration:stop_iteration+1, X_norm(stop_iteration+1:stop_iteration+2), '-.r', plot_props{:});" << std::endl;
3206  f << "legend([h1 h2 h3],{'$x^2$', '$2x-x^2$', 'After stop'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3207  f << "else" << std::endl;
3208  f << "legend([h1 h2],{'$x^2$', '$2x-x^2$'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3209  f << "end" << std::endl;
3210  f << "xlabel('Iteration SP2', 'Interpreter','latex');" << std::endl;
3211  f << "ylabel({['$\\|X_i-X_i^2\\|_{' norm_letter '}$']},'interpreter','latex');" << std::endl;
3212  f << "grid on" << std::endl;
3213  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3214  f << "set(gca,'FontSize',20);" << std::endl;
3215  f << "xlim([0 it]);" << std::endl;
3216  f << "ylim([-inf inf]);" << std::endl;
3217  f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
3218  f << "a = 16; S = logspace(-a, 1, a+2);" << std::endl;
3219  f << "set(gca,'YTick',S(1:2:end));" << std::endl;
3220 
3221  f << "hold off" << std::endl;
3222 
3223  f << "% print( fighandle, '-depsc2', 'norm_diff.eps');" << std::endl;
3224 
3225  f.close();
3226 }
3227 
3228 
3229 template<typename MatrixType>
3231 {
3232  std::ofstream f;
3233  f.open(filename, std::ios::out);
3234  if (!f.is_open())
3235  {
3236  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3237  return;
3238  }
3239 
3240  int it = info.total_it;
3241  f << "Thresh = [";
3242 
3243  for (int i = 0; i <= it; ++i)
3244  {
3245  f << (double)info.Iterations[i].threshold_X << " ";
3246  }
3247  f << "];" << std::endl;
3248 
3249  f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
3250  f << "it = " << it << ";" << std::endl;
3251  f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3252  f << "fighandle = figure; clf;" << std::endl;
3253  f << "semilogy(0:stop_iteration, Thresh(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3254  f << "hold on" << std::endl;
3255  f << "if stop_iteration ~= it" << std::endl;
3256  f << "semilogy(stop_iteration+1:it, Thresh(stop_iteration+2:it+1), '-^r', plot_props{:});" << std::endl;
3257  f << "semilogy(stop_iteration:stop_iteration+1, Thresh(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3258  f << "legend('before stop', 'after stop', 'Location','NorthWest');" << std::endl;
3259  f << "end" << std::endl;
3260  f << "grid on" << std::endl;
3261  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3262  f << "set(gca,'FontSize',20);" << std::endl;
3263  f << "xlim([0 it]);" << std::endl;
3264  f << "ylim([-inf inf]);" << std::endl;
3265  f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
3266  f << "hold off" << std::endl;
3267  f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3268  f << "ylabel('Threshold value', 'interpreter','latex');" << std::endl;
3269  f << "% print( fighandle, '-depsc2', 'threshold.eps');" << std::endl;
3270 
3271  f.close();
3272 }
3273 
3274 
3275 template<typename MatrixType>
3277 {
3278  std::ofstream f;
3279  f.open(filename, std::ios::out);
3280  if (!f.is_open())
3281  {
3282  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3283  return;
3284  }
3285 
3286  int it = info.total_it;
3287  f << "NNZ_X = [";
3288 
3289  for (int i = 0; i <= it; ++i)
3290  {
3291  f << (double)info.Iterations[i].NNZ_X << " ";
3292  }
3293  f << "];" << std::endl;
3294 
3295  f << "NNZ_X2 = [";
3296 
3297  for (int i = 0; i <= it; ++i)
3298  {
3299  f << (double)info.Iterations[i].NNZ_X2 << " ";
3300  }
3301  f << "];" << std::endl;
3302 
3303 
3304 
3305  f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
3306  f << "it = " << it << ";" << std::endl;
3307  f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3308  f << "fighandle = figure; clf;" << std::endl;
3309  f << "h2 = plot(0:stop_iteration, NNZ_X(1:stop_iteration+1), '-ob', plot_props{:});" << std::endl;
3310  f << "hold on" << std::endl;
3311  f << "h1 = plot(0:stop_iteration, NNZ_X2(1:stop_iteration+1), '-vm', plot_props{:});" << std::endl;
3312  f << "if stop_iteration ~= it" << std::endl;
3313  f << "plot(stop_iteration+1:it, NNZ_X(stop_iteration+2:it+1), '-vr', plot_props{:});" << std::endl;
3314  f << "plot(stop_iteration+1:it, NNZ_X2(stop_iteration+2:it+1), '-*r', plot_props{:});" << std::endl;
3315  f << "plot(stop_iteration:stop_iteration+1, NNZ_X(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3316  f << "plot(stop_iteration:stop_iteration+1, NNZ_X2(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3317  f << "end" << std::endl;
3318  f << "legend([h1, h2], {'$X^2$', '$X$'}, 'interpreter','latex');" << std::endl;
3319  f << "grid on" << std::endl;
3320  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3321  f << "set(gca,'FontSize',20);" << std::endl;
3322  f << "xlim([0 it]);" << std::endl;
3323  f << "ylim([0 inf]);" << std::endl;
3324  f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
3325  f << "hold off" << std::endl;
3326  f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3327  f << "ylabel('NNZ [\\%]', 'interpreter','latex');" << std::endl;
3328  f << "% print( fighandle, '-depsc2', 'nnz.eps');" << std::endl;
3329 
3330  f.close();
3331 }
3332 
3333 
3334 template<typename MatrixType>
3336 {
3337  std::ofstream f;
3338  f.open(filename, std::ios::out);
3339  if (!f.is_open())
3340  {
3341  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3342  return;
3343  }
3344 
3345  int it = info.total_it;
3346  f << "in= [";
3347 
3348  for (int i = 0; i <= it; ++i)
3349  {
3350  f << (double)(1 - info.Iterations[i].homo_bound_upp - info.Iterations[i].lumo_bound_upp) << " ";
3351  }
3352  f << "];" << std::endl;
3353 
3354  f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
3355  f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3356  f << "fighandle = figure; clf;" << std::endl;
3357  f << "plot(0:stop_iteration, 1./in(1:stop_iteration+1), '-*r', plot_props{:});" << std::endl;
3358  f << "grid on" << std::endl;
3359  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3360  f << "set(gca,'FontSize',20);" << std::endl;
3361  f << "xlim([0 it]);" << std::endl;
3362  f << "ylim([-inf inf]);" << std::endl;
3363  f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
3364  f << "hold off" << std::endl;
3365  f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3366  f << "ylabel('$\\kappa$', 'interpreter','latex');" << std::endl;
3367  f << "% print( fighandle, '-depsc2', 'cond.eps');" << std::endl;
3368 
3369  f.close();
3370 }
3371 
3372 
3373 template<typename MatrixType>
3375 {
3376  std::ofstream f;
3377  f.open(filename, std::ios::out);
3378  if (!f.is_open())
3379  {
3380  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3381  return;
3382  }
3383 
3384  int it = info.total_it;
3385  f << "homo_low= [";
3386 
3387  for (int i = 0; i <= it; ++i)
3388  {
3389  f << (double)info.Iterations[i].homo_bound_low << " ";
3390  }
3391  f << "];" << std::endl;
3392 
3393  f << "homo_upp= [";
3394 
3395  for (int i = 0; i <= it; ++i)
3396  {
3397  f << (double)info.Iterations[i].homo_bound_upp << " ";
3398  }
3399  f << "];" << std::endl;
3400 
3401  f << "lumo_low= [";
3402 
3403  for (int i = 0; i <= it; ++i)
3404  {
3405  f << (double)info.Iterations[i].lumo_bound_low << " ";
3406  }
3407  f << "];" << std::endl;
3408 
3409  f << "lumo_upp= [";
3410 
3411  for (int i = 0; i <= it; ++i)
3412  {
3413  f << (double)info.Iterations[i].lumo_bound_upp << " ";
3414  }
3415  f << "];" << std::endl;
3416 
3417 
3418 
3419  f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
3420  f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3421  f << "fighandle = figure; clf;" << std::endl;
3422  f << "semilogy(0:stop_iteration, homo_upp(1:stop_iteration+1), '-^b', plot_props{:});" << std::endl;
3423  f << "hold on" << std::endl;
3424  f << "semilogy(0:stop_iteration, homo_low(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3425  f << "semilogy(0:stop_iteration, lumo_low(1:stop_iteration+1), '-vr', plot_props{:});" << std::endl;
3426  f << "semilogy(0:stop_iteration, lumo_upp(1:stop_iteration+1), '-^r', plot_props{:});" << std::endl;
3427  f << "grid on" << std::endl;
3428  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3429  f << "set(gca,'FontSize',20);" << std::endl;
3430  f << "xlim([0 stop_iteration]);" << std::endl;
3431  f << "set(gca,'XTick',[0 5:5:stop_iteration]);" << std::endl;
3432  f << "ylim([-inf inf]);" << std::endl;
3433  f << "hold off" << std::endl;
3434  f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3435  f << "ylabel('Homo and lumo bounds', 'interpreter','latex');" << std::endl;
3436  f << "% print( fighandle, '-depsc2', 'eigs.eps');" << std::endl;
3437 
3438  f.close();
3439 }
3440 
3441 
3442 template<typename MatrixType>
3444 {
3445  std::ofstream f;
3446  f.open(filename, std::ios::out);
3447  if (!f.is_open())
3448  {
3449  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3450  return;
3451  }
3452 
3453  int it = info.total_it;
3454 
3455  f << "time_total = [";
3456  for (int i = 0; i <= it; ++i)
3457  {
3458  f << std::setprecision(16) << (double)info.Iterations[i].total_time << " ";
3459  }
3460  f << "];" << std::endl;
3461 
3462  f << "time_square = [";
3463  for (int i = 0; i <= it; ++i)
3464  {
3465  f << std::setprecision(16) << (double)info.Iterations[i].Xsquare_time << " ";
3466  }
3467  f << "];" << std::endl;
3468 
3469  f << "time_trunc = [";
3470  for (int i = 0; i <= it; ++i)
3471  {
3472  f << std::setprecision(16) << (double)info.Iterations[i].trunc_time << " ";
3473  }
3474  f << "];" << std::endl;
3475 
3476  if (info.compute_eigenvectors_in_this_SCF_cycle)
3477  {
3478  f << "time_eigenvectors_homo = [";
3479  for (int i = 0; i <= it; ++i)
3480  {
3481  f << std::setprecision(16) << (double)info.Iterations[i].orbital_homo_time << " ";
3482  }
3483  f << "];" << std::endl;
3484  f << "time_eigenvectors_lumo = [";
3485  for (int i = 0; i <= it; ++i)
3486  {
3487  f << std::setprecision(16) << (double)info.Iterations[i].orbital_lumo_time << " ";
3488  }
3489  f << "];" << std::endl;
3490 
3491  f << "time_solver_homo = [";
3492  for (int i = 0; i <= it; ++i)
3493  {
3494  f << std::setprecision(16) << (double)info.Iterations[i].homo_eig_solver_time << " ";
3495  }
3496  f << "];" << std::endl;
3497  f << "time_solver_lumo = [";
3498  for (int i = 0; i <= it; ++i)
3499  {
3500  f << std::setprecision(16) << (double)info.Iterations[i].lumo_eig_solver_time << " ";
3501  }
3502  f << "];" << std::endl;
3503 
3504 
3505  if (eigenvectors_method == EIG_SQUARE_INT)
3506  {
3507  // time for X^2, truncation, eigensolver for homo and lumo, additional time for homo and lumo and other time
3508  f << "X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_eigenvectors_homo-time_solver_homo;"
3509  " time_eigenvectors_lumo-time_solver_lumo; "
3510  " time_total - time_square - time_trunc - time_eigenvectors_homo - time_eigenvectors_lumo];" << std::endl;
3511  }
3512  else
3513  {
3514  f << "time_DX_homo = [";
3515  for (int i = 0; i <= it; ++i)
3516  {
3517  f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_homo_time << " ";
3518  }
3519  f << "];" << std::endl;
3520  f << "time_DX_lumo = [";
3521  for (int i = 0; i <= it; ++i)
3522  {
3523  f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_lumo_time << " ";
3524  }
3525  f << "];" << std::endl;
3526 
3527  // for projection total time of the iteration does not include computation of homo and lumo
3528  // time for X^2, eigensolver for homo and lumo, DX multiplication
3529  f << "X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_DX_homo; time_DX_lumo;"
3530  " time_eigenvectors_homo - time_DX_homo - time_solver_homo; time_eigenvectors_lumo - time_DX_lumo - time_solver_lumo;"
3531  " time_total - time_square - time_trunc];" << std::endl;
3532  }
3533  }
3534  else
3535  {
3536  f << "X = [time_square; time_trunc; time_total - time_square - time_trunc];" << std::endl;
3537  }
3538 
3539  f << "it = " << it << ";" << std::endl;
3540  f << "xtick = 0:it;" << std::endl;
3541  f << "fighandle = figure; clf;" << std::endl;
3542  f << "b=bar(xtick, X', 'stacked');" << std::endl;
3543  f << "axis tight;" << std::endl;
3544  f << "grid on" << std::endl;
3545  f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3546  f << "set(gca,'FontSize',20);" << std::endl;
3547  f << "xlim([0 it]);" << std::endl;
3548  f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
3549  f << "ylim([-inf inf]);" << std::endl;
3550  f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3551  f << "ylabel('Time [s]', 'interpreter','latex');" << std::endl;
3552  if (info.compute_eigenvectors_in_this_SCF_cycle)
3553  {
3554 /*
3555  * Legend with matlab's bar appear with default settings in reverse or order. Thus we force it to follow the order of the bar manually.
3556  */
3557  if (eigenvectors_method == EIG_SQUARE_INT)
3558  {
3559  f << "legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3560  }
3561  else
3562  {
3563  f << "legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'DX (lumo)', 'DX (homo)', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3564  }
3565  }
3566  else
3567  {
3568  f << "legend(flipud(b(:)), {'other', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3569  }
3570  f << "% print( fighandle, '-depsc2', 'time.eps');" << std::endl;
3571 
3572  f.close();
3573 }
3574 
3575 
3576 
3577 /******************************************************************/
3578 /*********************** PYTHON FUNCTIONS *************************/
3579 /******************************************************************/
3580 
3581 template<typename MatrixType>
3583 {
3584  std::ofstream f;
3585  f.open(filename, std::ios::out);
3586  if (!f.is_open())
3587  {
3588  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
3589  return;
3590  }
3591 
3592  f << "import numpy as np" << std::endl;
3593  f << "import pylab as pl" << std::endl;
3594  f << "import matplotlib.font_manager as font_manager" << std::endl;
3595  f << "from matplotlib import rc" << std::endl;
3596  f << "rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})" << std::endl;
3597  f << "rc('text', usetex=True)" << std::endl;
3598  f << std::endl;
3599 
3600  int it = info.total_it;
3601  f << "NNZ_X = np.array([";
3602 
3603  for (int i = 0; i <= it; ++i)
3604  {
3605  f << (double)info.Iterations[i].NNZ_X << ", ";
3606  }
3607  f << "])" << std::endl;
3608 
3609  f << "NNZ_X2 = np.array([";
3610 
3611  for (int i = 0; i <= it; ++i)
3612  {
3613  f << (double)info.Iterations[i].NNZ_X2 << ", ";
3614  }
3615  f << "])" << std::endl;
3616 
3617  f << "stop_iteration = " << it - info.additional_iterations << std::endl;
3618  f << "it = " << it << std::endl;
3619  f << "font_prop = font_manager.FontProperties(size=20)" << std::endl;
3620  f << "prop = {'markersize':8, 'fillstyle':'none', 'linewidth':2, 'markeredgewidth':2}" << std::endl;
3621  f << "fig1 = pl.figure(figsize = (8, 6), num='nnz')" << std::endl;
3622  f << "p1, = pl.plot(range(0, stop_iteration+1), NNZ_X2[0:stop_iteration+1], '-vm', **prop);" << std::endl;
3623  f << "p2, = pl.plot(range(0, stop_iteration+1), NNZ_X[0:stop_iteration+1], '-ob', **prop);" << std::endl;
3624  f << "if stop_iteration != it:" << std::endl;
3625  f << " pl.plot(range(stop_iteration+1,it+1), NNZ_X[stop_iteration+1:it+1], '-vr', **prop);" << std::endl;
3626  f << " pl.plot(range(stop_iteration+1,it+1), NNZ_X2[stop_iteration+1:it+1], '-*r', **prop);" << std::endl;
3627  f << " pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
3628  f << " pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X2[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
3629  f << std::endl;
3630  f << "pl.xlim(0, it);" << std::endl;
3631  f << "pl.ylim(ymin=0);" << std::endl;
3632  f << "pl.xlabel('Iteration SP2', fontproperties=font_prop);" << std::endl;
3633  f << "pl.ylabel('NNZ [\\%]', fontproperties=font_prop); " << std::endl;
3634  f << "pl.legend((p1, p2), ('$X^2$', '$X$'), loc='best', numpoints=1)" << std::endl;
3635  f << "pl.xticks(range(5, it, 5))" << std::endl;
3636  f << "pl.tick_params(labelsize=20)" << std::endl;
3637  f << "pl.grid(which='major', alpha=0.5, color='k', linestyle='--', linewidth=1) " << std::endl;
3638  f << "pl.savefig('nnz.eps', format='eps', dpi=1000)" << std::endl;
3639  f << "pl.show()" << std::endl;
3640 
3641  f.close();
3642 }
3643 
3644 
3645 #endif //HEADER_PURIFICATION_GENERAL
int EIG_LANCZOS_INT
Definition: purification_general.cc:57
real error_sub
Allowed error in invariant subspaces.
Definition: purification_general.h:419
real error_per_it
Error allowed in each iteration due to truncation.
Definition: purification_general.h:422
Treal upp() const
Definition: Interval.h:145
#define A
void set_compute_eigenvectors_in_each_iteration()
Definition: purification_general.h:239
VectorTypeInt really_good_iterations_homo
Iterations where homo eigenvector is actually computed.
Definition: purification_general.h:500
real total_time
Definition: puri_info.h:61
real lumo_bound_upp
Definition: puri_info.h:87
int eigensolver_maxiter
Definition: purification_general.h:474
VectorTypeReal EIG_ABS_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for lumo eigenvalue.
Definition: purification_general.h:441
#define LOG_CAT_WARNING
Definition: output.h:48
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)=0
Treal template_blas_pow(Treal x, Treal y)
PurificationGeneral()
Definition: purification_general.h:122
virtual bool is_initialized() const
Check is function initialize() is already called.
Definition: purification_general.h:152
double ergo_real
Definition: realtype.h:69
VectorTypeReal ITER_ERROR_VEC
(Eigenvectors) Maximum error introduced in each iteration.
Definition: purification_general.h:439
ergo_real real
Definition: purification_general.h:114
IntervalType spectrum_bounds
Outer bounds for the whole spectrum of F/Xi.
Definition: purification_general.h:462
mat::SizesAndBlocks rows
Definition: test.cc:51
bool computed_spectrum_bounds
Definition: purification_general.h:463
Functionality for writing output messages to a text file.
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition: purification_general.h:436
int use_new_stopping_criterion
True for new stopping criterion.
Definition: purification_general.h:401
int computeEigenvectors(const MatrixType &A, Treal tol, std::vector< Treal > &eigVal, std::vector< VectorType > &eigVec, int number_of_eigenvalues_to_compute, std::string method, std::vector< int > &num_iter, int maxit=200, bool do_deflation=false)
Function for choosing method for computing eigenvectors.
Definition: get_eigenvectors.h:237
mat::SizesAndBlocks cols
Definition: test.cc:52
virtual void set_init_params()=0
real homo_bound_upp
Definition: puri_info.h:85
real lumo_bound_low
Definition: puri_info.h:86
void enable_printf_output()
Definition: output.cc:66
void output_norms_and_traces(IterationInfo &iter_info) const
Definition: purification_general.h:785
intervalType IntervalType
Definition: random_matrices.h:67
Defined namespace eigvec containing functions for computing largest eigenvalues and corresponding eig...
VectorTypeReal EIG_ABS_GAP_LUMO_VEC
Definition: purification_general.h:441
Wrapper routines for different parts of the integral code, including conversion of matrices from/to t...
Treal low() const
Definition: Interval.h:144
int write_vector(const char *filename, const vector< real > &v)
Definition: files_dense.cc:206
virtual void stopping_criterion(IterationInfo &iter_info, int &stop, real &estim_order)
Choose stopping criterion and check it.
Definition: purification_general.h:1719
int use_prev_vector_as_initial_guess
Definition: purification_general.h:479
ergo_real real
Definition: purification_general.h:64
Basic OS access utilities.
Definition: Failure.h:81
virtual void purify_bounds(const int it)=0
IntervalType lumo_bounds_X0
Initial lumo bounds for X.
Definition: purification_general.h:453
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:111
real eigensolver_accuracy
Definition: purification_general.h:472
MatrixType X_homo
Definition: purification_general.h:257
real TOL_OVERLAPPING_BOUNDS
If the difference between inner bounds for homo and lumo is less then this tolerance, i.e.
Definition: purification_general.cc:45
void get_eigenvalue_of_F_from_eigv_of_Xi(real &eigVal, const VectorType &eigVec)
Definition: purification_general.h:3010
static real get_epsilon()
Get machine epsilon.
Definition: purification_general.h:198
void find_truncation_thresh_every_iter()
Definition: purification_general.h:2690
VectorTypeInt good_iterations_homo
Iterations where homo eigenvector can be computed.
Definition: purification_general.h:497
mat::Interval< ergo_real > intervalType
Definition: matrix_typedefs.h:78
Proxy structs used by the matrix API.
void get_interval_with_lambda(real &eigVal, VectorType &eigVec, bool &is_homo, bool &is_lumo, const int iter)
Definition: purification_general.h:3020
virtual void check_standard_stopping_criterion(const real XmX2_norm, int &stop)
Check stopping criterion (obsolete).
Definition: purification_general.h:1790
VectorTypeInt good_iterations_lumo
Iterations where homo eigenvector can be computed.
Definition: purification_general.h:498
string eigenvectors_iterative_method_str
Definition: purification_general.h:477
File contataining definitions of constants used in the recursive expansion.
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
void get_iterations_for_lumo_and_homo(int &chosen_iter_lumo, int &chosen_iter_homo)
Find the best iterations for computing eigenvectors.
Definition: purification_general.h:2628
IntervalType homo_bounds
(1-homo) bounds for Xi in iteration i
Definition: purification_general.h:449
IntervalType lumo_bounds
Lumo bounds for Xi in iteration i.
Definition: purification_general.h:450
int get_int_eig_iter_method(string eigenvectors_iterative_method)
Definition: purification_general.h:766
File containing declaration of functions for reading/writing dense matrices and vectors.
real constant_C
Asymptotic constant C needed for the new stopping criterion.
Definition: purification_general.h:424
real homo_bound_low
Definition: puri_info.h:84
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:405
real XmX2_mixed_norm
Definition: puri_info.h:75
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
generalVector VectorType
Definition: purification_general.h:119
MatrixType X_lumo
Save matrix Xi in certain iterations for computing eigenvectors (projection method).
Definition: purification_general.h:257
#define LOG_AREA_DENSFROMF
Definition: output.h:61
virtual ~PurificationGeneral()
Definition: purification_general.h:171
real THRESHOLD_EIG_TOLERANCE
Inner homo and lumo bounds may be too good, and it may happen that computed eigenvalue slightly outsi...
Definition: purification_general.cc:49
void gen_matlab_file_threshold(const char *filename) const
Create MATLAB .m file which plots the actual introduced error after truncation of the matrix X_i in e...
Definition: purification_general.h:3230
virtual void initialize(const MatrixType &F_, const IntervalType &lumo_bounds_, const IntervalType &homo_bounds_, int maxit_, real error_sub_, real error_eig_, int use_new_stopping_criterion_, NormType norm_truncation, NormType norm_stop_crit, int nocc_)
Set imporatant parameters for the recursive expansion.
Definition: purification_general.h:521
MatrixType MatrixTypeWrapper
Definition: purification_general.h:120
std::vector< int > VectorTypeInt
Definition: purification_general.h:117
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
File containing declaration of functions for reading/writing sparse matrices from/to binary files...
int get_int_eig_method(string eigenvectors_method)
Definition: purification_general.h:747
virtual void PurificationStart()
Start recursive expansion.
Definition: purification_general.h:825
VectorTypeReal EIG_REL_GAP_LUMO_VEC
Definition: purification_general.h:442
real trunc_time
Definition: puri_info.h:59
int EIG_PROJECTION_INT
Definition: purification_general.cc:55
void save_matrix_now(string str)
Definition: purification_general.h:2070
generalVector VectorType
Definition: GetDensFromFock.cc:62
#define ORDER
Order of convergence of recursive expansion.
Definition: constants.h:41
Treal template_blas_fabs(Treal x)
bool compute_eigenvectors_in_this_SCF_cycle
Definition: purification_general.h:481
static real get_max_double()
Get largest number.
Definition: purification_general.h:202
ergo_real real
Definition: test.cc:46
int additional_iterations
Do a few more iterations after convergence.
Definition: purification_general.h:402
void compute_spectrum_bounds()
Compute spectrum bounds.
Definition: purification_general.h:1630
File containing declarations of functions for reading/writing sparse matrices from/to mtx (MatrixMark...
void set_spectrum_bounds(real eigmin, real eigmax)
Set spectrum bounds.
Definition: purification_general.h:1606
int iter_for_homo
Definition: purification_general.h:494
real frob_diff_time
Definition: puri_info.h:66
virtual void clear()
Clear all matrices in the class.
Definition: purification_general.h:177
VectorTypeReal SIGMA_HOMO_VEC
Definition: purification_general.h:440
void unset_compute_eigenvectors_in_each_iteration()
Definition: purification_general.h:240
#define UNIT_one_eV
Definition: units.h:45
#define max(a, b)
Definition: integrator.cc:87
int nocc
Number of occupied orbitals.
Definition: purification_general.h:407
VectorType * eigVecHOMO
Definition: purification_general.h:485
intervalType IntervalType
Definition: purification_general.h:115
real Xsquare_time
Definition: puri_info.h:58
void gen_matlab_file_eigs(const char *filename) const
Create MATLAB .m file which plots the homo and lumo bounds in each recursive expansion iteration...
Definition: purification_general.h:3374
void discard_homo_eigenvector()
Definition: purification_general.h:1353
Treal template_blas_log(Treal x)
Definition of the main floating-point datatype used; the ergo_real type.
#define LOG_CAT_INFO
Definition: output.h:49
Definition: matInclude.h:139
NormType normPuriStopCrit
Norm used in the stopping criterion Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm.
Definition: purification_general.h:415
virtual bool puri_is_prepared() const
Check is function prepare_to_purification() is already called.
Definition: purification_general.h:156
double get_nnz_X(size_t &nnzX)
Get nnz of X in %.
Definition: purification_general.h:306
Time-measuring class.
Definition: utilities.h:80
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
real NNZ_X
Definition: puri_info.h:80
virtual void prepare_to_purification()
Prepare data for recursive expansion.
Definition: purification_general.h:863
void gen_python_file_nnz(const char *filename) const
Create PYTHON .py file which plots number of non-zero elements in matrices X_i and X_i^2 in each recu...
Definition: purification_general.h:3582
int min(int a, int b)
Definition: lin_trans.cc:66
int EIG_EMPTY
Definition: purification_general.cc:53
void compute_eigenvector(MatrixType const &M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
Definition: purification_general.h:2815
int iter_for_lumo
Definition: purification_general.h:495
real orbital_homo_time
Definition: puri_info.h:67
VectorTypeReal SIGMA_LUMO_VEC
(Eigenvectors) Approximation of shifts in each iteration.
Definition: purification_general.h:440
static real get_min_double()
Get smallest number.
Definition: purification_general.h:206
void print(int area, const char *routine)
Definition: utilities.h:111
virtual real apply_poly(const int it, real x)=0
IntervalType lumo_bounds_F
Initial homo bounds for F.
Definition: purification_general.h:456
bool compute_eigenvectors_in_each_iteration
Compute homo and lumo eigenpairs in every iteration and save eigenvectors in txt files.
Definition: purification_general.h:505
real purify_time
Definition: puri_info.h:60
void gen_matlab_file_norm_diff(const char *filename) const
Create MATLAB .m file which plots the idempotency error in each recursive expansion iteration...
Definition: purification_general.h:3134
void compute_eigenvectors_without_diagonalization_last_iter_proj()
Definition: purification_general.h:2399
MatrixType Xsq
Matrix X^2.
Definition: purification_general.h:250
void find_eig_gaps_every_iter()
Definition: purification_general.h:2743
virtual real total_subspace_error(int it)
Definition: purification_general.h:1871
Definition: matInclude.h:139
void map_bounds_to_0_1()
Get eigenvalue bounds for X0.
Definition: purification_general.h:1519
void gen_matlab_file_nnz(const char *filename) const
Create MATLAB .m file which plots the number of non-zero elements in matrices X_i and X_i^2 in each r...
Definition: purification_general.h:3276
double get_nnz_X()
Get nnz of X in %.
Definition: purification_general.h:310
virtual void purification_process()
Run recursive expansion.
Definition: purification_general.h:913
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
int eigenvectors_method
Chosen method for computing eigenvectors.
Definition: purification_general.h:467
int EIG_POWER_INT
Definition: purification_general.cc:56
int it
Definition: puri_info.h:56
void writeToTmpFile(MatrixType &A) const
Definition: purification_general.h:2952
VectorTypeReal EIG_REL_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for homo eigenvalue.
Definition: purification_general.h:442
VectorType * eigVecLUMO
Definition: purification_general.h:484
void compute_eigenvectors_without_diagonalization(int it, IterationInfo &iter_info)
Compute HOMO and LUMO eigenvalues and eigenvectors of the matrix F.
Definition: purification_general.h:2204
virtual void purify_X(const int it)=0
Constants for conversion between units for some common units like Angstrom, electron-volt (eV)...
real eigValHOMO
Definition: purification_general.h:492
real constantC
Definition: puri_info.h:93
VectorType eigVecHOMORef
Definition: purification_general.h:488
Definition: matInclude.h:138
int EIG_SQUARE_INT
Definition: purification_general.cc:54
double get_elapsed_wall_seconds()
Definition: utilities.h:107
MatrixType X
Matrix X.
Definition: purification_general.h:248
int write_matrix_to_mtx(const char *filename, const vector< int > &I, const vector< int > &J, const vector< real > &val, const int &N)
Definition: files_sparse.cc:151
int maxit
Maximum number of iterations.
Definition: purification_general.h:404
real eucl_acc
Tolerance used for computation of spectral norm.
Definition: purification_general.cc:42
virtual void truncate_matrix(real &thresh, int it)=0
bool puri_is_prepared_flag
Definition: purification_general.h:399
Header file with typedefs for matrix and vector types.
real mixed_acc
Tolerance used for computation of mixed norm.
Definition: purification_general.cc:43
virtual void check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace, int &stop, real &estim_order)
Check stopping criterion.
Definition: purification_general.h:1810
Definition: puri_info.h:136
void check_eigenvectors_at_the_end()
Definition: purification_general.h:1382
real mixed_diff_time
Definition: puri_info.h:65
NormType normPuriTrunc
Norm used for the truncation of matrices.
Definition: purification_general.h:411
IntervalType lumo_bounds_F_new
Definition: purification_general.h:459
bool try_eigv_on_next_iteration_if_fail
Definition: purification_general.h:482
virtual void estimate_number_of_iterations(int &estim_num_iter)=0
IntervalType homo_bounds_F
Initial lumo bounds for F.
Definition: purification_general.h:455
double get_nnz_Xsq()
Get nnz of X^2 in %.
Definition: purification_general.h:318
real homo_eig_solver_time
Definition: puri_info.h:71
Definition: puri_info.h:52
int scf_step
Definition: purification_general.h:503
void compute_X()
Get matrix X0 by mapping spectrum of F into [0,1] in reverse order.
Definition: purification_general.h:1500
real lumo_eig_solver_time
Definition: puri_info.h:72
void estimate_homo_lumo(const VectorTypeReal &XmX2_norm_mixed, const VectorTypeReal &XmX2_norm_frob, const VectorTypeReal &XmX2_trace)
Get homo and lumo bounds from traces and norms of Xi-Xi^2.
Definition: purification_general.h:1992
real eigValLUMO
Definition: purification_general.h:491
VectorTypeInt really_good_iterations_lumo
Iterations where lumo eigenvector is actually computed.
Definition: purification_general.h:501
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:246
virtual void apply_poly_to_eigs(const int it, real &homo, real &lumo)=0
void get_frob_norm_est(const VectorTypeReal &XmX2_norm_frob, const std::vector< real > &h_in, const std::vector< real > &l_in, VectorTypeReal &YmY2_norm_frob_est)
Definition: purification_general.h:1964
void gen_matlab_file_cond_num(const char *filename) const
Create MATLAB .m file which plots a condition number of a problem of computing the density matrix in ...
Definition: purification_general.h:3335
int get_exact_number_of_puri_iterations()
Definition: purification_general.h:1899
real XmX2_eucl
Definition: puri_info.h:76
real XmX2_fro_norm
Definition: puri_info.h:74
void gen_matlab_file_time(const char *filename) const
Create MATLAB .m file which creates a bar plot presenting time spent on various parts of the iteratio...
Definition: purification_general.h:3443
Definition: matInclude.h:139
real trace_diff_time
Definition: puri_info.h:64
VectorTypeInt VecPoly
Polynomials computed in the function estimated_number_of_iterations() VecPoly[i] = 1 if we use X=X^2 ...
Definition: purification_general.h:431
double get_nnz_Xsq(size_t &nnzXsq)
Get nnz of X^2 in %.
Definition: purification_general.h:314
mat::normType NormType
Definition: purification_general.h:116
virtual real compute_derivative(const int it, real x, real &DDf)=0
real order
Definition: puri_info.h:77
VectorType eigVecLUMORef
Definition: purification_general.h:487
real NNZ_X2
Definition: puri_info.h:81
void get_spectrum_bounds(real &eigmin, real &eigmax)
Get spectrum bounds.
Definition: purification_general.h:1616
real eucl_diff_time
Definition: puri_info.h:63
real XmX2_trace
Definition: puri_info.h:73
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
virtual void return_constant_C(const int it, real &Cval)=0
int number_of_eigenvalues
Definition: purification_general.h:447
void get_eigenvalue_estimates(const VectorTypeReal &XmX2_norm_mixed, const VectorTypeReal &XmX2_norm_frob, const VectorTypeReal &XmX2_trace)
Definition: purification_general.h:1940
void set_truncation_parameters()
Definition: purification_general.h:1578
static std::string writeAndReadAll()
Definition: FileWritable.cc:509
void set_eigenvectors_params(string eigenvectors_method_, string eigenvectors_iterative_method_, real eigensolver_accuracy_, int eigensolver_maxiter_, int scf_step_, int use_prev_vector_as_initial_guess_, int try_eigv_on_next_iteration_if_fail_, VectorType *eigVecLUMO_, VectorType *eigVecHOMO_)
Set parameters for computing eigenvectors.
Definition: purification_general.h:649
void output_current_memory_usage(int logArea, const char *contextString)
Definition: output.cc:186
virtual void save_other_iter_info(IterationInfo &iter_info, int it)=0
real gap
Definition: puri_info.h:79
MatrixType F
Matrix F.
Definition: purification_general.h:255
int get_est_number_of_puri_iterations()
Definition: purification_general.h:1892
File containing classes IterationInfo and PuriInfo.
void discard_lumo_eigenvector()
Definition: purification_general.h:1369
real orbital_lumo_time
Definition: puri_info.h:68
#define NUM_ADDITIONAL_ITERATIONS
Definition: purification_general.h:76
void compute_eigenvectors_without_diagonalization_on_F(const MatrixType &F, int eigensolver_maxiter_for_F)
Definition: purification_general.h:2139
void find_shifts_every_iter()
/brief Find shifts sigma which will be used for construction of the filtering polynomial for computin...
Definition: purification_general.h:2714
string eigenvectors_method_str
Definition: purification_general.h:476
bool initialized_flag
Definition: purification_general.h:398
Treal template_blas_sqrt(Treal x)
IntervalType homo_bounds_F_new
Definition: purification_general.h:458
real error_eig
Error in eigenvalues (used just in old stopping criterion).
Definition: purification_general.h:420
int eigenvectors_iterative_method
Chosen eigensolver.
Definition: purification_general.h:468
void determine_iteration_for_eigenvectors()
Definition: purification_general.h:2599
normType
Definition: matInclude.h:139
virtual void eigenvalue_bounds_estimation()
Estimate eigenvalues near homo-lumo gap.
Definition: purification_general.h:1316
real stopping_criterion_time
Definition: puri_info.h:62
void save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it)
Definition: purification_general.h:2972
void readFromTmpFile(MatrixType &A) const
Definition: purification_general.h:2962
real threshold_X
Definition: puri_info.h:57
IntervalType homo_bounds_X0
Initial lumo bounds for X.
Definition: purification_general.h:452
void propagate_values_in_each_iter(real value_unocc, real value_occ, VectorTypeReal &out_unocc, VectorTypeReal &out_occ, int nmax)
Definition: purification_general.h:1914