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
30 /** @file purification_general.h
31
32 @brief Recursive density matrix expansion (or density matrix
33 purification).
34
35 @author Anastasia Kruchinina <em>responsible</em>
36 */
37
38
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"
49 #include "integral_matrix_wrappers.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;
102 extern real TOL_OVERLAPPING_BOUNDS;
103 extern real THRESHOLD_EIG_TOLERANCE;
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
111 /** PurificationGeneral is an abstract class which provides an
112 * interface for SP2, SP2ACC and possibly other recursive expansions.
113 *
114 * \tparam MatrixType Type of a matrix (ex. symmMatrix). */
115 template<typename MatrixType>
116 class PurificationGeneral
117 {
118 public:
119 typedef ergo_real real;
120 typedef intervalType IntervalType;
121 typedef mat::normType NormType;
122 typedef std::vector<int> VectorTypeInt;
123 typedef std::vector<real> VectorTypeReal;
124 typedef generalVector VectorType;
125 typedef MatrixType MatrixTypeWrapper;
126
127
128 PuriInfo info; /**< Fill in during purification with useful information. */
129
130 MatrixType X; /**< Matrix X. */
131
132 MatrixType Xsq; /**< Matrix X^2. */
133
134
PurificationGeneral()135 PurificationGeneral()
136 {
137 initialized_flag = false;
138 puri_is_prepared_flag = false;
139 computed_spectrum_bounds = false;
140
141 // eigenvectors staff
142 number_of_occ_eigenvectors = -1;
143 number_of_unocc_eigenvectors = -1;
144 compute_eigenvectors_in_each_iteration = false;
145 jump_over_X_iter_proj_method = -1;
146 go_back_X_iter_proj_method = -1;
147
148 use_new_stopping_criterion = false;
149 }
150
151 /** Set imporatant parameters for the recursive expansion.
152 */
153 virtual void initialize(const MatrixType& F_, /**< [in] Effective Hamiltonian matrix. */
154 const IntervalType& lumo_bounds_, /**< [in] Bounds for lumo eigenvalue. */
155 const IntervalType& homo_bounds_, /**< [in] Bounds for homo eigenvalue. */
156 int maxit_, /**< [in] Maximum number of recursive expansion iterations. */
157 real error_sub_, /**< [in] Allowed error in invariant subspaces. */
158 real error_eig_, /**< [in] Error in eigenvalues (used just in old stopping criterion). */
159 bool use_new_stopping_criterion_, /**< [in] Set if want to use new stopping criterion. */
160 NormType norm_truncation, /**< [in] Truncation norm. */
161 NormType norm_stop_crit, /**< [in] Stopping criterion norm. */
162 int nocc_ /**< [in] Number of occupied orbitals. */
163 );
164
165
166
167 /** Start recursive expansion.
168 */
169 virtual void PurificationStart();
170
171
172 /** Clear all matrices in the class.
173 * Needed to be called if Chunks and Tasks are used,
174 * since we need to delete all ChunkIDs
175 * before exiting the program. */
clear()176 virtual void clear() { X.clear(); Xsq.clear(); }
177
178
179 /** Set spectrum bounds.
180 *
181 * Used if we know spectrum bounds or want to compute them outside
182 * the class. */
183 void set_spectrum_bounds(real eigmin, real eigmax);
184
185 /** Get spectrum bounds.
186 *
187 * Return spectrum bounds for the matrix F. If not computed before, compute them. */
188 void get_spectrum_bounds(real& eigmin, real& eigmax);
189
190 int get_exact_number_of_puri_iterations();
191 int get_est_number_of_puri_iterations();
192
193 virtual real total_subspace_error(int it);
194
195 /** Get machine epsilon. */
get_epsilon()196 static real get_epsilon()
197 { return mat::getMachineEpsilon<real>(); }
198
199 /** Get largest number. */
get_max_double()200 static real get_max_double()
201 { return std::numeric_limits<real>::max(); }
202
203 /** Get smallest number. */
get_min_double()204 static real get_min_double()
205 { return std::numeric_limits<real>::min(); }
206
207
208
209 /** Set parameters for computing eigenvectors. */
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
218 /** Set parameters for computing eigenvectors - long version. */
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
extract_computed_eigenpairs(std::vector<VectorType> & eigVecUNOCCref,std::vector<VectorType> & eigVecOCCref,std::vector<real> & eigValUNOCCref,std::vector<real> & eigValOCCref)232 void extract_computed_eigenpairs(
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
set_compute_eigenvectors_in_each_iteration()245 void set_compute_eigenvectors_in_each_iteration() { compute_eigenvectors_in_each_iteration = true; }
unset_compute_eigenvectors_in_each_iteration()246 void unset_compute_eigenvectors_in_each_iteration() { compute_eigenvectors_in_each_iteration = false; }
247
set_number_of_eigenvectors_to_compute(const int occ,const int unocc)248 void set_number_of_eigenvectors_to_compute(const int occ, const int unocc) {
249 number_of_occ_eigenvectors = std::max(occ, 0);
250 number_of_unocc_eigenvectors = std::max(unocc, 0);
251 }
get_number_of_occ_eigenvectors_to_compute()252 int get_number_of_occ_eigenvectors_to_compute() const {
253 return number_of_occ_eigenvectors; }
get_number_of_unocc_eigenvectors_to_compute()254 int get_number_of_unocc_eigenvectors_to_compute() const {
255 return number_of_unocc_eigenvectors; }
256
set_jump_over_X_iter_proj_method(int val)257 void set_jump_over_X_iter_proj_method(int val)
258 { jump_over_X_iter_proj_method = val;}
get_jump_over_X_iter_proj_method()259 int get_jump_over_X_iter_proj_method() const
260 { return jump_over_X_iter_proj_method;}
set_go_back_X_iter_proj_method(int val)261 void set_go_back_X_iter_proj_method(int val)
262 { go_back_X_iter_proj_method = val;}
get_go_back_X_iter_proj_method()263 int get_go_back_X_iter_proj_method() const
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
270 /** Create MATLAB .m file which plots the idempotency error in each recursive expansion iteration. */
271 void gen_matlab_file_norm_diff(const char *filename) const;
272 /** Create MATLAB .m file which plots the actual introduced error after truncation of the matrix X_i in each recursive expansion iteration. */
273 void gen_matlab_file_threshold(const char *filename) const;
274 /** Create MATLAB .m file which plots the number of non-zero elements in matrices X_i and X_i^2 in each recursive expansion iteration. */
275 void gen_matlab_file_nnz(const char *filename) const;
276 /** Create MATLAB .m file which plots the homo and lumo bounds in each recursive expansion iteration. */
277 void gen_matlab_file_eigs(const char *filename) const;
278 /** Create MATLAB .m file which creates a bar plot presenting time spent on various parts of the iteration (such as matrix square and computation of eigenvectors) in each recursive expansion iteration. */
279 void gen_matlab_file_time(const char *filename) const;
280 /** Create MATLAB .m file which plots a condition number of a problem of computing the density matrix in each recursive expansion iteration. The condition number is equal to inverse of the homo-lumo gap approximation. */
281 // void gen_matlab_file_cond_num(const char *filename) const;
282 /** Create PYTHON .py file which plots number of non-zero elements in matrices X_i and X_i^2 in each recursive expansion iteration. */
283 // void gen_python_file_nnz(const char *filename) const;
284
285
286
~PurificationGeneral()287 virtual ~PurificationGeneral() {}
288
289 protected:
290
291 MatrixType F; /**< Matrix F. Needed for computation of eigenvectors.*/
292
293 std::vector<MatrixType> vec_matrices_Xi; /**< Save matrices Xi in each iteration (if used projection method for computing eigenvectors). If ergo input file mat.write_to_file=1 than matrices are saved to the file to save memory. */
294
295
296 /** Check is function initialize() is already called.
297 */
is_initialized()298 virtual bool is_initialized() const { return initialized_flag; }
299
300 /** Check is function prepare_to_purification() is already called.
301 */
puri_is_prepared()302 virtual bool puri_is_prepared() const { return puri_is_prepared_flag; }
303
304 /** Prepare data for recursive expansion. */
305 virtual void prepare_to_purification();
306
307 /** Prepare data related to the eigenvectors computations. */
308 virtual void prepare_to_purification_eigenvectors();
309
310 /** Run recursive expansion. */
311 virtual void purification_process();
312
313 /** Estimate eigenvalues near homo-lumo gap. */
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
320 /** Compute spectrum bounds. */
321 virtual void compute_spectrum_bounds();
322
323 /** Get matrix X0 by mapping spectrum of F into [0,1] in reverse
324 * order.*/
325 virtual void compute_X();
326
327 /** Get eigenvalue bounds for X0. */
328 void map_bounds_to_0_1();
329
330 /** Check stopping criterion (obsolete).
331 *
332 * Use stopping criterion based on user defined threshold
333 * values. */
334 virtual void check_standard_stopping_criterion(const real XmX2_norm, int& stop);
335
336 /** Check stopping criterion.
337 *
338 * The new stopping criterion based on the order of convergence is
339 * used, see article "Parameterless Stopping Criteria for Recursive Density
340 * Matrix Expansions", J. Chem. Theory Comput., 2016, 12 (12), pp 5788–5802
341 * DOI: 10.1021/acs.jctc.6b00626 */
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
345 /** Choose stopping criterion and check it. */
346 virtual void stopping_criterion(IterationInfo& iter_info, int& stop, real& estim_order);
347
348 int get_int_eig_iter_method(string eigenvectors_iterative_method);
349 int get_int_eig_method(string eigenvectors_method);
350
351 /** Compute HOMO and LUMO eigenvalues and eigenvectors of the matrix F.
352 *
353 * The method uses the polynomial constructed during the recursive expansion,
354 * so no additional matrix multiplications are required. See article
355 * J. Chem. Theory Comput., Just Accepted Manuscript,
356 * DOI: 10.1021/acs.jctc.7b00968
357 */
358 void compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info);
359 void compute_eigenvectors_without_diagonalization_last_iter_proj();
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
365 /** Set parameters for computing eigenvectors. */
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
376 /** Get nnz of X in %. */
get_nnz_X(size_t & nnzX)377 double get_nnz_X(size_t& nnzX /**< Number of nz of X*/)
378 { nnzX = X.nnz(); return (double)(((double)nnzX) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
379
380 /** Get nnz of X in %. */
get_nnz_X()381 double get_nnz_X()
382 { return (double)(((double)X.nnz()) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
383
384 /** Get nnz of X^2 in %. */
get_nnz_Xsq(size_t & nnzXsq)385 double get_nnz_Xsq(size_t& nnzXsq /**< Number of nz of X^2*/)
386 { nnzXsq = Xsq.nnz(); return (double)(((double)nnzXsq) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
387
388 /** Get nnz of X^2 in %. */
get_nnz_Xsq()389 double get_nnz_Xsq()
390 { return (double)(((double)Xsq.nnz()) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
391
392 /** Get homo and lumo bounds from traces and norms of Xi-Xi^2. Called by get_eigenvalue_estimates().
393 *
394 * Used at the end of the recursive expansion. See article SIAM
395 * J. Sci. Comput., 36(2), B147–B170. */
396 void estimate_homo_lumo(const VectorTypeReal& XmX2_norm_mixed,
397 const VectorTypeReal& XmX2_norm_frob,
398 const VectorTypeReal& XmX2_trace);
399
400 /** Get homo and lumo bounds from traces and norms of Xi-Xi^2.
401 *
402 * Used at the end of the recursive expansion. See article SIAM
403 * J. Sci. Comput., 36(2), B147–B170. */
404 void get_eigenvalue_estimates(const VectorTypeReal& XmX2_norm_mixed,
405 const VectorTypeReal& XmX2_norm_frob,
406 const VectorTypeReal& XmX2_trace);
407
408
409 /** Determine in which iterations will be computed homo and lumo eigenvectors. Calling get_iterations_for_lumo_and_homo(). */
410 virtual void determine_iteration_for_eigenvectors();
411
412 /** \brief Find the best iterations for computing eigenvectors.
413 *
414 * The best iteration defined by both a small error in the eigenvector
415 * and a small number of iterations of an iterative method.
416 */
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;
425 void output_separate_total_times(PuriInfo &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
437 /** /brief Find shifts sigma which will be used for construction of
438 * the filtering polynomial for computing eigenvectors.
439 */
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
475 bool initialized_flag;
476 bool puri_is_prepared_flag;
477
478 bool use_new_stopping_criterion; /**< True for new stopping criterion */
479 int additional_iterations; /**< Do a few more iterations after convergence */
480
481 int maxit; /**< Maximum number of iterations */
482 int check_stopping_criterion_iter; /**< Iteration when to start to check stopping criterion. */
483
484 int nocc; /**< Number of occupied orbitals */
485
486
487
488 NormType normPuriTrunc; /**< Norm used for the truncation of matrices.
489 * Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm. */
490
491
492 NormType normPuriStopCrit; /**< Norm used in the stopping criterion
493 * Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm. */
494
495
496 real error_sub; /**< Allowed error in invariant subspaces. */
497 real error_eig; /**< Error in eigenvalues (used just in old stopping criterion). */
498
499 real error_per_it; /**< Error allowed in each iteration due to truncation. */
500
501 real constant_C; /**< Asymptotic constant C needed for the new stopping criterion. */
502
503 real gammaStopEstim; /**< Used on the stopping criterion for
504 * estimation of eigenvalues from
505 * purification */
506
507
508 VectorTypeInt VecPoly; /**< Polynomials computed in the function
509 * estimated_number_of_iterations()
510 * VecPoly[i] = 1 if we use X=X^2
511 * VecPoly[i] = 0 if we use X=2X-X^2
512 * (or their accelerated versions) */
513 VectorTypeReal VecGap; /**< Gap computed using inner homo and lumo bounds on each iteration. */
514
515
516 VectorTypeReal ITER_ERROR_VEC; /**< (Eigenvectors) Maximum error introduced in each iteration. */
517 VectorTypeReal SIGMA_HOMO_VEC, SIGMA_LUMO_VEC; /**< (Eigenvectors) Approximation of shifts in each iteration. */
518 VectorTypeReal EIG_ABS_GAP_LUMO_VEC, EIG_ABS_GAP_HOMO_VEC; /**< (Eigenvectors) Absolute and relative gap in filter for lumo eigenvalue. */
519 VectorTypeReal EIG_REL_GAP_LUMO_VEC, EIG_REL_GAP_HOMO_VEC; /**< (Eigenvectors) Absolute and relative gap in filter for homo eigenvalue. */
520
521
522 /*EIGENVECTORS STAFF*/
523
524 int number_of_occ_eigenvectors;
525 int number_of_unocc_eigenvectors;
526
527 int jump_over_X_iter_proj_method;
528 int go_back_X_iter_proj_method;
529
530
531 IntervalType homo_bounds; /**< (1-homo) bounds for Xi in iteration i */
532 IntervalType lumo_bounds; /**< Lumo bounds for Xi in iteration i */
533
534 IntervalType homo_bounds_X0; /**< Initial lumo bounds for X */
535 IntervalType lumo_bounds_X0; /**< Initial lumo bounds for X */
536
537 IntervalType homo_bounds_F; /**< Initial lumo bounds for F */
538 IntervalType lumo_bounds_F; /**< Initial homo bounds for F */
539
540 IntervalType homo_bounds_F_new;
541 IntervalType lumo_bounds_F_new;
542
543
544 IntervalType spectrum_bounds; /**< Outer bounds for the whole spectrum of F/Xi. */
545 bool computed_spectrum_bounds;
546
547
548 /* Eigenvectors */
549
550 int eigenvectors_method; /**< Chosen method for computing eigenvectors. */
551 int eigenvectors_iterative_method; /**< Chosen eigensolver. */
552
553 real eigensolver_accuracy; /**< Accuracy of the eigenvalue problem solver. When residual is less then this value, stop. */
554
555 int eigensolver_maxiter; /**< Maximum number of iterations for eigensolver. */
556
557 string eigenvectors_method_str;
558 string eigenvectors_iterative_method_str;
559
560 bool use_prev_vector_as_initial_guess;
561
562 bool compute_eigenvectors_in_this_SCF_cycle;
563 bool try_eigv_on_next_iteration_if_fail;
564
565 VectorType eigVecLUMORef;
566 VectorType eigVecHOMORef;
567
568 std::vector<VectorType> eigVecOCC; /**< Here we save eigenvectors corresponding to the occupied orbitals. */
569 std::vector<VectorType> eigVecUNOCC; /**< Here we save eigenvectors corresponding to the unoccupied orbitals. */
570 std::vector<real> eigValOCC; /**< Here we save eigenvalues corresponding to the occupied orbitals. */
571 std::vector<real> eigValUNOCC; /**< Here we save eigenvalues corresponding to the unoccupied orbitals. */
572
573
574 real eigValLUMO;
575 real eigValHOMO;
576
577 int iter_for_homo;
578 int iter_for_lumo;
579
580 VectorTypeInt good_iterations_homo; /**< Iterations where homo eigenvector can be computed. */
581 VectorTypeInt good_iterations_lumo; /**< Iterations where homo eigenvector can be computed. */
582
583 VectorTypeInt really_good_iterations_homo; /**< Iterations where homo eigenvector is actually computed. */
584 VectorTypeInt really_good_iterations_lumo; /**< Iterations where lumo eigenvector is actually computed. */
585
586 int scf_step;
587
588 bool compute_eigenvectors_in_each_iteration; /**< Compute homo and lumo eigenpairs in every iteration
589 * and save eigenvectors in txt files.
590 *
591 * NOTE: if we are computing eigenvector in every iteration, we will
592 * save eigenvector computed in the last iteration, and it is most
593 * probably the wrong eigenvector.
594 */
595 };
596
597 /**************************************/
598
599
600 /******************* INIT *******************/
601
602
603 template<typename MatrixType>
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_)604 void PurificationGeneral<MatrixType>::initialize(const MatrixType& F_,
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
640 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:
654 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "mixed");
655 break;
656
657 case mat::euclNorm:
658 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "eucl");
659 break;
660
661 case mat::frobNorm:
662 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "frob");
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:
682 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "mixed");
683 break;
684
685 case mat::euclNorm:
686 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "eucl");
687 break;
688
689 case mat::frobNorm:
690 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "frob");
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>
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_)737 void PurificationGeneral<MatrixType>::set_eigenvectors_params_basic(string eigenvectors_method_,
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>
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_)802 void PurificationGeneral<MatrixType>::set_eigenvectors_params(string eigenvectors_method_,
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>
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_,bool use_prev_vector_as_initial_guess_,const std::vector<VectorType> & eigVecOCCRef,const std::vector<VectorType> & eigVecUNOCCRef)815 void PurificationGeneral<MatrixType>::set_eigenvectors_params(string eigenvectors_method_,
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);
840 mat::SizesAndBlocks cols;
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);
853 mat::SizesAndBlocks cols;
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>
get_int_eig_method(string eigenvectors_method)867 int PurificationGeneral<MatrixType>::get_int_eig_method(string eigenvectors_method)
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>
get_int_eig_iter_method(string eigenvectors_iterative_method)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>
output_norms_and_traces(IterationInfo & iter_info)905 void PurificationGeneral<MatrixType>::output_norms_and_traces(IterationInfo& iter_info) const
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>
PurificationStart()953 void PurificationGeneral<MatrixType>::PurificationStart()
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>
prepare_to_purification_eigenvectors()1011 void PurificationGeneral<MatrixType>::prepare_to_purification_eigenvectors()
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>
prepare_to_purification()1065 void PurificationGeneral<MatrixType>::prepare_to_purification()
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>
purification_process()1100 void PurificationGeneral<MatrixType>::purification_process()
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);
1164 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
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);
1174 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
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);
1343 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
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);
1353 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,
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>
output_time_WriteAndReadAll()1571 void PurificationGeneral<MatrixType>::output_time_WriteAndReadAll() const
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>
output_separate_total_times(PuriInfo & info)1582 void PurificationGeneral<MatrixType>::output_separate_total_times(PuriInfo &info) const
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>
eigenvalue_bounds_estimation()1625 void PurificationGeneral<MatrixType>::eigenvalue_bounds_estimation()
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>
discard_homo_eigenvector()1668 void PurificationGeneral<MatrixType>::discard_homo_eigenvector()
1669 {
1670 eigVecOCC.clear();
1671 info.eigValHOMO = 0;
1672 info.homo_eigenvector_is_computed = false;
1673 }
1674
1675
1676 template<typename MatrixType>
discard_lumo_eigenvector()1677 void PurificationGeneral<MatrixType>::discard_lumo_eigenvector()
1678 {
1679 eigVecUNOCC.clear();
1680 info.eigValLUMO = 0;
1681 info.lumo_eigenvector_is_computed = false;
1682 }
1683
1684
1685 template<typename MatrixType>
check_eigenvectors_at_the_end()1686 void PurificationGeneral<MatrixType>::check_eigenvectors_at_the_end()
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>
compute_X()1804 void PurificationGeneral<MatrixType>::compute_X()
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>
map_bounds_to_0_1()1823 void PurificationGeneral<MatrixType>::map_bounds_to_0_1()
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>
truncate_matrix(real & threshold,int it)1882 void PurificationGeneral<MatrixType>::truncate_matrix(real& threshold, int it)
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>
set_truncation_parameters()1926 void PurificationGeneral<MatrixType>::set_truncation_parameters()
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>
set_spectrum_bounds(real eigmin,real eigmax)1954 void PurificationGeneral<MatrixType>::set_spectrum_bounds(real eigmin, real eigmax)
1955 {
1956 spectrum_bounds = IntervalType(eigmin, eigmax);
1957 computed_spectrum_bounds = true;
1958 }
1959
1960
1961 /****************** GET_SPECTRUM_BOUNDS *****************/
1962
1963 template<typename MatrixType>
get_spectrum_bounds(real & eigmin,real & eigmax)1964 void PurificationGeneral<MatrixType>::get_spectrum_bounds(real& eigmin, real& eigmax)
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>
compute_spectrum_bounds()1978 void PurificationGeneral<MatrixType>::compute_spectrum_bounds()
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>
stopping_criterion(IterationInfo & iter_info,int & stop,real & estim_order)2067 void PurificationGeneral<MatrixType>::stopping_criterion(IterationInfo& iter_info, int& stop, real& estim_order)
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>
check_standard_stopping_criterion(const real XmX2_norm,int & stop)2130 void PurificationGeneral<MatrixType>::check_standard_stopping_criterion(const real XmX2_norm, int& stop)
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>
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)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>
2212 typename PurificationGeneral<MatrixType>::real
total_subspace_error(int it)2213 PurificationGeneral<MatrixType>::total_subspace_error(int it)
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>
get_est_number_of_puri_iterations()2234 int PurificationGeneral<MatrixType>::get_est_number_of_puri_iterations()
2235 {
2236 return info.estim_total_it;
2237 }
2238
2239
2240 template<typename MatrixType>
get_exact_number_of_puri_iterations()2241 int PurificationGeneral<MatrixType>::get_exact_number_of_puri_iterations()
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>
get_eigenvalue_estimates(const VectorTypeReal & XmX2_norm_mixed,const VectorTypeReal & XmX2_norm_frob,const VectorTypeReal & XmX2_trace)2257 void PurificationGeneral<MatrixType>::get_eigenvalue_estimates(const VectorTypeReal& XmX2_norm_mixed,
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>
estimate_homo_lumo(const VectorTypeReal & XmX2_norm_mixed,const VectorTypeReal & XmX2_norm_frob,const VectorTypeReal & XmX2_trace)2270 void PurificationGeneral<MatrixType>::estimate_homo_lumo(const VectorTypeReal& XmX2_norm_mixed,
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>
save_matrix_now(string str)2350 void PurificationGeneral<MatrixType>::save_matrix_now(string str)
2351 {
2352 save_matrix_A_now(X, str);
2353 }
2354
2355
2356 template<typename MatrixType>
save_matrix_A_now(const MatrixType & A,string str)2357 void PurificationGeneral<MatrixType>::save_matrix_A_now(const MatrixType &A, string str)
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>
compute_eigenvectors_without_diagonalization_on_F(const MatrixType & F,int eigensolver_maxiter_for_F)2403 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F)
2404 {
2405 throw std::runtime_error("compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
2406 }
2407
2408
2409 template<typename MatrixType>
compute_eigenvectors_without_diagonalization_last_iter_proj()2410 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_last_iter_proj()
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>
compute_eigenvectors_without_diagonalization(int it,IterationInfo & iter_info)2417 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info)
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>
compute_eigenvectors_without_diagonalization_on_F(const MatrixType & F,int eigensolver_maxiter_for_F)2426 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F)
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
2452 mat::SizesAndBlocks rows;
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
2483 mat::SizesAndBlocks rows;
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>
compute_eigenvectors_without_diagonalization(int it,IterationInfo & iter_info)2506 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info)
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. */
2642 mat::SizesAndBlocks rows;
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>
compute_eigenvectors_without_diagonalization_last_iter_proj()2677 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_last_iter_proj()
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>
projection_method_one_puri_iter(int current_iteration)2757 void PurificationGeneral<MatrixType>::projection_method_one_puri_iter(int current_iteration)
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>
determine_iteration_for_eigenvectors()2881 void PurificationGeneral<MatrixType>::determine_iteration_for_eigenvectors()
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>
get_iterations_for_lumo_and_homo(int & chosen_iter_lumo,int & chosen_iter_homo)2899 void PurificationGeneral<MatrixType>::get_iterations_for_lumo_and_homo(int& chosen_iter_lumo,
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>
find_shifts_every_iter()2962 void PurificationGeneral<MatrixType>::find_shifts_every_iter()
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>
2992 void PurificationGeneral<MatrixType>::find_truncation_thresh_every_iter()
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>
3013 void PurificationGeneral<MatrixType>::find_eig_gaps_every_iter()
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>
compute_eigenvector(MatrixType const & M,std::vector<VectorType> & eigVec,std::vector<real> & eigVal,int it,bool is_homo)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>
compute_eigenvector(MatrixType const & M,std::vector<VectorType> & eigVec,std::vector<real> & eigVal,int it,bool is_homo)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. */
3122 mat::SizesAndBlocks rows;
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>
3293 void PurificationGeneral<MatrixType>::
writeToTmpFile(MatrixType & A)3294 writeToTmpFile(MatrixType& A) const
3295 {
3296 #ifndef USE_CHUNKS_AND_TASKS
3297 A.writeToFile();
3298 #endif
3299 }
3300
3301
3302 template<typename MatrixType>
3303 void PurificationGeneral<MatrixType>::
readFromTmpFile(MatrixType & A)3304 readFromTmpFile(MatrixType& A) const
3305 {
3306 #ifndef USE_CHUNKS_AND_TASKS
3307 A.readFromFile();
3308 #endif
3309 }
3310
3311
3312 template<typename MatrixType>
3313 void PurificationGeneral<MatrixType>::
save_selected_eigenvector_to_file(const VectorType & v,int num,bool is_homo,int it)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>
get_eigenvalue_of_F_from_eigv_of_Xi(real & eigVal,const VectorType & eigVec)3346 void PurificationGeneral<MatrixType>::get_eigenvalue_of_F_from_eigv_of_Xi(real& eigVal, const VectorType& eigVec)
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>
check_homo_lumo_eigenvalues(real & eigVal,VectorType & eigVec,bool & is_homo,bool & is_lumo,const int iter)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>
gen_matlab_file_norm_diff(const char * filename)3470 void PurificationGeneral<MatrixType>::gen_matlab_file_norm_diff(const char *filename) const
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>
gen_matlab_file_threshold(const char * filename)3566 void PurificationGeneral<MatrixType>::gen_matlab_file_threshold(const char *filename) const
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>
gen_matlab_file_nnz(const char * filename)3612 void PurificationGeneral<MatrixType>::gen_matlab_file_nnz(const char *filename) const
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>
gen_matlab_file_eigs(const char * filename)3710 void PurificationGeneral<MatrixType>::gen_matlab_file_eigs(const char *filename) const
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>
gen_matlab_file_time(const char * filename)3779 void PurificationGeneral<MatrixType>::gen_matlab_file_time(const char *filename) const
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
3982