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