00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00039 #ifndef HEADER_PURIFICATION_GENERAL
00040 #define HEADER_PURIFICATION_GENERAL
00041
00042 #include <iostream>
00043 #include <fstream>
00044 #include <sstream>
00045
00046 #include "matrix_typedefs.h"
00047 #include "realtype.h"
00048 #include "matrix_utilities.h"
00049 #include "integral_matrix_wrappers.h"
00050 #include "output.h"
00051 #include "matrix_proxy.h"
00052
00053 #include "puri_info.h"
00054 #include "constants.h"
00055 #include "utilities.h"
00056 #include "units.h"
00057
00058 #include "files_dense.h"
00059 #include "files_sparse.h"
00060 #include "files_sparse_bin.h"
00061
00062 #include "get_eigenvectors.h"
00063
00064 typedef ergo_real real;
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076 #define NUM_ADDITIONAL_ITERATIONS 0
00077
00078
00079
00080 #define PURI_OUTPUT_NNZ
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 extern real eucl_acc;
00096 extern real mixed_acc;
00097 extern real TOL_OVERLAPPING_BOUNDS;
00098 extern real THRESHOLD_EIG_TOLERANCE;
00099 extern int EIG_EMPTY;
00100 extern int EIG_SQUARE_INT;
00101 extern int EIG_PROJECTION_INT;
00102 extern int EIG_POWER_INT;
00103 extern int EIG_LANCZOS_INT;
00104
00105
00110 template<typename MatrixType>
00111 class PurificationGeneral
00112 {
00113 public:
00114 typedef ergo_real real;
00115 typedef intervalType IntervalType;
00116 typedef mat::normType NormType;
00117 typedef std::vector<int> VectorTypeInt;
00118 typedef std::vector<real> VectorTypeReal;
00119 typedef generalVector VectorType;
00120 typedef MatrixType MatrixTypeWrapper;
00121
00122 PurificationGeneral()
00123 {
00124 initialized_flag = false;
00125 puri_is_prepared_flag = false;
00126 computed_spectrum_bounds = false;
00127
00128
00129 number_of_eigenvalues = 0;
00130 compute_eigenvectors_in_each_iteration = false;
00131 eigVecHOMO = NULL;
00132 eigVecLUMO = NULL;
00133 }
00134
00137 virtual void initialize(const MatrixType& F_,
00138 const IntervalType& lumo_bounds_,
00139 const IntervalType& homo_bounds_,
00140 int maxit_,
00141 real error_sub_,
00142 real error_eig_,
00143 int use_new_stopping_criterion_,
00144 NormType norm_truncation,
00145 NormType norm_stop_crit,
00146 int nocc_
00147 );
00148
00149
00152 virtual bool is_initialized() const { return initialized_flag; }
00153
00156 virtual bool puri_is_prepared() const { return puri_is_prepared_flag; }
00157
00160 virtual void PurificationStart();
00161
00163 virtual void prepare_to_purification();
00164
00166 virtual void purification_process();
00167
00169 virtual void eigenvalue_bounds_estimation();
00170
00171 virtual ~PurificationGeneral() {}
00172
00177 virtual void clear() { X.clear(); Xsq.clear(); }
00178
00179
00184 void set_spectrum_bounds(real eigmin, real eigmax);
00185
00189 void get_spectrum_bounds(real& eigmin, real& eigmax);
00190
00191 int get_exact_number_of_puri_iterations();
00192 int get_est_number_of_puri_iterations();
00193
00194
00195 virtual real total_subspace_error(int it);
00196
00198 static real get_epsilon()
00199 { return mat::getMachineEpsilon<real>(); }
00200
00202 static real get_max_double()
00203 { return std::numeric_limits<real>::max(); }
00204
00206 static real get_min_double()
00207 { return std::numeric_limits<real>::min(); }
00208
00210 void gen_matlab_file_norm_diff(const char *filename) const;
00212 void gen_matlab_file_threshold(const char *filename) const;
00214 void gen_matlab_file_nnz(const char *filename) const;
00216 void gen_matlab_file_eigs(const char *filename) const;
00218 void gen_matlab_file_time(const char *filename) const;
00220 void gen_matlab_file_cond_num(const char *filename) const;
00221
00223 void gen_python_file_nnz(const char *filename) const;
00224
00225
00226
00228 void set_eigenvectors_params(string eigenvectors_method_,
00229 string eigenvectors_iterative_method_,
00230 real eigensolver_accuracy_,
00231 int eigensolver_maxiter_,
00232 int scf_step_,
00233 int use_prev_vector_as_initial_guess_,
00234 int try_eigv_on_next_iteration_if_fail_,
00235 VectorType *eigVecLUMO_,
00236 VectorType *eigVecHOMO_
00237 );
00238
00239 void set_compute_eigenvectors_in_each_iteration() { compute_eigenvectors_in_each_iteration = true; }
00240 void unset_compute_eigenvectors_in_each_iteration() { compute_eigenvectors_in_each_iteration = false; }
00241
00242
00243 void compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F);
00244
00245
00246 PuriInfo info;
00248 MatrixType X;
00250 MatrixType Xsq;
00253 protected:
00254
00255 MatrixType F;
00257 MatrixType X_homo, X_lumo;
00259 void save_matrix_now(string str);
00260
00262 void compute_spectrum_bounds();
00263
00266 void compute_X();
00267
00269 void map_bounds_to_0_1();
00270
00275 virtual void check_standard_stopping_criterion(const real XmX2_norm, int& stop);
00276
00283 virtual void check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace,
00284 int& stop, real& estim_order);
00285
00287 virtual void stopping_criterion(IterationInfo& iter_info, int& stop, real& estim_order);
00288
00289 int get_int_eig_iter_method(string eigenvectors_iterative_method);
00290 int get_int_eig_method(string eigenvectors_method);
00291
00299 void compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info);
00300 void compute_eigenvectors_without_diagonalization_last_iter_proj();
00301
00302 void compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo);
00303
00304
00306 double get_nnz_X(size_t& nnzX )
00307 { nnzX = X.nnz(); return (double)(((double)nnzX) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
00308
00310 double get_nnz_X()
00311 { return (double)(((double)X.nnz()) / ((double)X.get_ncols() * X.get_nrows()) * 100); }
00312
00314 double get_nnz_Xsq(size_t& nnzXsq )
00315 { nnzXsq = Xsq.nnz(); return (double)(((double)nnzXsq) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
00316
00318 double get_nnz_Xsq()
00319 { return (double)(((double)Xsq.nnz()) / ((double)Xsq.get_ncols() * Xsq.get_nrows()) * 100); }
00320
00325 void estimate_homo_lumo(const VectorTypeReal& XmX2_norm_mixed,
00326 const VectorTypeReal& XmX2_norm_frob,
00327 const VectorTypeReal& XmX2_trace);
00328
00329
00330 void get_frob_norm_est(const VectorTypeReal& XmX2_norm_frob,
00331 const std::vector<real>& h_in,
00332 const std::vector<real>& l_in,
00333 VectorTypeReal& YmY2_norm_frob_est);
00334
00335
00336
00337 void get_eigenvalue_estimates(const VectorTypeReal& XmX2_norm_mixed,
00338 const VectorTypeReal& XmX2_norm_frob,
00339 const VectorTypeReal& XmX2_trace);
00340
00341
00342 void propagate_values_in_each_iter(real value_unocc, real value_occ,
00343 VectorTypeReal& out_unocc,
00344 VectorTypeReal& out_occ,
00345 int nmax);
00346
00347
00348 void determine_iteration_for_eigenvectors();
00349
00350 void get_iterations_for_lumo_and_homo(int& chosen_iter_lumo,
00351 int& chosen_iter_homo);
00352
00353 void check_eigenvectors_at_the_end();
00354 void discard_homo_eigenvector();
00355 void discard_lumo_eigenvector();
00356
00357 void output_norms_and_traces(IterationInfo& iter_info) const;
00358
00359 void get_interval_with_lambda(real& eigVal, VectorType& eigVec, bool& is_homo, bool& is_lumo, const int iter);
00360 void get_eigenvalue_of_F_from_eigv_of_Xi(real& eigVal, const VectorType& eigVec);
00361
00362 void save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it);
00363
00364 void set_truncation_parameters();
00365
00366 void find_truncation_thresh_every_iter();
00367 void find_shifts_every_iter();
00368 void find_eig_gaps_every_iter();
00369
00370
00371 void writeToTmpFile(MatrixType& A) const;
00372 void readFromTmpFile(MatrixType& A) const;
00373
00374
00375
00376
00377
00378
00379
00380 virtual void set_init_params() = 0;
00381 virtual void truncate_matrix(real& thresh, int it) = 0;
00382 virtual void estimate_number_of_iterations(int& estim_num_iter) = 0;
00383 virtual void purify_X(const int it) = 0;
00384 virtual void purify_bounds(const int it) = 0;
00385 virtual void save_other_iter_info(IterationInfo& iter_info, int it) = 0;
00386 virtual void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it) = 0;
00387 virtual void return_constant_C(const int it, real& Cval) = 0;
00388
00389
00390 virtual real apply_poly(const int it, real x) = 0;
00391 virtual void apply_poly_to_eigs(const int it, real& homo, real& lumo) = 0;
00392 virtual real compute_derivative(const int it, real x, real& DDf) = 0;
00393
00394
00395
00396
00397
00398 bool initialized_flag;
00399 bool puri_is_prepared_flag;
00400
00401 int use_new_stopping_criterion;
00402 int additional_iterations;
00404 int maxit;
00405 int check_stopping_criterion_iter;
00407 int nocc;
00411 NormType normPuriTrunc;
00415 NormType normPuriStopCrit;
00419 real error_sub;
00420 real error_eig;
00422 real error_per_it;
00424 real constant_C;
00426 real gammaStopEstim;
00431 VectorTypeInt VecPoly;
00436 VectorTypeReal VecGap;
00439 VectorTypeReal ITER_ERROR_VEC;
00440 VectorTypeReal SIGMA_HOMO_VEC, SIGMA_LUMO_VEC;
00441 VectorTypeReal EIG_ABS_GAP_LUMO_VEC, EIG_ABS_GAP_HOMO_VEC;
00442 VectorTypeReal EIG_REL_GAP_LUMO_VEC, EIG_REL_GAP_HOMO_VEC;
00445
00446
00447 int number_of_eigenvalues;
00448
00449 IntervalType homo_bounds;
00450 IntervalType lumo_bounds;
00452 IntervalType homo_bounds_X0;
00453 IntervalType lumo_bounds_X0;
00455 IntervalType homo_bounds_F;
00456 IntervalType lumo_bounds_F;
00458 IntervalType homo_bounds_F_new;
00459 IntervalType lumo_bounds_F_new;
00460
00461
00462 IntervalType spectrum_bounds;
00463 bool computed_spectrum_bounds;
00464
00465
00466
00467 int eigenvectors_method;
00468 int eigenvectors_iterative_method;
00470
00471
00472 real eigensolver_accuracy;
00473
00474 int eigensolver_maxiter;
00475
00476 string eigenvectors_method_str;
00477 string eigenvectors_iterative_method_str;
00478
00479 int use_prev_vector_as_initial_guess;
00480
00481 bool compute_eigenvectors_in_this_SCF_cycle;
00482 bool try_eigv_on_next_iteration_if_fail;
00483
00484 VectorType *eigVecLUMO;
00485 VectorType *eigVecHOMO;
00486
00487 VectorType eigVecLUMORef;
00488 VectorType eigVecHOMORef;
00489
00490
00491 real eigValLUMO;
00492 real eigValHOMO;
00493
00494 int iter_for_homo;
00495 int iter_for_lumo;
00496
00497 VectorTypeInt good_iterations_homo;
00498 VectorTypeInt good_iterations_lumo;
00500 VectorTypeInt really_good_iterations_homo;
00501 VectorTypeInt really_good_iterations_lumo;
00503 int scf_step;
00504
00505 bool compute_eigenvectors_in_each_iteration;
00512 };
00513
00514
00515
00516
00517
00518
00519
00520 template<typename MatrixType>
00521 void PurificationGeneral<MatrixType>::initialize(const MatrixType& F_,
00522 const IntervalType& lumo_bounds_,
00523 const IntervalType& homo_bounds_,
00524 int maxit_,
00525 real error_sub_,
00526 real error_eig_,
00527 int use_new_stopping_criterion_,
00528 NormType norm_truncation_,
00529 NormType norm_stop_crit_,
00530 int nocc_
00531 )
00532 {
00533 X = F_;
00534 maxit = maxit_;
00535 assert(maxit >= 1);
00536 error_sub = error_sub_;
00537 error_eig = error_eig_;
00538 use_new_stopping_criterion = use_new_stopping_criterion_;
00539 normPuriTrunc = norm_truncation_;
00540 normPuriStopCrit = norm_stop_crit_;
00541 nocc = nocc_;
00542
00543 initialized_flag = true;
00544
00545
00546 lumo_bounds_F = lumo_bounds_;
00547 homo_bounds_F = homo_bounds_;
00548
00549
00550
00551 #ifdef ENABLE_PRINTF_OUTPUT
00552 enable_printf_output();
00553 #endif
00554
00555 size_t nnzX;
00556 double nnzXp = get_nnz_X(nnzX);
00557 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Creating purification object: N = %d"
00558 " , nocc = %d , NNZ = %lu <-> %.5lf %%",
00559 X.get_nrows(), nocc, nnzX, nnzXp);
00560
00561
00562 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen norm for the truncation: ");
00563 switch (normPuriTrunc)
00564 {
00565 case mat::mixedNorm:
00566 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "mixed");
00567 break;
00568
00569 case mat::euclNorm:
00570 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "eucl");
00571 break;
00572
00573 case mat::frobNorm:
00574 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "frob");
00575 break;
00576
00577 default:
00578 throw std::runtime_error("Unknown norm in PurificationGeneral");
00579 }
00580
00581
00582 #ifdef USE_CHUNKS_AND_TASKS
00583 if ((normPuriTrunc == mat::mixedNorm) || (normPuriTrunc == mat::euclNorm))
00584 {
00585 normPuriTrunc = mat::frobNorm;
00586 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change norm for the truncation to Frobenius.");
00587 }
00588 #endif
00589
00590 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen norm for the stopping criterion: ");
00591 switch (normPuriStopCrit)
00592 {
00593 case mat::mixedNorm:
00594 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "mixed");
00595 break;
00596
00597 case mat::euclNorm:
00598 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "eucl");
00599 break;
00600
00601 case mat::frobNorm:
00602 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "frob");
00603 break;
00604
00605 default:
00606 throw std::runtime_error("Unknown norm in PurificationGeneral");
00607 }
00608
00609
00610 #ifdef USE_CHUNKS_AND_TASKS
00611 if ((normPuriStopCrit == mat::mixedNorm) || (normPuriStopCrit == mat::euclNorm))
00612 {
00613 normPuriStopCrit = mat::frobNorm;
00614 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change norm the stopping criterion to Frobenius.");
00615 }
00616 #endif
00617
00618 if (this->use_new_stopping_criterion == 1)
00619 {
00620 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen the NEW stopping criterion.");
00621 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Allowed error in subspace %e", (double)error_sub);
00622 }
00623 else
00624 {
00625 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen the OLD stopping criterion.");
00626 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Allowed error in subspace %e , in eigenvalues %e", (double)error_sub, (double)error_eig);
00627 }
00628
00629 check_stopping_criterion_iter = -1;
00630 compute_eigenvectors_in_this_SCF_cycle = false;
00631
00632 set_init_params();
00633
00634 additional_iterations = NUM_ADDITIONAL_ITERATIONS;
00635 info.stopping_criterion = use_new_stopping_criterion;
00636 info.error_subspace = error_sub;
00637
00638 info.debug_output = 0;
00639 #ifdef DEBUG_PURI_OUTPUT
00640 info.debug_output = 1;
00641 #endif
00642 }
00643
00644
00645
00646
00647
00648 template<typename MatrixType>
00649 void PurificationGeneral<MatrixType>::set_eigenvectors_params(string eigenvectors_method_,
00650 string eigenvectors_iterative_method_,
00651 real eigensolver_accuracy_,
00652 int eigensolver_maxiter_,
00653 int scf_step_,
00654 int use_prev_vector_as_initial_guess_,
00655 int try_eigv_on_next_iteration_if_fail_,
00656 VectorType *eigVecLUMO_,
00657 VectorType *eigVecHOMO_
00658 )
00659 {
00660
00661 number_of_eigenvalues = 1;
00662 if (number_of_eigenvalues > 1)
00663 {
00664 throw std::runtime_error("Error in set_eigenvectors_params() : cannot compute more than 1 eigenpair.");
00665 }
00666
00667
00668 eigVecLUMO = eigVecLUMO_;
00669 eigVecHOMO = eigVecHOMO_;
00670
00671 eigensolver_accuracy = eigensolver_accuracy_;
00672 eigensolver_maxiter = eigensolver_maxiter_;
00673 assert(eigensolver_maxiter >= 1);
00674 scf_step = scf_step_;
00675 eigenvectors_method_str = eigenvectors_method_;
00676 eigenvectors_iterative_method_str = eigenvectors_iterative_method_;
00677 eigenvectors_method = get_int_eig_method(eigenvectors_method_);
00678 eigenvectors_iterative_method = get_int_eig_iter_method(eigenvectors_iterative_method_);
00679 use_prev_vector_as_initial_guess = use_prev_vector_as_initial_guess_;
00680 try_eigv_on_next_iteration_if_fail = try_eigv_on_next_iteration_if_fail_;
00681
00682 info.lumo_eigenvector_is_computed = false;
00683 info.homo_eigenvector_is_computed = false;
00684
00685 iter_for_homo = -1;
00686 iter_for_lumo = -1;
00687
00688
00689 if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && (eigenvectors_method == EIG_EMPTY))
00690 {
00691
00692 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "No given method for computing eigenvectors."
00693 "Eigenvectors will not be computed in this SCF cycle. Set eigenvectors to NULL.");
00694
00695 delete eigVecLUMO;
00696 delete eigVecHOMO;
00697 eigVecLUMO = NULL;
00698 eigVecHOMO = NULL;
00699 }
00700 else
00701 {
00702 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen method to compute eigenvectors: %s", eigenvectors_method_str.c_str());
00703 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen iterative method to compute eigenvectors: %s", eigenvectors_iterative_method_str.c_str());
00704 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen eigensolver accuracy: %g", (double)eigensolver_accuracy);
00705 }
00706
00707
00708 if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && use_prev_vector_as_initial_guess)
00709 {
00710 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use eigenvectors from the previous SCF cycle as an initial guess for this SCF cycle");
00711 }
00712
00713
00714 #ifndef USE_CHUNKS_AND_TASKS
00715 if (compute_eigenvectors_in_each_iteration)
00716 {
00717
00718 if ((eigVecLUMO != NULL) && use_prev_vector_as_initial_guess)
00719 {
00720 mat::SizesAndBlocks cols;
00721 if (X.is_empty())
00722 {
00723 throw std::runtime_error("Error in set_eigenvectors_params() : cannot save initial guess for LUMO!");
00724 }
00725 X.getCols(cols);
00726 eigVecLUMORef.resetSizesAndBlocks(cols);
00727 eigVecLUMORef = *eigVecLUMO;
00728 }
00729
00730 if ((eigVecHOMO != NULL) && use_prev_vector_as_initial_guess)
00731 {
00732 mat::SizesAndBlocks cols;
00733 if (X.is_empty())
00734 {
00735 throw std::runtime_error("Error in set_eigenvectors_params() : cannot save initial guess for HOMO!");
00736 }
00737 X.getCols(cols);
00738 eigVecHOMORef.resetSizesAndBlocks(cols);
00739 eigVecHOMORef = *eigVecHOMO;
00740 }
00741 }
00742 #endif
00743 }
00744
00745
00746 template<typename MatrixType>
00747 int PurificationGeneral<MatrixType>::get_int_eig_method(string eigenvectors_method)
00748 {
00749 if (eigenvectors_method == "square")
00750 {
00751 return EIG_SQUARE_INT;
00752 }
00753 if (eigenvectors_method == "projection")
00754 {
00755 return EIG_PROJECTION_INT;
00756 }
00757 if (eigenvectors_method == "")
00758 {
00759 return EIG_EMPTY;
00760 }
00761 throw std::runtime_error("Error in get_int_eig_method(): unknown method to compute eigenvectors");
00762 }
00763
00764
00765 template<typename MatrixType>
00766 int PurificationGeneral<MatrixType>::get_int_eig_iter_method(string eigenvectors_iterative_method)
00767 {
00768 if (eigenvectors_iterative_method == "power")
00769 {
00770 return EIG_POWER_INT;
00771 }
00772 if (eigenvectors_iterative_method == "lanczos")
00773 {
00774 return EIG_LANCZOS_INT;
00775 }
00776 if (eigenvectors_iterative_method == "")
00777 {
00778 return EIG_EMPTY;
00779 }
00780 throw std::runtime_error("Error in get_int_eig_iter_method(): unknown iterative method to compute eigenvectors");
00781 }
00782
00783
00784 template<typename MatrixType>
00785 void PurificationGeneral<MatrixType>::output_norms_and_traces(IterationInfo& iter_info) const
00786 {
00787 real XmX2_fro_norm = iter_info.XmX2_fro_norm;
00788 real XmX2_eucl = iter_info.XmX2_eucl;
00789 real XmX2_mixed_norm = iter_info.XmX2_mixed_norm;
00790 real XmX2_trace = iter_info.XmX2_trace;
00791 int it = iter_info.it;
00792
00793 assert(it >= 0);
00794
00795 #ifndef USE_CHUNKS_AND_TASKS
00796 if (normPuriStopCrit == mat::euclNorm)
00797 {
00798 assert(XmX2_eucl >= 0);
00799 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "||X-X^2||_2 = %e", (double)XmX2_eucl);
00800 }
00801
00802 if (normPuriStopCrit == mat::mixedNorm)
00803 {
00804 assert(XmX2_fro_norm >= 0);
00805 assert(XmX2_mixed_norm >= 0);
00806 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);
00807 }
00808 #endif
00809
00810 if (normPuriStopCrit == mat::frobNorm)
00811 {
00812 assert(XmX2_fro_norm >= 0);
00813 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "||X-X^2||_F = %e", (double)XmX2_fro_norm);
00814 }
00815
00816 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "trace(X-X^2) = %e", (double)XmX2_trace);
00817 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "trace(X) = %e", (double)X.trace());
00818 }
00819
00820
00821
00822
00823
00824 template<typename MatrixType>
00825 void PurificationGeneral<MatrixType>::PurificationStart()
00826 {
00827 Util::TimeMeter total_time;
00828
00829 prepare_to_purification();
00830
00831 #ifdef DEBUG_PURI_OUTPUT
00832 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Starting recursive expansion");
00833 #endif
00834 purification_process();
00835
00836 eigenvalue_bounds_estimation();
00837
00838
00839
00840
00841 if (info.converged == 1)
00842 {
00843 if (compute_eigenvectors_in_this_SCF_cycle && (eigenvectors_method == EIG_PROJECTION_INT))
00844 {
00845 compute_eigenvectors_without_diagonalization_last_iter_proj();
00846 }
00847 }
00848 else
00849 {
00850 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot compute eigenvectors using projection method since the purification did not converge");
00851 }
00852
00853
00854 check_eigenvectors_at_the_end();
00855
00856 total_time.print(LOG_AREA_DENSFROMF, "Recursive expansion");
00857 double total_time_stop = total_time.get_elapsed_wall_seconds();
00858 info.total_time = total_time_stop;
00859 }
00860
00861
00862 template<typename MatrixType>
00863 void PurificationGeneral<MatrixType>::prepare_to_purification()
00864 {
00865 if (!is_initialized())
00866 {
00867 throw std::runtime_error("Error in prepare_to_purification() : function is called for an uninitialized class.");
00868 }
00869
00870 if (!computed_spectrum_bounds)
00871 {
00872 Util::TimeMeter total_time_spectrum_bounds;
00873 compute_spectrum_bounds();
00874 total_time_spectrum_bounds.print(LOG_AREA_DENSFROMF, "compute_spectrum_bounds");
00875 double total_time_spectrum_bounds_stop = total_time_spectrum_bounds.get_elapsed_wall_seconds();
00876 info.time_spectrum_bounds = total_time_spectrum_bounds_stop;
00877 }
00878
00879 info.set_spectrum_bounds(spectrum_bounds.low(), spectrum_bounds.upp());
00880 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Spectrum of F: \t [ %.12lf , %.12lf ]", (double)spectrum_bounds.low(), (double)spectrum_bounds.upp());
00881
00882 map_bounds_to_0_1();
00883 set_truncation_parameters();
00884
00885 if ((eigVecLUMO != NULL) || (eigVecHOMO != NULL))
00886 {
00887
00888 if (1 - homo_bounds.upp() - lumo_bounds.upp() >= TOL_OVERLAPPING_BOUNDS)
00889 {
00890 compute_eigenvectors_in_this_SCF_cycle = true;
00891 info.compute_eigenvectors_in_this_SCF_cycle = compute_eigenvectors_in_this_SCF_cycle;
00892 F = X;
00893
00894
00895 writeToTmpFile(F);
00896 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "calling determine_iteration_for_eigenvectors()");
00897 determine_iteration_for_eigenvectors();
00898 }
00899 else
00900 {
00901 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Homo and lumo inner bounds are too close (< %g), "
00902 "homo and lumo eigenvectors will not be computed", (double)TOL_OVERLAPPING_BOUNDS);
00903 }
00904 }
00905
00906 compute_X();
00907
00908 puri_is_prepared_flag = true;
00909 }
00910
00911
00912 template<typename MatrixType>
00913 void PurificationGeneral<MatrixType>::purification_process()
00914 {
00915 if (!puri_is_prepared())
00916 {
00917 throw std::runtime_error("Error in purification_process() : "
00918 "first expect a call for function prepare_to_purification().");
00919 }
00920
00921 int it;
00922 int stop = 0;
00923 real thresh;
00924 real Xsquare_time_stop = -1, total_time_stop = -1, trunc_time_stop = -1, purify_time_stop = -1;
00925 real frob_diff_time_stop = -1, eucl_diff_time_stop = -1, trace_diff_time_stop = -1, mixed_diff_time_stop = -1, stopping_criterion_time_stop = -1;
00926 int maxit_tmp = maxit;
00927 real estim_order = -1;
00928 real XmX2_trace = -1;
00929 real XmX2_fro_norm = -1;
00930 real XmX2_mixed_norm = -1;
00931 real XmX2_eucl = -1;
00932
00933 int already_stop_before = 0;
00934
00935 info.Iterations.clear();
00936 info.Iterations.reserve(100);
00937
00938
00939 IterationInfo iter_info;
00940
00941
00942 it = 0;
00943
00944
00945 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "BEFORE ITERATIONS:");
00946
00947
00948
00949 #ifdef SAVE_MATRIX_IN_EACH_ITERATION
00950 {
00951 ostringstream str;
00952 str << it;
00953 save_matrix_now(str.str());
00954 }
00955 #endif
00956
00957 Util::TimeMeter total_time;
00958
00959 #ifdef PURI_OUTPUT_NNZ
00960 double nnzX = get_nnz_X();
00961 #endif
00962
00963
00964 Util::TimeMeter trunc_time;
00965 truncate_matrix(thresh, it);
00966 trunc_time.print(LOG_AREA_DENSFROMF, "truncate_matrix");
00967 trunc_time_stop = trunc_time.get_elapsed_wall_seconds();
00968
00969 #ifdef PURI_OUTPUT_NNZ
00970 if((double)thresh >= 0)
00971 {
00972 double nnzXafter = get_nnz_X();
00973 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e , nnz before %.2lf %% , nnz after %.2lf %%", (double)thresh, nnzX, nnzXafter);
00974 }
00975 else
00976 {
00977 double nnzXafter = get_nnz_X();
00978 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "nnz before %.2lf %% , nnz after %.2lf %%", nnzX, nnzXafter);
00979 }
00980 #else
00981 if((double)thresh >= 0)
00982 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e", (double)thresh);
00983 #endif
00984
00985 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before X square");
00986
00987 Util::TimeMeter Xsquare_time;
00988 Xsq = (real)1.0 * X * X;
00989 Xsquare_time.print(LOG_AREA_DENSFROMF, "square");
00990 Xsquare_time_stop = Xsquare_time.get_elapsed_wall_seconds();
00991
00992 #ifdef PURI_OUTPUT_NNZ
00993 size_t nnzXsq;
00994 double nnzXsqp = get_nnz_Xsq(nnzXsq);
00995 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NNZ Xsq = %lu <-> %.2lf %%", nnzXsq, nnzXsqp);
00996 #endif
00997
00998
00999
01000 Util::TimeMeter frob_diff_time;
01001 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
01002 frob_diff_time.print(LOG_AREA_DENSFROMF, "Frobenius norm of X-X^2");
01003 frob_diff_time_stop = frob_diff_time.get_elapsed_wall_seconds();
01004
01005 if (normPuriStopCrit == mat::mixedNorm)
01006 {
01007 #ifndef USE_CHUNKS_AND_TASKS
01008 Util::TimeMeter mixed_diff_time;
01009 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq, mixed_acc);
01010 mixed_diff_time.print(LOG_AREA_DENSFROMF, "Mixed norm of X-X^2");
01011 mixed_diff_time_stop = mixed_diff_time.get_elapsed_wall_seconds();
01012 #endif
01013 }
01014
01015
01016 Util::TimeMeter trace_diff_time;
01017 XmX2_trace = X.trace() - Xsq.trace();
01018 trace_diff_time.print(LOG_AREA_DENSFROMF, "Trace of X-X^2");
01019 trace_diff_time_stop = trace_diff_time.get_elapsed_wall_seconds();
01020
01021 if (normPuriStopCrit == mat::euclNorm)
01022 {
01023 #ifndef USE_CHUNKS_AND_TASKS
01024 Util::TimeMeter eucl_diff_time;
01025 XmX2_eucl = MatrixType::eucl_diff(X, Xsq, eucl_acc);
01026 eucl_diff_time.print(LOG_AREA_DENSFROMF, "Spectral norm of X-X^2");
01027 eucl_diff_time_stop = eucl_diff_time.get_elapsed_wall_seconds();
01028 #endif
01029 }
01030
01031
01032 iter_info.it = it;
01033 iter_info.XmX2_fro_norm = XmX2_fro_norm;
01034 iter_info.XmX2_eucl = XmX2_eucl;
01035 iter_info.XmX2_mixed_norm = XmX2_mixed_norm;
01036 iter_info.XmX2_trace = XmX2_trace;
01037 #ifdef DEBUG_PURI_OUTPUT
01038 output_norms_and_traces(iter_info);
01039 #endif
01040
01041
01042 if (compute_eigenvectors_in_this_SCF_cycle)
01043 {
01044 compute_eigenvectors_without_diagonalization(it, iter_info);
01045 }
01046
01047 ostringstream str_out;
01048 str_out << "Iteration " << it;
01049 total_time.print(LOG_AREA_DENSFROMF, str_out.str().c_str());
01050 total_time_stop = total_time.get_elapsed_wall_seconds();
01051
01052
01053 {
01054 iter_info.gap = 1 - homo_bounds.upp() - lumo_bounds.upp();
01055 iter_info.threshold_X = thresh;
01056 iter_info.Xsquare_time = Xsquare_time_stop;
01057 iter_info.trunc_time = trunc_time_stop;
01058 iter_info.purify_time = 0;
01059 iter_info.NNZ_X = get_nnz_X();
01060 iter_info.NNZ_X2 = get_nnz_Xsq();
01061 iter_info.total_time = total_time_stop;
01062 iter_info.homo_bound_low = homo_bounds.low();
01063 iter_info.lumo_bound_low = lumo_bounds.low();
01064 iter_info.homo_bound_upp = homo_bounds.upp();
01065 iter_info.lumo_bound_upp = lumo_bounds.upp();
01066 stopping_criterion_time_stop = 0;
01067 iter_info.stopping_criterion_time = stopping_criterion_time_stop;
01068 iter_info.eucl_diff_time = eucl_diff_time_stop;
01069 iter_info.frob_diff_time = frob_diff_time_stop;
01070 iter_info.mixed_diff_time = mixed_diff_time_stop;
01071 iter_info.trace_diff_time = trace_diff_time_stop;
01072
01073 save_other_iter_info(iter_info, it);
01074 }
01075
01076
01077 info.Iterations.push_back(iter_info);
01078
01079
01080 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before iteration 1");
01081 Util::TimeMeter timeMeterWriteAndReadAll;
01082 std::string sizesStr = mat::FileWritable::writeAndReadAll();
01083 timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
01084 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
01085
01086
01087 it = 1;
01088 while (it <= maxit_tmp)
01089 {
01090 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "ITERATION %d :", it);
01091
01092 IterationInfo iter_info;
01093
01094 Util::TimeMeter total_time;
01095
01096 Util::TimeMeter purify_time;
01097 purify_X(it);
01098 purify_time.print(LOG_AREA_DENSFROMF, "purify_X");
01099 purify_time_stop = purify_time.get_elapsed_wall_seconds();
01100
01101
01102 #ifdef SAVE_MATRIX_IN_EACH_ITERATION
01103 {
01104 ostringstream str;
01105 str << it;
01106 save_matrix_now(str.str());
01107 }
01108 #endif
01109
01110 #ifdef PURI_OUTPUT_NNZ
01111 double nnzX = get_nnz_X();
01112 #endif
01113
01114
01115 Util::TimeMeter trunc_time;
01116 truncate_matrix(thresh, it);
01117 trunc_time.print(LOG_AREA_DENSFROMF, "truncate_matrix");
01118 trunc_time_stop = trunc_time.get_elapsed_wall_seconds();
01119
01120 #ifdef PURI_OUTPUT_NNZ
01121 double nnzXafter = get_nnz_X();
01122 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e , nnz before %.2lf %% , nnz after %.2lf %%", thresh, nnzX, nnzXafter);
01123 #else
01124 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Actual introduced error %e", thresh);
01125 #endif
01126
01127
01128 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before X square");
01129
01130 Util::TimeMeter Xsquare_time;
01131 Xsq = (real)1.0 * X * X;
01132 Xsquare_time.print(LOG_AREA_DENSFROMF, "square");
01133 Xsquare_time_stop = Xsquare_time.get_elapsed_wall_seconds();
01134 #ifdef PURI_OUTPUT_NNZ
01135 size_t nnzXsq;
01136 double nnzXsqp = get_nnz_Xsq(nnzXsq);
01137 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NNZ Xsq = %lu <-> %.2lf %%", nnzXsq, nnzXsqp);
01138 #endif
01139
01140
01141
01142 purify_bounds(it);
01143
01144
01145 Util::TimeMeter frob_diff_time;
01146 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
01147 frob_diff_time.print(LOG_AREA_DENSFROMF, "Frobenius norm of X-X^2");
01148 frob_diff_time_stop = frob_diff_time.get_elapsed_wall_seconds();
01149
01150 if (normPuriStopCrit == mat::mixedNorm)
01151 {
01152 #ifndef USE_CHUNKS_AND_TASKS
01153 Util::TimeMeter mixed_diff_time;
01154 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq, mixed_acc);
01155
01156 mixed_diff_time.print(LOG_AREA_DENSFROMF, "Mixed norm of X-X^2");
01157 mixed_diff_time_stop = mixed_diff_time.get_elapsed_wall_seconds();
01158 #endif
01159 }
01160
01161 if (normPuriStopCrit == mat::euclNorm)
01162 {
01163 #ifndef USE_CHUNKS_AND_TASKS
01164 Util::TimeMeter eucl_diff_time;
01165 XmX2_eucl = MatrixType::eucl_diff(X, Xsq, eucl_acc);
01166 eucl_diff_time.print(LOG_AREA_DENSFROMF, "Spectral norm of X-X^2");
01167 eucl_diff_time_stop = eucl_diff_time.get_elapsed_wall_seconds();
01168 #endif
01169 }
01170
01171
01172 Util::TimeMeter trace_diff_time;
01173 XmX2_trace = X.trace() - Xsq.trace();
01174 trace_diff_time.print(LOG_AREA_DENSFROMF, "Trace of X-X^2");
01175 trace_diff_time_stop = trace_diff_time.get_elapsed_wall_seconds();
01176
01177
01178
01179 if (XmX2_trace < -1e10)
01180 {
01181 throw std::runtime_error("Error in purification_process() : trace of X-X^2 is negative, seems as a"
01182 " misconvergence of the recursive expansion.");
01183 }
01184
01185 iter_info.it = it;
01186 iter_info.XmX2_fro_norm = XmX2_fro_norm;
01187 iter_info.XmX2_eucl = XmX2_eucl;
01188 iter_info.XmX2_mixed_norm = XmX2_mixed_norm;
01189 iter_info.XmX2_trace = XmX2_trace;
01190 #ifdef DEBUG_PURI_OUTPUT
01191 output_norms_and_traces(iter_info);
01192 #endif
01193
01194
01195
01196 {
01197 iter_info.threshold_X = thresh;
01198 iter_info.Xsquare_time = Xsquare_time_stop;
01199 iter_info.trunc_time = trunc_time_stop;
01200 iter_info.purify_time = purify_time_stop;
01201
01202 iter_info.eucl_diff_time = eucl_diff_time_stop;
01203 iter_info.frob_diff_time = frob_diff_time_stop;
01204 iter_info.mixed_diff_time = mixed_diff_time_stop;
01205 iter_info.trace_diff_time = trace_diff_time_stop;
01206 }
01207
01208
01209
01210 if (compute_eigenvectors_in_this_SCF_cycle)
01211 {
01212 compute_eigenvectors_without_diagonalization(it, iter_info);
01213 }
01214
01215
01216
01217
01218 if (it >= check_stopping_criterion_iter)
01219 {
01220 Util::TimeMeter stopping_criterion_time;
01221 stopping_criterion(iter_info, stop, estim_order);
01222 stopping_criterion_time.print(LOG_AREA_DENSFROMF, "stopping_criterion");
01223 stopping_criterion_time_stop = stopping_criterion_time.get_elapsed_wall_seconds();
01224 iter_info.stopping_criterion_time = stopping_criterion_time_stop;
01225 }
01226 else
01227 {
01228 stop = 0;
01229 estim_order = -1;
01230 }
01231
01232
01233
01234 if ((already_stop_before == 0) && ((stop == 1) || (it == maxit)))
01235 {
01236 if (stop == 1)
01237 {
01238 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "PURIFICATION CONVERGED after %d iterations", it);
01239 info.converged = 1;
01240 }
01241 else
01242 {
01243 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "NOT CONVERGED. Reached the maximum number of iterations %d", maxit);
01244 info.converged = 0;
01245 }
01246
01247 assert(maxit_tmp <= (int)VecPoly.size());
01248 maxit_tmp = it + additional_iterations;
01249 already_stop_before = 1;
01250 }
01251
01252 ostringstream str_out;
01253 str_out << "Iteration " << it;
01254 total_time.print(LOG_AREA_DENSFROMF, str_out.str().c_str());
01255 total_time_stop = total_time.get_elapsed_wall_seconds();
01256
01257
01258
01259
01260
01261
01262 iter_info.NNZ_X = get_nnz_X();
01263 iter_info.NNZ_X2 = get_nnz_Xsq();
01264
01265 iter_info.homo_bound_low = homo_bounds.low();
01266 iter_info.homo_bound_upp = homo_bounds.upp();
01267 iter_info.lumo_bound_low = lumo_bounds.low();
01268 iter_info.lumo_bound_upp = lumo_bounds.upp();
01269
01270 iter_info.total_time = total_time_stop;
01271 iter_info.constantC = constant_C;
01272 if (use_new_stopping_criterion)
01273 {
01274 iter_info.order = estim_order;
01275 }
01276
01277 save_other_iter_info(iter_info, it);
01278
01279
01280
01281 info.Iterations.push_back(iter_info);
01282
01283 it++;
01284 }
01285
01286
01287
01288 double nnzD = get_nnz_X();
01289 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Number of non-zeros in D is %.2lf %%", nnzD);
01290
01291
01292 output_current_memory_usage(LOG_AREA_DENSFROMF, "After the last iteration");
01293 Util::TimeMeter timeMeterWriteAndReadAll_end;
01294 sizesStr = mat::FileWritable::writeAndReadAll();
01295 timeMeterWriteAndReadAll_end.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
01296 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
01297
01298
01299 info.total_it = maxit_tmp;
01300 info.additional_iterations = additional_iterations;
01301
01302 real acc_err_sub = total_subspace_error(maxit_tmp - additional_iterations);
01303 #ifdef DEBUG_PURI_OUTPUT
01304 if (acc_err_sub != -1)
01305 {
01306 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "TOTAL accumulated subspace error is %e", acc_err_sub);
01307 }
01308 #endif
01309 info.accumulated_error_subspace = acc_err_sub;
01310 }
01311
01312
01313
01314
01315 template<typename MatrixType>
01316 void PurificationGeneral<MatrixType>::eigenvalue_bounds_estimation()
01317 {
01318 if (info.converged == 1)
01319 {
01320
01321 VectorTypeReal norms_mixed, norms_frob, traces;
01322
01323 if (normPuriStopCrit == mat::mixedNorm)
01324 {
01325 info.get_vec_mixed_norms(norms_mixed);
01326 }
01327 info.get_vec_frob_norms(norms_frob);
01328 info.get_vec_traces(traces);
01329 get_eigenvalue_estimates(norms_mixed, norms_frob, traces);
01330
01331
01332 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Estimated bounds for the eigenvalues for the Fock matrix:");
01333 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO: [ %.12lf , %.12lf ]", (double)lumo_bounds_F_new.low(), (double)lumo_bounds_F_new.upp());
01334 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO: [ %.12lf , %.12lf ]", (double)homo_bounds_F_new.low(), (double)homo_bounds_F_new.upp());
01335 }
01336 else
01337 {
01338 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot estimate eigenvalues since the purification did not converge");
01339 }
01340
01341
01342 info.homo_estim_low_F = homo_bounds_F_new.low();
01343 info.homo_estim_upp_F = homo_bounds_F_new.upp();
01344 info.lumo_estim_low_F = lumo_bounds_F_new.low();
01345 info.lumo_estim_upp_F = lumo_bounds_F_new.upp();
01346 }
01347
01348
01349
01350
01351
01352 template<typename MatrixType>
01353 void PurificationGeneral<MatrixType>::discard_homo_eigenvector()
01354 {
01355 if ((eigVecHOMO == NULL) || eigVecHOMO->is_empty())
01356 {
01357 return;
01358 }
01359 info.eigValHOMO = 0;
01360 info.homo_eigenvector_is_computed = false;
01361
01362
01363
01364 eigVecHOMO->clear_structure();
01365 }
01366
01367
01368 template<typename MatrixType>
01369 void PurificationGeneral<MatrixType>::discard_lumo_eigenvector()
01370 {
01371 if ((eigVecLUMO == NULL) || eigVecLUMO->is_empty())
01372 {
01373 return;
01374 }
01375 info.eigValLUMO = 0;
01376 info.lumo_eigenvector_is_computed = false;
01377 eigVecLUMO->clear_structure();
01378 }
01379
01380
01381 template<typename MatrixType>
01382 void PurificationGeneral<MatrixType>::check_eigenvectors_at_the_end()
01383 {
01384 if (info.converged != 1)
01385 {
01386 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Discard all computed eigenvectors since the purification did not converge");
01387 discard_homo_eigenvector();
01388 discard_lumo_eigenvector();
01389 return;
01390 }
01391
01392
01393
01394 if (compute_eigenvectors_in_this_SCF_cycle)
01395 {
01396 int ITER_DIFF = 2;
01397
01398 if (!eigVecHOMO->is_empty() && (info.total_it < iter_for_homo))
01399 {
01400 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "HOMO eigenvector was not computed. Iteration for homo: %d, total number of iterations: %d",
01401 iter_for_homo, info.total_it);
01402 discard_homo_eigenvector();
01403 }
01404 else
01405 {
01406
01407 if (!eigVecHOMO->is_empty() && (info.total_it - iter_for_homo < ITER_DIFF) && info.homo_eigenvector_is_computed)
01408 {
01409 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "HOMO eigenvector was computed in one of the last recursive expansion iterations (%d of total %d). "
01410 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
01411 "thus result may not be accurate!", iter_for_homo, info.total_it);
01412 }
01413 }
01414
01415
01416 if (!eigVecLUMO->is_empty() && (info.total_it < iter_for_lumo))
01417 {
01418 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "LUMO eigenvector was not computed. Iteration for lumo: %d, total number of iterations: %d",
01419 iter_for_lumo, info.total_it);
01420 discard_lumo_eigenvector();
01421 }
01422 else
01423 {
01424
01425 if (!eigVecLUMO->is_empty() && (info.total_it - iter_for_lumo < ITER_DIFF) && info.lumo_eigenvector_is_computed)
01426 {
01427 do_output(LOG_CAT_WARNING, LOG_AREA_DENSFROMF, "LUMO eigenvector was computed in one of the last recursive expansion iterations (%d of total %d). "
01428 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
01429 "thus result may not be accurate!", iter_for_lumo, info.total_it);
01430 }
01431 }
01432
01433
01434
01435 real low_lumo_F_bound = info.lumo_estim_low_F;
01436 real upp_lumo_F_bound = info.lumo_estim_upp_F;
01437 real low_homo_F_bound = info.homo_estim_low_F;
01438 real upp_homo_F_bound = info.homo_estim_upp_F;
01439
01440
01441
01442 real flex_tolerance = THRESHOLD_EIG_TOLERANCE;
01443
01444 if (info.homo_eigenvector_is_computed)
01445 {
01446 if ((low_homo_F_bound - flex_tolerance <= eigValHOMO) && (eigValHOMO <= upp_homo_F_bound + flex_tolerance))
01447 {
01448 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO eigenvalue is %lf , HOMO bounds are [ %lf , %lf ]",
01449 (double)eigValHOMO, (double)low_homo_F_bound, (double)upp_homo_F_bound);
01450 info.eigValHOMO = eigValHOMO;
01451 }
01452 else
01453 {
01454 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO eigenvalue is outside HOMO bounds [ %lf , %lf ],"
01455 " discard computed HOMO eigenvector.",
01456 (double)low_homo_F_bound, (double)upp_homo_F_bound);
01457 discard_homo_eigenvector();
01458 }
01459 }
01460 else
01461 {
01462 discard_homo_eigenvector();
01463 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "HOMO eigenvector is not computed.");
01464 }
01465
01466 if (info.lumo_eigenvector_is_computed)
01467 {
01468 if ((low_lumo_F_bound - flex_tolerance <= eigValLUMO) && (eigValLUMO <= upp_lumo_F_bound + flex_tolerance))
01469 {
01470 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO eigenvalue is %lf , LUMO bounds are [ %lf , %lf ]",
01471 (double)eigValLUMO, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
01472 info.eigValLUMO = eigValLUMO;
01473 }
01474 else
01475 {
01476 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed LUMO eigenvalue is outside LUMO bounds [ %lf , %lf ],"
01477 " discard computed LUMO eigenvector.",
01478 (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
01479 discard_lumo_eigenvector();
01480 }
01481 }
01482 else
01483 {
01484 discard_lumo_eigenvector();
01485 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO eigenvector is not computed.");
01486 }
01487
01488 if (info.homo_eigenvector_is_computed && info.lumo_eigenvector_is_computed)
01489 {
01490 real gap = eigValLUMO - eigValHOMO;
01491 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO-LUMO gap is %lf = %lf eV", (double)gap, (double)(gap / UNIT_one_eV));
01492 }
01493 }
01494 }
01495
01496
01497
01498
01499 template<typename MatrixType>
01500 void PurificationGeneral<MatrixType>::compute_X()
01501 {
01502 if (spectrum_bounds.empty())
01503 {
01504 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval spectrum_bounds is empty in compute_X().");
01505 }
01506
01507 #ifdef DEBUG_PURI_OUTPUT
01508 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Put eigenvalues of F to the interval [0,1] in reverse order.");
01509 #endif
01510
01511 real eigmin = spectrum_bounds.low();
01512 real eigmax = spectrum_bounds.upp();
01513 X.add_identity(-eigmax);
01514 X *= ((real)1.0 / (eigmin - eigmax));
01515 }
01516
01517
01518 template<typename MatrixType>
01519 void PurificationGeneral<MatrixType>::map_bounds_to_0_1()
01520 {
01521 if (spectrum_bounds.empty())
01522 {
01523 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval spectrum_bounds is empty in map_bounds_to_0_1().");
01524 }
01525
01526 #ifdef DEBUG_PURI_OUTPUT
01527 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Transform homo and lumo bounds...");
01528 #endif
01529
01530 real eigmin = spectrum_bounds.low();
01531 real eigmax = spectrum_bounds.upp();
01532
01533
01534
01535 homo_bounds = homo_bounds_F;
01536 lumo_bounds = lumo_bounds_F;
01537
01538 homo_bounds.intersect(spectrum_bounds);
01539 lumo_bounds.intersect(spectrum_bounds);
01540
01541 #ifdef DEBUG_PURI_OUTPUT
01542 if (homo_bounds.empty())
01543 {
01544 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
01545 }
01546 if (lumo_bounds.empty())
01547 {
01548 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
01549 }
01550 #endif
01551
01552 if (!mat::Interval<real>::intersect(homo_bounds, lumo_bounds).empty())
01553 {
01554 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Bounds for homo and lumo of F are overlapping.");
01555 }
01556
01557
01558
01559 homo_bounds = (homo_bounds - eigmax) / (eigmin - eigmax);
01560 homo_bounds_X0 = homo_bounds;
01561 homo_bounds = IntervalType(1 - homo_bounds.upp(), 1 - homo_bounds.low());
01562
01563
01564
01565 lumo_bounds = (lumo_bounds - eigmax) / (eigmin - eigmax);
01566 lumo_bounds_X0 = lumo_bounds;
01567
01568 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()));
01569 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "LUMO bounds of X: \t [ %.12lf , %.12lf ]", (double)lumo_bounds.low(), (double)lumo_bounds.upp());
01570
01571 }
01572
01573
01574
01575
01576
01577 template<typename MatrixType>
01578 void PurificationGeneral<MatrixType>::set_truncation_parameters()
01579 {
01580 int estim_num_iter = 0;
01581
01582 estimate_number_of_iterations(estim_num_iter);
01583 info.estim_total_it = estim_num_iter;
01584
01585 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "ESTIMATED NUMBER OF ITERATIONS IS %d", estim_num_iter);
01586
01587 if (estim_num_iter < maxit)
01588 {
01589
01590 maxit = estim_num_iter;
01591 }
01592
01593
01594 error_per_it = error_sub / estim_num_iter;
01595
01596 #ifdef DEBUG_OUTPUT
01597 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "The total allowed subspace error is %e", error_sub);
01598 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Then error due to truncation allowed in each iteration is %e", error_per_it);
01599 #endif
01600 }
01601
01602
01603
01604
01605 template<typename MatrixType>
01606 void PurificationGeneral<MatrixType>::set_spectrum_bounds(real eigmin, real eigmax)
01607 {
01608 spectrum_bounds = IntervalType(eigmin, eigmax);
01609 computed_spectrum_bounds = true;
01610 }
01611
01612
01613
01614
01615 template<typename MatrixType>
01616 void PurificationGeneral<MatrixType>::get_spectrum_bounds(real& eigmin, real& eigmax)
01617 {
01618 if (!computed_spectrum_bounds)
01619 {
01620 compute_spectrum_bounds();
01621 }
01622 eigmin = spectrum_bounds.low();
01623 eigmax = spectrum_bounds.upp();
01624 }
01625
01626
01627
01628
01629 template<typename MatrixType>
01630 void PurificationGeneral<MatrixType>::compute_spectrum_bounds()
01631 {
01632
01633 real eigminG, eigmaxG, eigmin, eigmax;
01634
01635 Util::TimeMeter total_time_spectrum_bounds;
01636 X.gershgorin(eigminG, eigmaxG);
01637 total_time_spectrum_bounds.print(LOG_AREA_DENSFROMF, "gershgorin");
01638
01639
01640 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Gershgorin bounds: [ %.12lf , %.12lf ]", (double)eigminG, (double)eigmaxG);
01641
01642
01643
01644
01645
01646
01647
01648 real smallNumberToExpandWith = template_blas_sqrt(mat::getMachineEpsilon<real>());
01649 eigminG -= smallNumberToExpandWith;
01650 eigmaxG += smallNumberToExpandWith;
01651 #ifdef DEBUG_PURI_OUTPUT
01652 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "EXPANDED Gershgorin bounds: [ %.12lf , %.12lf ]", (double)eigminG, (double)eigmaxG);
01653 #endif
01654
01655 eigmin = eigminG;
01656 eigmax = eigmaxG;
01657
01658
01659 #if 1 // 0 - without Lanczos, 1 - with Lanczos
01660 #ifndef USE_CHUNKS_AND_TASKS
01661
01662
01663 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Trying to impove bounds using Lanczos algorithm...");
01664
01665 real acc = 1e3 * template_blas_sqrt(get_epsilon());
01666 MatrixType Xshifted(X);
01667 Xshifted.add_identity((real)(-1.0) * eigminG);
01668
01669 int maxIter = 100;
01670 try
01671 {
01672 eigmax = Xshifted.eucl(acc, maxIter) + eigminG + acc;
01673 }
01674 catch (mat::AcceptableMaxIter& e)
01675 {
01676
01677 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Lanczos failed to find extreme upper eigenvalue within maxiter... using Gershgorin bound");
01678
01679 eigmax = eigmaxG;
01680 }
01681
01682
01683
01684
01685 Xshifted.add_identity((real)(1.0) * eigminG);
01686 Xshifted.add_identity((real)(-1.0) * eigmaxG);
01687
01688 try
01689 {
01690 eigmin = -Xshifted.eucl(acc, maxIter) + eigmaxG - acc;
01691 }
01692 catch (mat::AcceptableMaxIter& e)
01693 {
01694
01695 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Lanczos failed to find extreme lower eigenvalue within maxiter... using Gershgorin bound");
01696
01697 eigmin = eigminG;
01698 }
01699 #endif // USE_CHUNKS_AND_TASKS
01700 #endif
01701
01702
01703 if (eigmin == eigmax)
01704 {
01705 real tol = 1e-2;
01706 eigmin -= tol;
01707 eigmax += tol;
01708 }
01709
01710 spectrum_bounds = IntervalType(eigmin, eigmax);
01711
01712 computed_spectrum_bounds = true;
01713 }
01714
01715
01716
01717
01718 template<typename MatrixType>
01719 void PurificationGeneral<MatrixType>::stopping_criterion(IterationInfo& iter_info, int& stop, real& estim_order)
01720 {
01721 int it = iter_info.it;
01722 real XmX2_norm_it = -1, XmX2_norm_itm2 = -1;
01723
01724 if (it < check_stopping_criterion_iter)
01725 {
01726 return;
01727 }
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737 if(use_new_stopping_criterion && check_stopping_criterion_iter == -1 && it < 2)
01738 return;
01739 if (use_new_stopping_criterion)
01740 {
01741 assert(it >= 2);
01742
01743 if (normPuriStopCrit == mat::euclNorm)
01744 {
01745 XmX2_norm_it = iter_info.XmX2_eucl;
01746 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_eucl;
01747 }
01748 else
01749 if (normPuriStopCrit == mat::frobNorm)
01750 {
01751 XmX2_norm_it = iter_info.XmX2_fro_norm;
01752 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_fro_norm;
01753 }
01754 else
01755 if (normPuriStopCrit == mat::mixedNorm)
01756 {
01757 XmX2_norm_it = iter_info.XmX2_mixed_norm;
01758 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_mixed_norm;
01759 }
01760 else
01761 {
01762 throw std::runtime_error("Error in stopping_criterion() : unknown matrix norm.");
01763 }
01764
01765 real XmX2_trace = iter_info.XmX2_trace;
01766 check_new_stopping_criterion(it, XmX2_norm_it, XmX2_norm_itm2, XmX2_trace, stop, estim_order);
01767 }
01768 else
01769 {
01770 if (normPuriStopCrit == mat::euclNorm)
01771 {
01772 XmX2_norm_it = iter_info.XmX2_eucl;
01773 }
01774 if (normPuriStopCrit == mat::frobNorm)
01775 {
01776 XmX2_norm_it = iter_info.XmX2_fro_norm;
01777 }
01778 if (normPuriStopCrit == mat::mixedNorm)
01779 {
01780 XmX2_norm_it = iter_info.XmX2_mixed_norm;
01781 }
01782
01783 estim_order = -1;
01784 check_standard_stopping_criterion(XmX2_norm_it, stop);
01785 }
01786 }
01787
01788
01789 template<typename MatrixType>
01790 void PurificationGeneral<MatrixType>::check_standard_stopping_criterion(const real XmX2_norm, int& stop)
01791 {
01792
01793 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Checking standard stopping criterion... ");
01794
01795 bool homoLumo_converged = (homo_bounds.upp() < error_eig &&
01796 lumo_bounds.upp() < error_eig);
01797 bool XmX2norm_converged = XmX2_norm < error_eig;
01798 if (homoLumo_converged || XmX2norm_converged)
01799 {
01800 stop = 1;
01801 }
01802
01803 #ifdef DEBUG_PURI_OUTPUT
01804 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "stop = %d", stop);
01805 #endif
01806 }
01807
01808
01809 template<typename MatrixType>
01810 void PurificationGeneral<MatrixType>::check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace,
01811 int& stop, real& estim_order)
01812 {
01813
01814 if (it < 2)
01815 {
01816 estim_order = -1;
01817 return;
01818 }
01819
01820
01821 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Checking stopping criterion... ");
01822
01823
01824 real C;
01825 return_constant_C(it, C);
01826 this->constant_C = C;
01827 if (C == -1)
01828 {
01829 estim_order = -1;
01830 return;
01831 }
01832
01833 estim_order = template_blas_log(XmX2_norm_it / C) / template_blas_log(XmX2_norm_itm2);
01834
01835 if ((VecPoly[it - 1] != VecPoly[it]) &&
01836 (XmX2_norm_itm2 < 1) &&
01837 (XmX2_norm_it >= C * template_blas_pow(XmX2_norm_itm2, (real)ORDER)))
01838 {
01839 stop = 1;
01840 }
01841
01842
01843 if ((stop != 1) && (XmX2_norm_it < get_epsilon() * template_blas_sqrt(template_blas_sqrt(get_epsilon()))))
01844 {
01845 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "************************************************************************************************************");
01846 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.");
01847 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "************************************************************************************************************");
01848 stop = 1;
01849 }
01850
01851
01852
01853 #ifdef DEBUG_PURI_OUTPUT
01854 if ((VecPoly[it - 1] != VecPoly[it]) && (XmX2_norm_itm2 < 1))
01855 {
01856 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));
01857 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Order of convergence = %lf, stop = %d", (double)estim_order, stop);
01858 }
01859 else
01860 {
01861 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Order of convergence cannot be computed");
01862 }
01863 #endif
01864 }
01865
01866
01867
01868
01869 template<typename MatrixType>
01870 typename PurificationGeneral<MatrixType>::real
01871 PurificationGeneral<MatrixType>::total_subspace_error(int it)
01872 {
01873 assert(it <= (int)VecGap.size());
01874
01875 real error = 0;
01876 real normE;
01877 for (int i = 0; i <= it; ++i)
01878 {
01879 if (VecGap[i] == -1)
01880 {
01881 return -1;
01882 }
01883 normE = info.Iterations[i].threshold_X;
01884 error += normE / (VecGap[i] - normE);
01885 }
01886
01887 return error;
01888 }
01889
01890
01891 template<typename MatrixType>
01892 int PurificationGeneral<MatrixType>::get_est_number_of_puri_iterations()
01893 {
01894 return info.estim_total_it;
01895 }
01896
01897
01898 template<typename MatrixType>
01899 int PurificationGeneral<MatrixType>::get_exact_number_of_puri_iterations()
01900 {
01901 if (info.converged == 1)
01902 {
01903 return info.total_it;
01904 }
01905 else
01906 {
01907 return -1;
01908 }
01909 }
01910
01911
01912
01913 template<typename MatrixType>
01914 void PurificationGeneral<MatrixType>::propagate_values_in_each_iter(real value_unocc, real value_occ,
01915 VectorTypeReal& out_unocc,
01916 VectorTypeReal& out_occ,
01917 int nmax)
01918 {
01919 out_occ.clear();
01920 out_occ.resize(nmax);
01921 out_unocc.clear();
01922 out_unocc.resize(nmax);
01923
01924 out_unocc[0] = value_unocc;
01925 out_occ[0] = value_occ;
01926 real occ, unocc;
01927
01928 for (int i = 1; i < nmax; ++i)
01929 {
01930 occ = out_occ[i - 1];
01931 unocc = out_unocc[i - 1];
01932 apply_poly_to_eigs(i, occ, unocc);
01933 out_occ[i] = occ;
01934 out_unocc[i] = unocc;
01935 }
01936 }
01937
01938
01939 template<typename MatrixType>
01940 void PurificationGeneral<MatrixType>::get_eigenvalue_estimates(const VectorTypeReal& XmX2_norm_mixed,
01941 const VectorTypeReal& XmX2_norm_frob,
01942 const VectorTypeReal& XmX2_trace)
01943 {
01944 estimate_homo_lumo(XmX2_norm_mixed, XmX2_norm_frob, XmX2_trace);
01945
01946 VectorTypeReal h_in, l_in;
01947
01948 real eigmax = spectrum_bounds.upp();
01949 real eigmin = spectrum_bounds.low();
01950 real homo_in_0 = 1 - (homo_bounds_F_new.low() - eigmax) / (eigmin - eigmax);
01951 real lumo_in_0 = (lumo_bounds_F_new.upp() - eigmax) / (eigmin - eigmax);
01952
01953 int n = get_exact_number_of_puri_iterations();
01954 assert(n > 0);
01955
01956 propagate_values_in_each_iter(lumo_in_0, homo_in_0, l_in, h_in, n);
01957
01958 VectorTypeReal YmY2_norm_frob_est;
01959 get_frob_norm_est(XmX2_norm_frob, h_in, l_in, YmY2_norm_frob_est);
01960 }
01961
01962
01963 template<typename MatrixType>
01964 void PurificationGeneral<MatrixType>::get_frob_norm_est(const VectorTypeReal& XmX2_norm_frob,
01965 const VectorTypeReal& h_in,
01966 const VectorTypeReal& l_in,
01967 VectorTypeReal& YmY2_norm_frob_est)
01968 {
01969 int n = get_exact_number_of_puri_iterations();
01970
01971 assert(n > 0);
01972
01973 YmY2_norm_frob_est.clear();
01974 YmY2_norm_frob_est.resize(n);
01975
01976 real th, tl;
01977 for (int i = 0; i < n; ++i)
01978 {
01979 th = h_in[i] - h_in[i] * h_in[i];
01980 tl = l_in[i] - l_in[i] * l_in[i];
01981 YmY2_norm_frob_est[i] = template_blas_sqrt(XmX2_norm_frob[i] * XmX2_norm_frob[i] - th * th - tl * tl);
01982 }
01983 }
01984
01985
01986
01987
01988
01989
01990
01991 template<typename MatrixType>
01992 void PurificationGeneral<MatrixType>::estimate_homo_lumo(const VectorTypeReal& XmX2_norm_mixed,
01993 const VectorTypeReal& XmX2_norm_frob,
01994 const VectorTypeReal& XmX2_trace)
01995 {
01996 homo_bounds_F_new = intervalType(-1e22, 1e22);
01997 lumo_bounds_F_new = intervalType(-1e22, 1e22);
01998
01999
02000 int total_it = info.total_it - info.additional_iterations;
02001
02002
02003 VectorTypeReal bounds_from_it(4);
02004 VectorTypeReal final_bounds(4, 1);
02005
02006
02007 real STOP_NORM = gammaStopEstim - gammaStopEstim * gammaStopEstim;
02008 real vi, wi, mi;
02009 real temp_value;
02010
02011 VectorTypeReal XmX2_norm_out;
02012 if (XmX2_norm_mixed.size() == XmX2_norm_frob.size())
02013 {
02014 XmX2_norm_out = XmX2_norm_mixed;
02015 }
02016 else
02017 {
02018 XmX2_norm_out = XmX2_norm_frob;
02019 }
02020
02021 for (int it = total_it; it >= 0; it--)
02022 {
02023 vi = XmX2_norm_frob[it];
02024 wi = XmX2_trace[it];
02025 mi = XmX2_norm_out[it];
02026
02027 if (vi >= STOP_NORM)
02028 {
02029 break;
02030 }
02031
02032 if (wi <= template_blas_sqrt(get_epsilon()))
02033 {
02034 continue;
02035 }
02036
02037
02038 temp_value = vi * vi / wi;
02039 bounds_from_it[0] = 0.5 * (1 - template_blas_sqrt(1 - 4 * temp_value));
02040 bounds_from_it[1] = 0.5 * (1 - template_blas_sqrt(1 - 4 * mi));
02041
02042
02043 bounds_from_it[2] = bounds_from_it[0];
02044 bounds_from_it[3] = bounds_from_it[1];
02045
02046
02047 apply_inverse_poly_vector(it, bounds_from_it);
02048
02049
02050 final_bounds[0] = std::min(final_bounds[0], bounds_from_it[0]);
02051 final_bounds[1] = std::min(final_bounds[1], bounds_from_it[1]);
02052
02053 final_bounds[2] = std::min(final_bounds[2], bounds_from_it[2]);
02054 final_bounds[3] = std::min(final_bounds[3], bounds_from_it[3]);
02055 }
02056
02057
02058 real maxeig = spectrum_bounds.upp();
02059 real mineig = spectrum_bounds.low();
02060 lumo_bounds_F_new = IntervalType(maxeig * (1 - final_bounds[1]) + mineig * final_bounds[1],
02061 maxeig * (1 - final_bounds[0]) + mineig * final_bounds[0]);
02062 homo_bounds_F_new = IntervalType(mineig * (1 - final_bounds[2]) + maxeig * final_bounds[2],
02063 mineig * (1 - final_bounds[3]) + maxeig * final_bounds[3]);
02064 }
02065
02066
02067
02068
02069 template<typename MatrixType>
02070 void PurificationGeneral<MatrixType>::save_matrix_now(string str)
02071 {
02072 #ifdef USE_CHUNKS_AND_TASKS
02073 vector<int> Itmp, I, Jtmp, J;
02074 vector<real> Vtmp, V;
02075 X.get_all_values(Itmp, Jtmp, Vtmp);
02076
02077 size_t nnz = 0;
02078
02079 for (size_t i = 0; i < Itmp.size(); i++)
02080 {
02081 nnz += (Vtmp[i] != 0);
02082 }
02083
02084 I.reserve(nnz);
02085 J.reserve(nnz);
02086 V.reserve(nnz);
02087
02088 for (size_t i = 0; i < Itmp.size(); i++)
02089 {
02090 if (Vtmp[i] != 0)
02091 {
02092 I.push_back(Itmp[i]);
02093 J.push_back(Jtmp[i]);
02094 V.push_back(Vtmp[i]);
02095 }
02096 }
02097
02098 string name = "X_" + str + ".mtx";
02099 if (write_matrix_to_mtx(name.c_str(), I, J, V, X.get_nrows()) == -1)
02100 {
02101 throw std::runtime_error("Error in save_matrix_now : error in write_matrix_to_mtx.\n");
02102 }
02103 #endif
02104 }
02105
02106
02107
02108
02109
02110
02111
02112
02113 #ifdef USE_CHUNKS_AND_TASKS
02114
02115 template<typename MatrixType>
02116 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F)
02117 {
02118 throw std::runtime_error("compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
02119 }
02120
02121
02122 template<typename MatrixType>
02123 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_last_iter_proj()
02124 {
02125 throw std::runtime_error("compute_eigenvectors_without_diagonalization_last_iter_proj() is not implemented for CHTMatrix.");
02126 }
02127
02128
02129 template<typename MatrixType>
02130 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info)
02131 {
02132 throw std::runtime_error("compute_eigenvectors_without_diagonalization() is not implemented for CHTMatrix.");
02133 }
02134
02135
02136 #else
02137
02138 template<typename MatrixType>
02139 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_on_F(const MatrixType& F, int eigensolver_maxiter_for_F)
02140 {
02141 real start_shift = homo_bounds_F.upp();
02142 real end_shift = lumo_bounds_F.low();
02143 real dist = (end_shift - start_shift) / 15;
02144 real acc_eigv = eigensolver_accuracy;
02145 real eigval;
02146 int eig_num = 1;
02147
02148 MatrixType Fsq;
02149
02150 Fsq = (real)1.0 * F * F;
02151
02152 real sigma;
02153
02154
02155 for (sigma = start_shift; sigma < end_shift; sigma += dist)
02156 {
02157 MatrixType M(Fsq);
02158 M.add_identity(sigma * sigma);
02159 M += ((real) - 2 * sigma) * F;
02160 M = ((real) - 1.0) * M;
02161
02162 vector<real> eigValTmp(1);
02163 vector<int> num_iter_solver(1, -1);
02164
02165 mat::SizesAndBlocks rows;
02166 F.getRows(rows);
02167 vector<VectorType> eigVec(1, VectorType(rows));
02168 eigVec[0].rand();
02169 eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
02170 eigenvectors_iterative_method_str, num_iter_solver,
02171 eigensolver_maxiter_for_F);
02172
02173 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
02174
02175 printf("sigma = %lf , eigval = %lf , iters = %d\n", (double)sigma, (double)eigval, num_iter_solver[0]);
02176 }
02177
02178 sigma = end_shift;
02179 {
02180 MatrixType M(Fsq);
02181 M.add_identity(sigma * sigma);
02182 M += ((real) - 2 * sigma) * F;
02183 M = ((real) - 1.0) * M;
02184
02185 vector<real> eigValTmp(1);
02186 vector<int> num_iter_solver(1, -1);
02187
02188 mat::SizesAndBlocks rows;
02189 F.getRows(rows);
02190 vector<VectorType> eigVec(1, VectorType(rows));
02191 eigVec[0].rand();
02192 eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
02193 eigenvectors_iterative_method_str, num_iter_solver,
02194 eigensolver_maxiter_for_F);
02195
02196 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
02197
02198 printf("sigma = %lf , eigval = %lf , iters = %d\n", (double)sigma, (double)eigval, num_iter_solver[0]);
02199 }
02200 }
02201
02202
02203 template<typename MatrixType>
02204 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization(int it, IterationInfo& iter_info)
02205 {
02206 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
02207
02208
02209 bool compute_homo_now = false;
02210 bool compute_lumo_now = false;
02211
02212
02213 if (compute_eigenvectors_in_each_iteration)
02214 {
02215 if (eigenvectors_method == EIG_SQUARE_INT)
02216 {
02217
02218 if (eigVecHOMO != NULL)
02219 {
02220 for (size_t i = 0; i < good_iterations_homo.size(); ++i)
02221 {
02222 if (good_iterations_homo[i] == it)
02223 {
02224 compute_homo_now = true;
02225 break;
02226 }
02227 }
02228 }
02229
02230
02231 if (eigVecLUMO != NULL)
02232 {
02233 for (size_t i = 0; i < good_iterations_lumo.size(); ++i)
02234 {
02235 if (good_iterations_lumo[i] == it)
02236 {
02237 compute_lumo_now = true;
02238 break;
02239 }
02240 }
02241 }
02242 }
02243 else
02244 {
02245
02246 if (eigVecHOMO != NULL)
02247 {
02248 compute_homo_now = true;
02249 }
02250 if (eigVecLUMO != NULL)
02251 {
02252 compute_lumo_now = true;
02253 }
02254 }
02255 }
02256 else
02257 {
02258
02259 if (eigVecHOMO != NULL)
02260 {
02261 if (it == iter_for_homo)
02262 {
02263 compute_homo_now = true;
02264 }
02265 }
02266 if (eigVecLUMO != NULL)
02267 {
02268 if (it == iter_for_lumo)
02269 {
02270 compute_lumo_now = true;
02271 }
02272 }
02273 }
02274
02275 if (compute_homo_now && !info.homo_eigenvector_is_computed)
02276 {
02277 if (eigenvectors_method == EIG_SQUARE_INT)
02278 {
02279 Util::TimeMeter homo_total_time;
02280
02281 MatrixType M(Xsq);
02282 writeToTmpFile(Xsq);
02283
02284 real sigma = SIGMA_HOMO_VEC[it];
02285 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "choose shift %lf", (double)sigma);
02286
02287
02288
02289 M.add_identity(sigma * sigma);
02290 M += ((real) - 2 * sigma) * X;
02291 M = ((real) - 1.0) * M;
02292
02293 Util::TimeMeter homo_solver_time;
02294 compute_eigenvector(M, eigVecHOMO, it, true);
02295 homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02296 homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
02297
02298 homo_total_time.print(LOG_AREA_DENSFROMF, "compute homo eigenvector");
02299 homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
02300
02301 iter_info.homo_eig_solver_time = homo_solver_time_stop;
02302 iter_info.orbital_homo_time = homo_total_time_stop;
02303
02304
02305 M.clear();
02306 readFromTmpFile(Xsq);
02307 }
02308 else
02309 {
02310 if (compute_eigenvectors_in_each_iteration)
02311 {
02312
02313 ostringstream name;
02314 name << "homo_" << it;
02315 Util::TimeMeter homo_time_save_matrix;
02316 save_matrix_now(name.str());
02317 homo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving homo matrix into mtx");
02318 }
02319 else
02320 {
02321
02322 Util::TimeMeter homo_time_save_matrix;
02323 writeToTmpFile(X);
02324 X_homo = X;
02325 readFromTmpFile(X);
02326 homo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving homo matrix using writeToFile");
02327 }
02328 }
02329 }
02330
02331
02332
02333 if (compute_lumo_now && !info.lumo_eigenvector_is_computed)
02334 {
02335 if (eigenvectors_method == EIG_SQUARE_INT)
02336 {
02337 Util::TimeMeter lumo_total_time;
02338
02339 MatrixType M(Xsq);
02340 writeToTmpFile(Xsq);
02341
02342 real sigma = SIGMA_LUMO_VEC[it];
02343 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "choose shift %lf", (double)sigma);
02344
02345
02346
02347 M.add_identity(sigma * sigma);
02348 M += ((real) - 2 * sigma) * X;
02349 M = ((real) - 1.0) * M;
02350
02351 Util::TimeMeter lumo_solver_time;
02352 compute_eigenvector(M, eigVecLUMO, it, false);
02353 lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02354 lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
02355
02356 lumo_total_time.print(LOG_AREA_DENSFROMF, "compute lumo eigenvector");
02357 lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
02358
02359 iter_info.lumo_eig_solver_time = lumo_solver_time_stop;
02360 iter_info.orbital_lumo_time = lumo_total_time_stop;
02361
02362 M.clear();
02363 readFromTmpFile(Xsq);
02364 }
02365 else
02366 {
02367 if (compute_eigenvectors_in_each_iteration)
02368 {
02369
02370 ostringstream name;
02371 name << "lumo_" << it;
02372 Util::TimeMeter lumo_time_save_matrix;
02373 save_matrix_now(name.str());
02374 lumo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving lumo matrix into mtx");
02375 }
02376 else
02377 {
02378
02379 Util::TimeMeter lumo_time_save_matrix;
02380 writeToTmpFile(X);
02381 X_lumo = X;
02382 readFromTmpFile(X);
02383 lumo_time_save_matrix.print(LOG_AREA_DENSFROMF, "saving lumo matrix using writeToFile");
02384 }
02385 }
02386 }
02387
02388
02389 if (compute_eigenvectors_in_each_iteration)
02390 {
02391
02392 info.homo_eigenvector_is_computed = false;
02393 info.lumo_eigenvector_is_computed = false;
02394 }
02395 }
02396
02397
02398 template<typename MatrixType>
02399 void PurificationGeneral<MatrixType>::compute_eigenvectors_without_diagonalization_last_iter_proj()
02400 {
02401 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
02402 real DX_mult_time_homo_stop, DX_mult_time_lumo_stop;
02403
02404 output_current_memory_usage(LOG_AREA_DENSFROMF, "Before computing eigenvectors:");
02405 Util::TimeMeter timeMeterWriteAndReadAll;
02406 std::string sizesStr = mat::FileWritable::writeAndReadAll();
02407 timeMeterWriteAndReadAll.print(LOG_AREA_DENSFROMF, "FileWritable::writeAndReadAll");
02408 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, ((std::string)"writeAndReadAll sizesStr: '" + sizesStr).c_str());
02409
02410
02411 Xsq.clear();
02412 assert(eigVecHOMO != NULL);
02413 assert(eigVecLUMO != NULL);
02414
02415 if (compute_eigenvectors_in_each_iteration)
02416 {
02417 int total_it = info.total_it;
02418
02419
02420
02421
02422
02423 for (int it = 0; it <= total_it; ++it)
02424 {
02425
02426 MatrixType Xi;
02427 ostringstream name;
02428 name << "X_lumo_" << it << ".mtx";
02429 mat::SizesAndBlocks rows;
02430 mat::SizesAndBlocks cols;
02431 X.getRows(rows);
02432 X.getCols(cols);
02433
02434
02435
02436 Util::TimeMeter homo_total_time;
02437
02438 MatrixType DXi(X);
02439
02440 Util::TimeMeter DX_mult_time_homo;
02441 MatrixType TMP(DXi);
02442
02443 MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, Xi, 0, DXi);
02444 DX_mult_time_homo.print(LOG_AREA_DENSFROMF, "computing D*X");
02445 DX_mult_time_homo_stop = DX_mult_time_homo.get_elapsed_wall_seconds();
02446
02447
02448 MatrixType Zh(X);
02449 Zh -= DXi;
02450
02451 Util::TimeMeter homo_solver_time;
02452 compute_eigenvector(Zh, eigVecHOMO, it, true);
02453 homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02454 homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
02455
02456 homo_total_time.print(LOG_AREA_DENSFROMF, "computing homo eigenvector");
02457 homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
02458
02459 Util::TimeMeter lumo_total_time;
02460
02461
02462 DXi -= Xi;
02463 DXi = (real)(-1) * DXi;
02464
02465 Util::TimeMeter lumo_solver_time;
02466 compute_eigenvector(DXi, eigVecLUMO, it, false);
02467 lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02468 lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
02469
02470 lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
02471 lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
02472
02473 info.Iterations[it].DX_mult_homo_time = DX_mult_time_homo_stop;
02474 info.Iterations[it].DX_mult_lumo_time = 0;
02475
02476 info.Iterations[it].homo_eig_solver_time = homo_solver_time_stop;
02477 info.Iterations[it].lumo_eig_solver_time = lumo_solver_time_stop;
02478
02479 info.Iterations[it].orbital_homo_time = homo_total_time_stop;
02480 info.Iterations[it].orbital_lumo_time = lumo_total_time_stop;
02481 }
02482
02483 output_current_memory_usage(LOG_AREA_DENSFROMF, "After computing eigenvectors:");
02484 return;
02485 }
02486
02487
02488
02489 if (info.total_it >= iter_for_homo)
02490 {
02491
02492 Util::TimeMeter X_homo_read;
02493 readFromTmpFile(X_homo);
02494 X_homo_read.print(LOG_AREA_DENSFROMF, "reading X matrix (for homo) using readFromFile");
02495
02496
02497 Util::TimeMeter homo_total_time;
02498
02499 MatrixType DXi(X);
02500
02501
02502 Util::TimeMeter DX_mult_time_homo;
02503 MatrixType TMP(DXi);
02504
02505 MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, X_homo, 0, DXi);
02506 DX_mult_time_homo.print(LOG_AREA_DENSFROMF, "computing D*X (for homo)");
02507 DX_mult_time_homo_stop = DX_mult_time_homo.get_elapsed_wall_seconds();
02508
02509
02510
02511 MatrixType Zh(X);
02512 Zh -= DXi;
02513 Util::TimeMeter homo_solver_time;
02514 compute_eigenvector(Zh, eigVecHOMO, iter_for_homo, true);
02515 homo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02516 homo_solver_time_stop = homo_solver_time.get_elapsed_wall_seconds();
02517
02518 info.Iterations[iter_for_homo].homo_eig_solver_time = homo_solver_time_stop;
02519 info.Iterations[iter_for_homo].DX_mult_homo_time = DX_mult_time_homo_stop;
02520 Zh.clear();
02521
02522 homo_total_time.print(LOG_AREA_DENSFROMF, "computing homo eigenvector");
02523 homo_total_time_stop = homo_total_time.get_elapsed_wall_seconds();
02524 info.Iterations[iter_for_homo].orbital_homo_time = homo_total_time_stop;
02525
02526
02527
02528 if (iter_for_homo == iter_for_lumo)
02529 {
02530 Util::TimeMeter lumo_total_time;
02531
02532 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Reuse matrix D*X_i for lumo computations");
02533
02534
02535 DXi -= X_homo;
02536 DXi = (real)(-1) * DXi;
02537
02538 Util::TimeMeter lumo_solver_time;
02539 compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo, false);
02540 lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02541 lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
02542
02543
02544 lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
02545 lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
02546
02547 info.Iterations[iter_for_lumo].DX_mult_lumo_time = 0;
02548 info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop;
02549 info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
02550 }
02551 }
02552
02553 X_homo.clear();
02554
02555
02556
02557 if ((info.total_it >= iter_for_lumo) && (iter_for_homo != iter_for_lumo))
02558 {
02559 Util::TimeMeter X_lumo_read;
02560 readFromTmpFile(X_lumo);
02561 X_lumo_read.print(LOG_AREA_DENSFROMF, "reading X matrix (for lumo) using readFromFile");
02562
02563 Util::TimeMeter lumo_total_time;
02564
02565 MatrixType DXi(X);
02566
02567 Util::TimeMeter DX_mult_time_lumo;
02568 MatrixType TMP(DXi);
02569
02570 MatrixType::ssmmUpperTriangleOnly((real)1.0, TMP, X_lumo, 0, DXi);
02571 DX_mult_time_lumo.print(LOG_AREA_DENSFROMF, "computing D*X (for lumo)");
02572 DX_mult_time_lumo_stop = DX_mult_time_lumo.get_elapsed_wall_seconds();
02573
02574
02575 DXi -= X_lumo;
02576 DXi = (real)(-1) * DXi;
02577
02578 Util::TimeMeter lumo_solver_time;
02579 compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo, false);
02580 lumo_solver_time.print(LOG_AREA_DENSFROMF, "compute_eigenvector()");
02581 lumo_solver_time_stop = lumo_solver_time.get_elapsed_wall_seconds();
02582 info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop;
02583 info.Iterations[iter_for_lumo].DX_mult_lumo_time = DX_mult_time_lumo_stop;
02584
02585 lumo_total_time.print(LOG_AREA_DENSFROMF, "computing lumo eigenvector");
02586 lumo_total_time_stop = lumo_total_time.get_elapsed_wall_seconds();
02587 info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
02588
02589 X_lumo.clear();
02590 }
02591 }
02592
02593
02594 #endif // USE_CHUNKS_AND_TASKS
02595
02596
02597
02598 template<typename MatrixType>
02599 void PurificationGeneral<MatrixType>::determine_iteration_for_eigenvectors()
02600 {
02601 if (eigenvectors_method == EIG_SQUARE_INT)
02602 {
02603 get_iterations_for_lumo_and_homo(iter_for_lumo, iter_for_homo);
02604 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvector for HOMO will be computed on the iteration %d. ", iter_for_homo);
02605 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvector for LUMO will be computed on the iteration %d. ", iter_for_lumo);
02606 }
02607 else if (eigenvectors_method == EIG_PROJECTION_INT)
02608 {
02609
02610
02611
02612 iter_for_lumo = 0;
02613 iter_for_homo = 0;
02614 }
02615 else
02616 {
02617 throw std::runtime_error("Error in determine_iteration_for_eigenvectors: unknown method for computing eigenvectors.");
02618 }
02619 }
02620
02621
02627 template<typename MatrixType>
02628 void PurificationGeneral<MatrixType>::get_iterations_for_lumo_and_homo(int& chosen_iter_lumo,
02629 int& chosen_iter_homo)
02630 {
02631 int maxiter = maxit;
02632
02633
02634 real homo0 = 1 - homo_bounds.upp();
02635 real lumo0 = lumo_bounds.upp();
02636 real homoi = homo0, lumoi = lumo0;
02637 real dummy = 0;
02638 real Dh_homo, Dh_lumo, Dgh_homo, Dgh_lumo,
02639 Dgh_homo_max = get_min_double(), Dgh_lumo_max = get_min_double();
02640
02641 chosen_iter_lumo = -1;
02642 chosen_iter_homo = -1;
02643
02644 good_iterations_homo.clear();
02645 good_iterations_lumo.clear();
02646
02647 find_shifts_every_iter();
02648
02649 for (int i = 1; i <= maxiter; ++i)
02650 {
02651 homoi = apply_poly(i, homoi);
02652 lumoi = apply_poly(i, lumoi);
02653
02654 Dh_homo = compute_derivative(i, homo0, dummy);
02655 Dh_lumo = compute_derivative(i, lumo0, dummy);
02656
02657
02658 Dgh_homo = 2 * (homoi - SIGMA_HOMO_VEC[i]) * Dh_homo;
02659 Dgh_lumo = 2 * (lumoi - SIGMA_LUMO_VEC[i]) * Dh_lumo;
02660
02661 if (homoi >= SIGMA_HOMO_VEC[i])
02662 {
02663 good_iterations_homo.push_back(i);
02664 if (Dgh_homo >= Dgh_homo_max)
02665 {
02666 Dgh_homo_max = Dgh_homo;
02667 chosen_iter_homo = i;
02668 }
02669 }
02670
02671 if (lumoi <= SIGMA_LUMO_VEC[i])
02672 {
02673 good_iterations_lumo.push_back(i);
02674 if (template_blas_fabs(Dgh_lumo) >= Dgh_lumo_max)
02675 {
02676 Dgh_lumo_max = template_blas_fabs(Dgh_lumo);
02677 chosen_iter_lumo = i;
02678 }
02679 }
02680 }
02681
02682 if ((chosen_iter_homo == -1) || (chosen_iter_lumo == -1))
02683 {
02684 throw "Error in get_iterations_for_lumo_and_homo() : something went wrong, cannot choose iteration to compute HOMO or LUMO eigenvector.";
02685 }
02686 }
02687
02688
02689 template<typename MatrixType>
02690 void PurificationGeneral<MatrixType>::find_truncation_thresh_every_iter()
02691 {
02692 int maxiter = this->maxit;
02693
02694 this->ITER_ERROR_VEC.clear();
02695 this->ITER_ERROR_VEC.resize(maxiter + 1);
02696
02697 real error = error_per_it;
02698 if (error_per_it == 0)
02699 {
02700 error = this->get_epsilon();
02701 }
02702
02703 for (int i = 1; i <= maxiter; ++i)
02704 {
02705 this->ITER_ERROR_VEC[i] = (error * this->VecGap[i]) / (1 + error);
02706 }
02707 }
02708
02709
02713 template<typename MatrixType>
02714 void PurificationGeneral<MatrixType>::find_shifts_every_iter()
02715 {
02716 int maxiter = maxit;
02717
02718 SIGMA_HOMO_VEC.resize(maxiter + 1);
02719 SIGMA_LUMO_VEC.resize(maxiter + 1);
02720
02721
02722
02723 real homo = 1 - homo_bounds.upp();
02724 real lumo = lumo_bounds.upp();
02725 real homo_out = 1 - homo_bounds.low();
02726 real lumo_out = lumo_bounds.low();
02727
02728 for (int i = 1; i <= maxiter; ++i)
02729 {
02730 homo = apply_poly(i, homo);
02731 lumo = apply_poly(i, lumo);
02732
02733 homo_out = apply_poly(i, homo_out);
02734 lumo_out = apply_poly(i, lumo_out);
02735
02736 SIGMA_HOMO_VEC[i] = (homo_out + lumo) / 2;
02737 SIGMA_LUMO_VEC[i] = (lumo_out + homo) / 2;
02738 }
02739 }
02740
02741
02742 template<typename MatrixType>
02743 void PurificationGeneral<MatrixType>::find_eig_gaps_every_iter()
02744 {
02745 int maxiter = maxit;
02746
02747 EIG_REL_GAP_HOMO_VEC.resize(maxiter + 1);
02748 EIG_REL_GAP_LUMO_VEC.clear();
02749 EIG_REL_GAP_LUMO_VEC.resize(maxiter + 1);
02750 EIG_ABS_GAP_HOMO_VEC.clear();
02751 EIG_ABS_GAP_HOMO_VEC.resize(maxiter + 1);
02752 EIG_ABS_GAP_LUMO_VEC.clear();
02753 EIG_ABS_GAP_LUMO_VEC.resize(maxiter + 1);
02754
02755
02756
02757 real homo0 = homo_bounds_X0.low();
02758 real lumo0 = lumo_bounds_X0.upp();
02759 real one = 1.0;
02760 real zero = 0.0;
02761
02762 real homo_map, lumo_map, one_map, zero_map, sigma;
02763
02764 real homo = homo0, lumo = lumo0;
02765
02766 for (int i = 1; i <= maxiter; ++i)
02767 {
02768 homo = apply_poly(i, homo);
02769 lumo = apply_poly(i, lumo);
02770
02771 sigma = SIGMA_HOMO_VEC[i];
02772
02773 homo_map = (homo - sigma) * (homo - sigma);
02774 lumo_map = (lumo - sigma) * (lumo - sigma);
02775 one_map = (one - sigma) * (one - sigma);
02776 zero_map = (zero - sigma) * (zero - sigma);
02777
02778 EIG_ABS_GAP_HOMO_VEC[i] = min(lumo_map - homo_map, one_map - homo_map);
02779 EIG_REL_GAP_HOMO_VEC[i] = EIG_ABS_GAP_HOMO_VEC[i] /
02780 max(zero_map - homo_map, one_map - homo_map);
02781 }
02782
02783 homo = homo0, lumo = lumo0;
02784
02785 for (int i = 1; i <= maxiter; ++i)
02786 {
02787 homo = apply_poly(i, homo);
02788 lumo = apply_poly(i, lumo);
02789
02790 sigma = SIGMA_LUMO_VEC[i];
02791
02792 homo_map = (homo - sigma) * (homo - sigma);
02793 lumo_map = (lumo - sigma) * (lumo - sigma);
02794 zero_map = (zero - sigma) * (zero - sigma);
02795 one_map = (one - sigma) * (one - sigma);
02796
02797 EIG_ABS_GAP_LUMO_VEC[i] = min(homo_map - lumo_map, zero_map - lumo_map);
02798 EIG_REL_GAP_LUMO_VEC[i] = EIG_ABS_GAP_LUMO_VEC[i] /
02799 max(zero_map - lumo_map, one_map - lumo_map);
02800 }
02801
02802 }
02803
02804
02805 #ifdef USE_CHUNKS_AND_TASKS
02806 template<typename MatrixType>
02807 void PurificationGeneral<MatrixType>::compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
02808 {
02809 throw "compute_eigenvector() is not implemented for CHTMatrix.";
02810 }
02811
02812
02813 #else
02814 template<typename MatrixType>
02815 void PurificationGeneral<MatrixType>::compute_eigenvector(MatrixType const& M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
02816 {
02817 assert(eigVecHOMOorLUMO != NULL);
02818 real acc_eigv = eigensolver_accuracy;
02819
02820 #ifdef DEBUG_PURI_OUTPUT
02821 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Starting compute_eigenvector()");
02822 #endif
02823
02824
02825 if (compute_eigenvectors_in_each_iteration && use_prev_vector_as_initial_guess)
02826 {
02827
02828 if (is_homo)
02829 {
02830 *eigVecHOMO = eigVecHOMORef;
02831 }
02832 else
02833 {
02834 *eigVecLUMO = eigVecLUMORef;
02835 }
02836 }
02837
02838
02839 vector<real> eigValTmp(number_of_eigenvalues);
02840 mat::SizesAndBlocks rows;
02841
02842
02843 X.getRows(rows);
02844
02845
02846
02847
02848
02849 vector<VectorType> eigVec(number_of_eigenvalues, VectorType(rows));
02850 if (use_prev_vector_as_initial_guess)
02851 {
02852 use_prev_vector_as_initial_guess = 0;
02853 if (is_homo)
02854 {
02855 if (!eigVecHOMO->is_empty())
02856 {
02857 eigVec[0] = *eigVecHOMO;
02858 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use HOMO eigenvector computed in the previous SCF cycle as initial guess");
02859 use_prev_vector_as_initial_guess = 1;
02860 }
02861 }
02862 else
02863 {
02864 if (!eigVecLUMO->is_empty())
02865 {
02866 eigVec[0] = *eigVecLUMO;
02867 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Use LUMO eigenvector computed in the previous SCF cycle as initial guess");
02868 use_prev_vector_as_initial_guess = 1;
02869 }
02870 }
02871 }
02872 else
02873 {
02874 eigVec[0].rand();
02875 }
02876
02877
02878 vector<int> num_iter_solver(number_of_eigenvalues, -1);
02879
02880 Util::TimeMeter computeEigenvectors_time;
02881
02882 int eig_num = 1;
02883
02884 eigvec::computeEigenvectors(M, acc_eigv, eigValTmp, eigVec, eig_num,
02885 eigenvectors_iterative_method_str, num_iter_solver,
02886 eigensolver_maxiter);
02887 double eigv_elapsed_wall_sec = computeEigenvectors_time.get_elapsed_wall_seconds();
02888 computeEigenvectors_time.print(LOG_AREA_DENSFROMF, "eigensolver");
02889
02890 if (num_iter_solver.empty())
02891 {
02892 throw std::runtime_error("Error in compute_eigenvector() : (num_iter_solver.empty())");
02893 }
02894
02895
02896
02897 bool is_homo_tmp = false, is_lumo_tmp = false;
02898
02899 if (num_iter_solver[0] == eigensolver_maxiter)
02900 {
02901 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigensolver did not converge within maxiter = %d iterations.", eigensolver_maxiter);
02902 }
02903 else
02904 {
02905 is_homo_tmp = is_homo, is_lumo_tmp = !is_homo;
02906 real eigVal;
02907 get_interval_with_lambda(eigVal, eigVec[0], is_homo_tmp, is_lumo_tmp, it);
02908 if (is_homo_tmp)
02909 {
02910 really_good_iterations_homo.push_back(it);
02911 *eigVecHOMO = eigVec[0];
02912 info.homo_eigenvector_is_computed = true;
02913 info.homo_eigenvector_is_computed_in_iter = it;
02914 info.homo_eigensolver_iter = num_iter_solver[0];
02915 info.homo_eigensolver_time = eigv_elapsed_wall_sec;
02916 eigValHOMO = eigVal;
02917 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() for HOMO in iteration %d "
02918 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
02919 }
02920 else if (is_lumo_tmp)
02921 {
02922 really_good_iterations_lumo.push_back(it);
02923 *eigVecLUMO = eigVec[0];
02924 info.lumo_eigenvector_is_computed = true;
02925 info.lumo_eigenvector_is_computed_in_iter = it;
02926 info.lumo_eigensolver_iter = num_iter_solver[0];
02927 info.lumo_eigensolver_time = eigv_elapsed_wall_sec;
02928 eigValLUMO = eigVal;
02929 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() for LUMO in iteration %d "
02930 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
02931 }
02932 else
02933 {
02934 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "compute_eigenvector() in iteration %d "
02935 ": number of eigensolver iterations is %d", it, num_iter_solver[0]);
02936 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error in compute_eigenvector: wrong eigenvalue ");
02937 }
02938 }
02939
02940 if (compute_eigenvectors_in_each_iteration)
02941 {
02942 save_eigenvectors_to_file(is_homo_tmp, is_lumo_tmp, it);
02943 }
02944 }
02945
02946
02947 #endif // USE_CHUNKS_AND_TASKS
02948
02949
02950 template<typename MatrixType>
02951 void PurificationGeneral<MatrixType>::
02952 writeToTmpFile(MatrixType& A) const
02953 {
02954 #ifndef USE_CHUNKS_AND_TASKS
02955 A.writeToFile();
02956 #endif
02957 }
02958
02959
02960 template<typename MatrixType>
02961 void PurificationGeneral<MatrixType>::
02962 readFromTmpFile(MatrixType& A) const
02963 {
02964 #ifndef USE_CHUNKS_AND_TASKS
02965 A.readFromFile();
02966 #endif
02967 }
02968
02969
02970 template<typename MatrixType>
02971 void PurificationGeneral<MatrixType>::
02972 save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it)
02973 {
02974 if (is_lumo || is_homo)
02975 {
02976 std::vector<real> fullVector;
02977 string eig_name;
02978 if (is_homo)
02979 {
02980 eigVecHOMO->fullvector(fullVector);
02981 eig_name = "homo";
02982 }
02983 else
02984 {
02985 eigVecLUMO->fullvector(fullVector);
02986 eig_name = "lumo";
02987 }
02988
02989
02990 ostringstream name;
02991 if (scf_step != -1)
02992 {
02993 name << eig_name << "_" << it << "_scf_step_" << scf_step << ".txt";
02994 }
02995 else
02996 {
02997 name << eig_name << "_" << it << ".txt";
02998 }
02999 if (write_vector(name.str().c_str(), fullVector) == -1)
03000 {
03001 throw std::runtime_error("Error in save_eigenvectors_to_file() : error in write_vector.");
03002 }
03003 name.str("");
03004 }
03005 }
03006
03007
03008
03009 template<typename MatrixType>
03010 void PurificationGeneral<MatrixType>::get_eigenvalue_of_F_from_eigv_of_Xi(real& eigVal, const VectorType& eigVec)
03011 {
03012 readFromTmpFile(F);
03013 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
03014 writeToTmpFile(F);
03015 eigVal = approx_eig;
03016 }
03017
03018
03019 template<typename MatrixType>
03020 void PurificationGeneral<MatrixType>::get_interval_with_lambda(real& eigVal, VectorType& eigVec, bool& is_homo, bool& is_lumo, const int iter)
03021 {
03022 assert(is_homo || is_lumo);
03023
03024 bool is_homo_init = is_homo;
03025 bool is_lumo_init = is_lumo;
03026
03027
03028
03029
03030
03031
03032
03033
03034
03035
03036
03037 real low_lumo_F_bound, low_homo_F_bound;
03038 real upp_lumo_F_bound, upp_homo_F_bound;
03039
03040 if (eigenvectors_method == EIG_SQUARE_INT)
03041
03042 {
03043 low_lumo_F_bound = lumo_bounds_F.low();
03044 upp_lumo_F_bound = lumo_bounds_F.upp();
03045 low_homo_F_bound = homo_bounds_F.low();
03046 upp_homo_F_bound = homo_bounds_F.upp();
03047 }
03048 else if (eigenvectors_method == EIG_PROJECTION_INT)
03049
03050 {
03051 low_lumo_F_bound = info.lumo_estim_low_F;
03052 upp_lumo_F_bound = info.lumo_estim_upp_F;
03053 low_homo_F_bound = info.homo_estim_low_F;
03054 upp_homo_F_bound = info.homo_estim_upp_F;
03055 }
03056 else
03057 {
03058 throw std::runtime_error("Error in get_interval_with_lambda() : unexpected eigenvectors_method value.");
03059 }
03060
03061 #ifdef DEBUG_PURI_OUTPUT
03062 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Check Rayleigh quotient...");
03063 #endif
03064
03065 readFromTmpFile(F);
03066 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
03067 writeToTmpFile(F);
03068
03069 eigVal = approx_eig;
03070
03071 real flex_tolerance = THRESHOLD_EIG_TOLERANCE;
03072
03073
03074 if ((approx_eig <= upp_homo_F_bound + flex_tolerance) && (approx_eig >= low_homo_F_bound - flex_tolerance))
03075 {
03076 is_homo = true;
03077 is_lumo = false;
03078
03079 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed HOMO eigenvalue of F is %lf, "
03080 "HOMO bounds are [ %lf , %lf ]", (double)approx_eig, (double)low_homo_F_bound, (double)upp_homo_F_bound);
03081
03082 iter_for_homo = iter;
03083
03084 if (is_lumo_init && (eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
03085 {
03086 iter_for_lumo = iter_for_lumo + 1;
03087 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute LUMO in the next iteration.");
03088 }
03089 }
03090
03091 else if ((approx_eig <= upp_lumo_F_bound + flex_tolerance) && (approx_eig >= low_lumo_F_bound - flex_tolerance))
03092 {
03093 is_homo = false;
03094 is_lumo = true;
03095
03096 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Computed LUMO eigenvalue of F is %lf, "
03097 "LUMO interval [ %lf , %lf ]", (double)approx_eig, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
03098
03099 iter_for_lumo = iter;
03100
03101 if (is_homo_init && (eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
03102 {
03103 iter_for_homo = iter_for_homo + 1;
03104 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute HOMO in the next iteration.");
03105 }
03106 }
03107 else
03108 {
03109 is_homo = false;
03110 is_lumo = false;
03111
03112 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvalue is outside of both intervals for homo and lumo.");
03113 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Eigenvalue is %lf, HOMO interval [ %lf , %lf ], LUMO interval [ %lf , %lf ]",
03114 (double)approx_eig, (double)low_homo_F_bound, (double)upp_homo_F_bound, (double)low_lumo_F_bound, (double)upp_lumo_F_bound);
03115
03116
03117
03118
03119 if ((eigenvectors_method == EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
03120 {
03121 iter_for_homo = iter_for_homo + 1;
03122 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "We will try to compute HOMO (or LUMO) in the next iteration.");
03123 }
03124 }
03125 }
03126
03127
03128
03129
03130
03131
03132
03133 template<typename MatrixType>
03134 void PurificationGeneral<MatrixType>::gen_matlab_file_norm_diff(const char *filename) const
03135 {
03136 std::ofstream f;
03137 f.open(filename, std::ios::out);
03138 if (!f.is_open())
03139 {
03140 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03141 return;
03142 }
03143
03144 int it = info.total_it;
03145
03146 f << "POLY = [";
03147 for (int i = 0; i <= it; ++i)
03148 {
03149 f << VecPoly[i] << " ";
03150 }
03151 f << "];" << std::endl;
03152
03153
03154 if (normPuriStopCrit == mat::euclNorm)
03155 {
03156 f << "X_norm = [";
03157 for (int i = 0; i <= it; ++i)
03158 {
03159 f << (double)info.Iterations[i].XmX2_eucl << " ";
03160 }
03161 f << "];" << std::endl;
03162 f << " norm_letter = '2';" << std::endl;
03163 }
03164 else if (normPuriStopCrit == mat::frobNorm)
03165 {
03166 f << "X_norm = [";
03167 for (int i = 0; i <= it; ++i)
03168 {
03169 f << (double)info.Iterations[i].XmX2_fro_norm << " ";
03170 }
03171 f << "];" << std::endl;
03172 f << " norm_letter = 'F';" << std::endl;
03173 }
03174 else if (normPuriStopCrit == mat::mixedNorm)
03175 {
03176 f << "X_norm = [";
03177 for (int i = 0; i <= it; ++i)
03178 {
03179 f << (double)info.Iterations[i].XmX2_mixed_norm << " ";
03180 }
03181 f << "];" << std::endl;
03182 f << " norm_letter = 'M';" << std::endl;
03183 }
03184 else
03185 {
03186 throw "Wrong norm in PurificationGeneral::gen_matlab_file_norm_diff()";
03187 }
03188
03189 f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
03190 f << "it = " << it << ";" << std::endl;
03191 f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
03192 f << "fighandle = figure; clf;" << std::endl;
03193 f << "MARKER = ['o', '+'];" << std::endl;
03194 f << "semilogy(0:stop_iteration, X_norm(1:stop_iteration+1), '-b', plot_props{:});" << std::endl;
03195 f << "hold on" << std::endl;
03196 f << "for i = 1:stop_iteration+1" << std::endl;
03197 f << " if POLY(i) == 1" << std::endl;
03198 f << " h1 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
03199 f << " else" << std::endl;
03200 f << " h2 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
03201 f << " end" << std::endl;
03202 f << "end" << std::endl;
03203 f << "if stop_iteration ~= it" << std::endl;
03204 f << "h3 = semilogy(stop_iteration+1:it, X_norm(stop_iteration+2:it+1), '-.vr', plot_props{:});" << std::endl;
03205 f << "semilogy(stop_iteration:stop_iteration+1, X_norm(stop_iteration+1:stop_iteration+2), '-.r', plot_props{:});" << std::endl;
03206 f << "legend([h1 h2 h3],{'$x^2$', '$2x-x^2$', 'After stop'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
03207 f << "else" << std::endl;
03208 f << "legend([h1 h2],{'$x^2$', '$2x-x^2$'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
03209 f << "end" << std::endl;
03210 f << "xlabel('Iteration SP2', 'Interpreter','latex');" << std::endl;
03211 f << "ylabel({['$\\|X_i-X_i^2\\|_{' norm_letter '}$']},'interpreter','latex');" << std::endl;
03212 f << "grid on" << std::endl;
03213 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03214 f << "set(gca,'FontSize',20);" << std::endl;
03215 f << "xlim([0 it]);" << std::endl;
03216 f << "ylim([-inf inf]);" << std::endl;
03217 f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
03218 f << "a = 16; S = logspace(-a, 1, a+2);" << std::endl;
03219 f << "set(gca,'YTick',S(1:2:end));" << std::endl;
03220
03221 f << "hold off" << std::endl;
03222
03223 f << "% print( fighandle, '-depsc2', 'norm_diff.eps');" << std::endl;
03224
03225 f.close();
03226 }
03227
03228
03229 template<typename MatrixType>
03230 void PurificationGeneral<MatrixType>::gen_matlab_file_threshold(const char *filename) const
03231 {
03232 std::ofstream f;
03233 f.open(filename, std::ios::out);
03234 if (!f.is_open())
03235 {
03236 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03237 return;
03238 }
03239
03240 int it = info.total_it;
03241 f << "Thresh = [";
03242
03243 for (int i = 0; i <= it; ++i)
03244 {
03245 f << (double)info.Iterations[i].threshold_X << " ";
03246 }
03247 f << "];" << std::endl;
03248
03249 f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
03250 f << "it = " << it << ";" << std::endl;
03251 f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
03252 f << "fighandle = figure; clf;" << std::endl;
03253 f << "semilogy(0:stop_iteration, Thresh(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
03254 f << "hold on" << std::endl;
03255 f << "if stop_iteration ~= it" << std::endl;
03256 f << "semilogy(stop_iteration+1:it, Thresh(stop_iteration+2:it+1), '-^r', plot_props{:});" << std::endl;
03257 f << "semilogy(stop_iteration:stop_iteration+1, Thresh(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
03258 f << "legend('before stop', 'after stop', 'Location','NorthWest');" << std::endl;
03259 f << "end" << std::endl;
03260 f << "grid on" << std::endl;
03261 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03262 f << "set(gca,'FontSize',20);" << std::endl;
03263 f << "xlim([0 it]);" << std::endl;
03264 f << "ylim([-inf inf]);" << std::endl;
03265 f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
03266 f << "hold off" << std::endl;
03267 f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
03268 f << "ylabel('Threshold value', 'interpreter','latex');" << std::endl;
03269 f << "% print( fighandle, '-depsc2', 'threshold.eps');" << std::endl;
03270
03271 f.close();
03272 }
03273
03274
03275 template<typename MatrixType>
03276 void PurificationGeneral<MatrixType>::gen_matlab_file_nnz(const char *filename) const
03277 {
03278 std::ofstream f;
03279 f.open(filename, std::ios::out);
03280 if (!f.is_open())
03281 {
03282 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03283 return;
03284 }
03285
03286 int it = info.total_it;
03287 f << "NNZ_X = [";
03288
03289 for (int i = 0; i <= it; ++i)
03290 {
03291 f << (double)info.Iterations[i].NNZ_X << " ";
03292 }
03293 f << "];" << std::endl;
03294
03295 f << "NNZ_X2 = [";
03296
03297 for (int i = 0; i <= it; ++i)
03298 {
03299 f << (double)info.Iterations[i].NNZ_X2 << " ";
03300 }
03301 f << "];" << std::endl;
03302
03303
03304
03305 f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
03306 f << "it = " << it << ";" << std::endl;
03307 f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
03308 f << "fighandle = figure; clf;" << std::endl;
03309 f << "h2 = plot(0:stop_iteration, NNZ_X(1:stop_iteration+1), '-ob', plot_props{:});" << std::endl;
03310 f << "hold on" << std::endl;
03311 f << "h1 = plot(0:stop_iteration, NNZ_X2(1:stop_iteration+1), '-vm', plot_props{:});" << std::endl;
03312 f << "if stop_iteration ~= it" << std::endl;
03313 f << "plot(stop_iteration+1:it, NNZ_X(stop_iteration+2:it+1), '-vr', plot_props{:});" << std::endl;
03314 f << "plot(stop_iteration+1:it, NNZ_X2(stop_iteration+2:it+1), '-*r', plot_props{:});" << std::endl;
03315 f << "plot(stop_iteration:stop_iteration+1, NNZ_X(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
03316 f << "plot(stop_iteration:stop_iteration+1, NNZ_X2(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
03317 f << "end" << std::endl;
03318 f << "legend([h1, h2], {'$X^2$', '$X$'}, 'interpreter','latex');" << std::endl;
03319 f << "grid on" << std::endl;
03320 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03321 f << "set(gca,'FontSize',20);" << std::endl;
03322 f << "xlim([0 it]);" << std::endl;
03323 f << "ylim([0 inf]);" << std::endl;
03324 f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
03325 f << "hold off" << std::endl;
03326 f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
03327 f << "ylabel('NNZ [\\%]', 'interpreter','latex');" << std::endl;
03328 f << "% print( fighandle, '-depsc2', 'nnz.eps');" << std::endl;
03329
03330 f.close();
03331 }
03332
03333
03334 template<typename MatrixType>
03335 void PurificationGeneral<MatrixType>::gen_matlab_file_cond_num(const char *filename) const
03336 {
03337 std::ofstream f;
03338 f.open(filename, std::ios::out);
03339 if (!f.is_open())
03340 {
03341 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03342 return;
03343 }
03344
03345 int it = info.total_it;
03346 f << "in= [";
03347
03348 for (int i = 0; i <= it; ++i)
03349 {
03350 f << (double)(1 - info.Iterations[i].homo_bound_upp - info.Iterations[i].lumo_bound_upp) << " ";
03351 }
03352 f << "];" << std::endl;
03353
03354 f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
03355 f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
03356 f << "fighandle = figure; clf;" << std::endl;
03357 f << "plot(0:stop_iteration, 1./in(1:stop_iteration+1), '-*r', plot_props{:});" << std::endl;
03358 f << "grid on" << std::endl;
03359 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03360 f << "set(gca,'FontSize',20);" << std::endl;
03361 f << "xlim([0 it]);" << std::endl;
03362 f << "ylim([-inf inf]);" << std::endl;
03363 f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
03364 f << "hold off" << std::endl;
03365 f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
03366 f << "ylabel('$\\kappa$', 'interpreter','latex');" << std::endl;
03367 f << "% print( fighandle, '-depsc2', 'cond.eps');" << std::endl;
03368
03369 f.close();
03370 }
03371
03372
03373 template<typename MatrixType>
03374 void PurificationGeneral<MatrixType>::gen_matlab_file_eigs(const char *filename) const
03375 {
03376 std::ofstream f;
03377 f.open(filename, std::ios::out);
03378 if (!f.is_open())
03379 {
03380 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03381 return;
03382 }
03383
03384 int it = info.total_it;
03385 f << "homo_low= [";
03386
03387 for (int i = 0; i <= it; ++i)
03388 {
03389 f << (double)info.Iterations[i].homo_bound_low << " ";
03390 }
03391 f << "];" << std::endl;
03392
03393 f << "homo_upp= [";
03394
03395 for (int i = 0; i <= it; ++i)
03396 {
03397 f << (double)info.Iterations[i].homo_bound_upp << " ";
03398 }
03399 f << "];" << std::endl;
03400
03401 f << "lumo_low= [";
03402
03403 for (int i = 0; i <= it; ++i)
03404 {
03405 f << (double)info.Iterations[i].lumo_bound_low << " ";
03406 }
03407 f << "];" << std::endl;
03408
03409 f << "lumo_upp= [";
03410
03411 for (int i = 0; i <= it; ++i)
03412 {
03413 f << (double)info.Iterations[i].lumo_bound_upp << " ";
03414 }
03415 f << "];" << std::endl;
03416
03417
03418
03419 f << "stop_iteration = " << it - info.additional_iterations << ";" << std::endl;
03420 f << "plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
03421 f << "fighandle = figure; clf;" << std::endl;
03422 f << "semilogy(0:stop_iteration, homo_upp(1:stop_iteration+1), '-^b', plot_props{:});" << std::endl;
03423 f << "hold on" << std::endl;
03424 f << "semilogy(0:stop_iteration, homo_low(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
03425 f << "semilogy(0:stop_iteration, lumo_low(1:stop_iteration+1), '-vr', plot_props{:});" << std::endl;
03426 f << "semilogy(0:stop_iteration, lumo_upp(1:stop_iteration+1), '-^r', plot_props{:});" << std::endl;
03427 f << "grid on" << std::endl;
03428 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03429 f << "set(gca,'FontSize',20);" << std::endl;
03430 f << "xlim([0 stop_iteration]);" << std::endl;
03431 f << "set(gca,'XTick',[0 5:5:stop_iteration]);" << std::endl;
03432 f << "ylim([-inf inf]);" << std::endl;
03433 f << "hold off" << std::endl;
03434 f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
03435 f << "ylabel('Homo and lumo bounds', 'interpreter','latex');" << std::endl;
03436 f << "% print( fighandle, '-depsc2', 'eigs.eps');" << std::endl;
03437
03438 f.close();
03439 }
03440
03441
03442 template<typename MatrixType>
03443 void PurificationGeneral<MatrixType>::gen_matlab_file_time(const char *filename) const
03444 {
03445 std::ofstream f;
03446 f.open(filename, std::ios::out);
03447 if (!f.is_open())
03448 {
03449 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03450 return;
03451 }
03452
03453 int it = info.total_it;
03454
03455 f << "time_total = [";
03456 for (int i = 0; i <= it; ++i)
03457 {
03458 f << std::setprecision(16) << (double)info.Iterations[i].total_time << " ";
03459 }
03460 f << "];" << std::endl;
03461
03462 f << "time_square = [";
03463 for (int i = 0; i <= it; ++i)
03464 {
03465 f << std::setprecision(16) << (double)info.Iterations[i].Xsquare_time << " ";
03466 }
03467 f << "];" << std::endl;
03468
03469 f << "time_trunc = [";
03470 for (int i = 0; i <= it; ++i)
03471 {
03472 f << std::setprecision(16) << (double)info.Iterations[i].trunc_time << " ";
03473 }
03474 f << "];" << std::endl;
03475
03476 if (info.compute_eigenvectors_in_this_SCF_cycle)
03477 {
03478 f << "time_eigenvectors_homo = [";
03479 for (int i = 0; i <= it; ++i)
03480 {
03481 f << std::setprecision(16) << (double)info.Iterations[i].orbital_homo_time << " ";
03482 }
03483 f << "];" << std::endl;
03484 f << "time_eigenvectors_lumo = [";
03485 for (int i = 0; i <= it; ++i)
03486 {
03487 f << std::setprecision(16) << (double)info.Iterations[i].orbital_lumo_time << " ";
03488 }
03489 f << "];" << std::endl;
03490
03491 f << "time_solver_homo = [";
03492 for (int i = 0; i <= it; ++i)
03493 {
03494 f << std::setprecision(16) << (double)info.Iterations[i].homo_eig_solver_time << " ";
03495 }
03496 f << "];" << std::endl;
03497 f << "time_solver_lumo = [";
03498 for (int i = 0; i <= it; ++i)
03499 {
03500 f << std::setprecision(16) << (double)info.Iterations[i].lumo_eig_solver_time << " ";
03501 }
03502 f << "];" << std::endl;
03503
03504
03505 if (eigenvectors_method == EIG_SQUARE_INT)
03506 {
03507
03508 f << "X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_eigenvectors_homo-time_solver_homo;"
03509 " time_eigenvectors_lumo-time_solver_lumo; "
03510 " time_total - time_square - time_trunc - time_eigenvectors_homo - time_eigenvectors_lumo];" << std::endl;
03511 }
03512 else
03513 {
03514 f << "time_DX_homo = [";
03515 for (int i = 0; i <= it; ++i)
03516 {
03517 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_homo_time << " ";
03518 }
03519 f << "];" << std::endl;
03520 f << "time_DX_lumo = [";
03521 for (int i = 0; i <= it; ++i)
03522 {
03523 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_lumo_time << " ";
03524 }
03525 f << "];" << std::endl;
03526
03527
03528
03529 f << "X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_DX_homo; time_DX_lumo;"
03530 " time_eigenvectors_homo - time_DX_homo - time_solver_homo; time_eigenvectors_lumo - time_DX_lumo - time_solver_lumo;"
03531 " time_total - time_square - time_trunc];" << std::endl;
03532 }
03533 }
03534 else
03535 {
03536 f << "X = [time_square; time_trunc; time_total - time_square - time_trunc];" << std::endl;
03537 }
03538
03539 f << "it = " << it << ";" << std::endl;
03540 f << "xtick = 0:it;" << std::endl;
03541 f << "fighandle = figure; clf;" << std::endl;
03542 f << "b=bar(xtick, X', 'stacked');" << std::endl;
03543 f << "axis tight;" << std::endl;
03544 f << "grid on" << std::endl;
03545 f << "set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
03546 f << "set(gca,'FontSize',20);" << std::endl;
03547 f << "xlim([0 it]);" << std::endl;
03548 f << "set(gca,'XTick',[0 5:5:it]);" << std::endl;
03549 f << "ylim([-inf inf]);" << std::endl;
03550 f << "xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
03551 f << "ylabel('Time [s]', 'interpreter','latex');" << std::endl;
03552 if (info.compute_eigenvectors_in_this_SCF_cycle)
03553 {
03554
03555
03556
03557 if (eigenvectors_method == EIG_SQUARE_INT)
03558 {
03559 f << "legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
03560 }
03561 else
03562 {
03563 f << "legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'DX (lumo)', 'DX (homo)', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
03564 }
03565 }
03566 else
03567 {
03568 f << "legend(flipud(b(:)), {'other', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
03569 }
03570 f << "% print( fighandle, '-depsc2', 'time.eps');" << std::endl;
03571
03572 f.close();
03573 }
03574
03575
03576
03577
03578
03579
03580
03581 template<typename MatrixType>
03582 void PurificationGeneral<MatrixType>::gen_python_file_nnz(const char *filename) const
03583 {
03584 std::ofstream f;
03585 f.open(filename, std::ios::out);
03586 if (!f.is_open())
03587 {
03588 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Error: cannot open file");
03589 return;
03590 }
03591
03592 f << "import numpy as np" << std::endl;
03593 f << "import pylab as pl" << std::endl;
03594 f << "import matplotlib.font_manager as font_manager" << std::endl;
03595 f << "from matplotlib import rc" << std::endl;
03596 f << "rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})" << std::endl;
03597 f << "rc('text', usetex=True)" << std::endl;
03598 f << std::endl;
03599
03600 int it = info.total_it;
03601 f << "NNZ_X = np.array([";
03602
03603 for (int i = 0; i <= it; ++i)
03604 {
03605 f << (double)info.Iterations[i].NNZ_X << ", ";
03606 }
03607 f << "])" << std::endl;
03608
03609 f << "NNZ_X2 = np.array([";
03610
03611 for (int i = 0; i <= it; ++i)
03612 {
03613 f << (double)info.Iterations[i].NNZ_X2 << ", ";
03614 }
03615 f << "])" << std::endl;
03616
03617 f << "stop_iteration = " << it - info.additional_iterations << std::endl;
03618 f << "it = " << it << std::endl;
03619 f << "font_prop = font_manager.FontProperties(size=20)" << std::endl;
03620 f << "prop = {'markersize':8, 'fillstyle':'none', 'linewidth':2, 'markeredgewidth':2}" << std::endl;
03621 f << "fig1 = pl.figure(figsize = (8, 6), num='nnz')" << std::endl;
03622 f << "p1, = pl.plot(range(0, stop_iteration+1), NNZ_X2[0:stop_iteration+1], '-vm', **prop);" << std::endl;
03623 f << "p2, = pl.plot(range(0, stop_iteration+1), NNZ_X[0:stop_iteration+1], '-ob', **prop);" << std::endl;
03624 f << "if stop_iteration != it:" << std::endl;
03625 f << " pl.plot(range(stop_iteration+1,it+1), NNZ_X[stop_iteration+1:it+1], '-vr', **prop);" << std::endl;
03626 f << " pl.plot(range(stop_iteration+1,it+1), NNZ_X2[stop_iteration+1:it+1], '-*r', **prop);" << std::endl;
03627 f << " pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
03628 f << " pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X2[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
03629 f << std::endl;
03630 f << "pl.xlim(0, it);" << std::endl;
03631 f << "pl.ylim(ymin=0);" << std::endl;
03632 f << "pl.xlabel('Iteration SP2', fontproperties=font_prop);" << std::endl;
03633 f << "pl.ylabel('NNZ [\\%]', fontproperties=font_prop); " << std::endl;
03634 f << "pl.legend((p1, p2), ('$X^2$', '$X$'), loc='best', numpoints=1)" << std::endl;
03635 f << "pl.xticks(range(5, it, 5))" << std::endl;
03636 f << "pl.tick_params(labelsize=20)" << std::endl;
03637 f << "pl.grid(which='major', alpha=0.5, color='k', linestyle='--', linewidth=1) " << std::endl;
03638 f << "pl.savefig('nnz.eps', format='eps', dpi=1000)" << std::endl;
03639 f << "pl.show()" << std::endl;
03640
03641 f.close();
03642 }
03643
03644
03645 #endif //HEADER_PURIFICATION_GENERAL