39#ifndef HEADER_PURIFICATION_GENERAL
40#define HEADER_PURIFICATION_GENERAL
76#define NUM_ADDITIONAL_ITERATIONS 0
79#define DEBUG_PURI_OUTPUT
80#define PURI_OUTPUT_NNZ
115template<
typename MatrixType>
159 bool use_new_stopping_criterion_,
197 {
return mat::getMachineEpsilon<real>(); }
201 {
return std::numeric_limits<real>::max(); }
205 {
return std::numeric_limits<real>::min(); }
211 string eigenvectors_iterative_method_,
212 real eigensolver_accuracy_,
213 int eigensolver_maxiter_,
215 bool try_eigv_on_next_iteration_if_fail_
220 string eigenvectors_iterative_method_,
221 real eigensolver_accuracy_,
222 int eigensolver_maxiter_,
224 bool try_eigv_on_next_iteration_if_fail_,
225 bool use_prev_vector_as_initial_guess_,
226 const std::vector<VectorType> &eigVecOCCRef,
227 const std::vector<VectorType> &eigVecUNOCCRef
233 std::vector<VectorType> & eigVecUNOCCref,
234 std::vector<VectorType> & eigVecOCCref,
235 std::vector<real> & eigValUNOCCref,
236 std::vector<real> & eigValOCCref
343 int& stop,
real& estim_order);
367 string eigenvectors_iterative_method_,
368 real eigensolver_accuracy_,
369 int eigensolver_maxiter_,
371 bool try_eigv_on_next_iteration_if_fail_,
372 bool use_prev_vector_as_initial_guess_
378 { nnzX =
X.nnz();
return (
double)(((double)nnzX) / ((double)
X.get_ncols() *
X.get_nrows()) * 100); }
382 {
return (
double)(((double)
X.nnz()) / ((
double)
X.get_ncols() *
X.get_nrows()) * 100); }
386 { nnzXsq =
Xsq.nnz();
return (
double)(((double)nnzXsq) / ((double)
Xsq.get_ncols() *
Xsq.get_nrows()) * 100); }
390 {
return (
double)(((double)
Xsq.nnz()) / ((
double)
Xsq.get_ncols() *
Xsq.get_nrows()) * 100); }
418 int& chosen_iter_homo);
445 void find_truncation_thresh_every_iter();
446 void find_eig_gaps_every_iter();
603template<
typename MatrixType>
610 bool use_new_stopping_criterion_,
617 assert(X.get_nrows() == X.get_ncols());
620 error_sub = error_sub_;
621 assert(error_sub>=0);
622 error_eig = error_eig_;
623 if( ! use_new_stopping_criterion_)
624 assert(error_eig>=0);
625 use_new_stopping_criterion = use_new_stopping_criterion_;
626 normPuriTrunc = norm_truncation_;
627 normPuriStopCrit = norm_stop_crit_;
629 assert(nocc >= 0 && nocc <= X.get_nrows());
631 initialized_flag =
true;
634 lumo_bounds_F = lumo_bounds_;
635 homo_bounds_F = homo_bounds_;
639#ifdef ENABLE_PRINTF_OUTPUT
644 double nnzXp = get_nnz_X(nnzX);
646 " , nocc = %d , NNZ = %lu <-> %.5lf %%",
647 X.get_nrows(), nocc, nnzX, nnzXp);
651 switch (normPuriTrunc)
666 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
670#ifdef USE_CHUNKS_AND_TASKS
679 switch (normPuriStopCrit)
694 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
698#ifdef USE_CHUNKS_AND_TASKS
706 if (use_new_stopping_criterion)
717 check_stopping_criterion_iter = -1;
718 compute_eigenvectors_in_this_SCF_cycle =
false;
723 info.stopping_criterion = (use_new_stopping_criterion ==
true);
724 info.error_subspace = error_sub;
726 info.debug_output = 0;
727#ifdef DEBUG_PURI_OUTPUT
728 info.debug_output = 1;
736template<
typename MatrixType>
738 string eigenvectors_iterative_method_,
739 real eigensolver_accuracy_,
740 int eigensolver_maxiter_,
742 bool try_eigv_on_next_iteration_if_fail_,
743 bool use_prev_vector_as_initial_guess_)
747 if(number_of_occ_eigenvectors < 0)
749 number_of_occ_eigenvectors = 1;
751 if(number_of_unocc_eigenvectors < 0)
753 number_of_unocc_eigenvectors = 1;
755 int total_number_of_eigenvectors_to_compute = number_of_occ_eigenvectors + number_of_unocc_eigenvectors;
757 if( jump_over_X_iter_proj_method < 0 )
758 jump_over_X_iter_proj_method = 1;
759 if( go_back_X_iter_proj_method < 0)
760 go_back_X_iter_proj_method = 10;
762 eigensolver_accuracy = eigensolver_accuracy_;
763 eigensolver_maxiter = eigensolver_maxiter_;
764 assert(eigensolver_maxiter >= 1);
765 scf_step = scf_step_;
766 eigenvectors_method_str = eigenvectors_method_;
767 eigenvectors_iterative_method_str = eigenvectors_iterative_method_;
768 eigenvectors_method = get_int_eig_method(eigenvectors_method_);
769 eigenvectors_iterative_method = get_int_eig_iter_method(eigenvectors_iterative_method_);
770 try_eigv_on_next_iteration_if_fail = try_eigv_on_next_iteration_if_fail_;
771 use_prev_vector_as_initial_guess = use_prev_vector_as_initial_guess_;
773 info.lumo_eigenvector_is_computed =
false;
774 info.homo_eigenvector_is_computed =
false;
780 if ((total_number_of_eigenvectors_to_compute > 0) && (eigenvectors_method ==
EIG_EMPTY))
784 "Eigenvectors will not be computed in this SCF cycle.");
788 number_of_occ_eigenvectors = 0;
789 number_of_unocc_eigenvectors = 0;
790 total_number_of_eigenvectors_to_compute = 0;
801template<
typename MatrixType>
803 string eigenvectors_iterative_method_,
804 real eigensolver_accuracy_,
805 int eigensolver_maxiter_,
807 bool try_eigv_on_next_iteration_if_fail_)
809 set_eigenvectors_params_basic(eigenvectors_method_, eigenvectors_iterative_method_, eigensolver_accuracy_, eigensolver_maxiter_, scf_step_, try_eigv_on_next_iteration_if_fail_, 0);
814template<
typename MatrixType>
816 string eigenvectors_iterative_method_,
817 real eigensolver_accuracy_,
818 int eigensolver_maxiter_,
820 bool try_eigv_on_next_iteration_if_fail_,
821 bool use_prev_vector_as_initial_guess_,
822 const std::vector<VectorType> &eigVecOCCRef,
823 const std::vector<VectorType> &eigVecUNOCCRef
826 set_eigenvectors_params_basic(eigenvectors_method_, eigenvectors_iterative_method_, eigensolver_accuracy_, eigensolver_maxiter_, scf_step_, try_eigv_on_next_iteration_if_fail_, use_prev_vector_as_initial_guess_);
829 if (use_prev_vector_as_initial_guess)
835#ifndef USE_CHUNKS_AND_TASKS
837 if (number_of_unocc_eigenvectors >= 1 && use_prev_vector_as_initial_guess)
839 assert(eigVecUNOCCRef.size() >= 1);
843 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for LUMO!");
846 eigVecLUMORef.resetSizesAndBlocks(
cols);
847 eigVecLUMORef = eigVecUNOCCRef[0];
850 if (number_of_occ_eigenvectors >= 1 && use_prev_vector_as_initial_guess)
852 assert(eigVecOCCRef.size() >= 1);
856 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for HOMO!");
859 eigVecHOMORef.resetSizesAndBlocks(
cols);
860 eigVecHOMORef = eigVecOCCRef[0];
866template<
typename MatrixType>
869 if (eigenvectors_method ==
"square")
873 if (eigenvectors_method ==
"projection")
877 if (eigenvectors_method ==
"")
881 throw std::runtime_error(
"Error in get_int_eig_method(): unknown method to compute eigenvectors");
885template<
typename MatrixType>
888 if (eigenvectors_iterative_method ==
"power")
892 if (eigenvectors_iterative_method ==
"lanczos")
896 if (eigenvectors_iterative_method ==
"")
900 throw std::runtime_error(
"Error in get_int_eig_iter_method(): unknown iterative method to compute eigenvectors");
904template<
typename MatrixType>
914 int it = iter_info.
it;
918#if defined USE_CHUNKS_AND_TASKS && defined TEST_INFTY
919 assert(XmX2_infty_norm >= 0);
923#ifndef USE_CHUNKS_AND_TASKS
926 assert(XmX2_eucl >= 0);
932 assert(XmX2_fro_norm >= 0);
933 assert(XmX2_mixed_norm >= 0);
940 assert(XmX2_fro_norm >= 0);
952template<
typename MatrixType>
958 prepare_to_purification();
961 prepare_to_purification_eigenvectors();
964#ifdef DEBUG_PURI_OUTPUT
967 purification_process();
970 eigenvalue_bounds_estimation();
976 if(compute_eigenvectors_in_this_SCF_cycle)
982 if (info.converged == 1)
984 compute_eigenvectors_without_diagonalization_last_iter_proj();
985 total_time_proj.
print(
LOG_AREA_DENSFROMF,
"Computation of eigenvalues and eigenvectors using projection method");
994 check_eigenvectors_at_the_end();
1006 info.total_time = total_time_stop;
1010template<
typename MatrixType>
1015 if (!puri_is_prepared())
1017 throw std::runtime_error(
"Error in prepare_to_purification_eigenvectors() : "
1018 "first expect a call for function prepare_to_purification().");
1022 int num_occ_eigs = get_number_of_occ_eigenvectors_to_compute();
1023 int num_unocc_eigs = get_number_of_unocc_eigenvectors_to_compute();
1024 int max_num_eigs = std::max(num_occ_eigs,num_unocc_eigs);
1030 throw std::invalid_argument(
"Error in prepare_to_purification_eigenvectors() : use projection method when requesting many eigenpairs. Set input parameter scf.eigenvectors_method = \"projection\".");
1032 if(max_num_eigs > 1 && info.method == 2 )
1034 throw std::invalid_argument(
"Error in prepare_to_purification_eigenvectors() : cannot compute more than 2 eigenpairs using SP2ACC recursive expansion. Please, use non-accelerated SP2: scf.purification_with_acceleration = 0.");
1037 if(max_num_eigs > X.get_nrows())
1038 throw std::invalid_argument(
"Error in prepare_to_purification_eigenvectors() : requested number of eigenvectors is larger than the matrix size.");
1041 do_output(
LOG_CAT_INFO,
LOG_AREA_DENSFROMF,
"(EIGV,PARAMS) go_back_X_iter_proj_method = %d, jump_over_X_iter_proj_method = %d", get_go_back_X_iter_proj_method(), get_jump_over_X_iter_proj_method());
1044 if ((num_occ_eigs > 0) || (num_unocc_eigs > 0))
1049 compute_eigenvectors_in_this_SCF_cycle =
true;
1050 info.compute_eigenvectors_in_this_SCF_cycle = compute_eigenvectors_in_this_SCF_cycle;
1051 determine_iteration_for_eigenvectors();
1064template<
typename MatrixType>
1067 if (!is_initialized())
1069 throw std::runtime_error(
"Error in prepare_to_purification() : function is called for an uninitialized class.");
1073 if (!computed_spectrum_bounds)
1076 compute_spectrum_bounds();
1079 info.time_spectrum_bounds = total_time_spectrum_bounds_stop;
1082 info.set_spectrum_bounds(spectrum_bounds.low(), spectrum_bounds.upp());
1085 map_bounds_to_0_1();
1086 set_truncation_parameters();
1095 puri_is_prepared_flag =
true;
1099template<
typename MatrixType>
1102 if (!puri_is_prepared())
1104 throw std::runtime_error(
"Error in purification_process() : "
1105 "first expect a call for function prepare_to_purification().");
1111 double Xsquare_time_stop = -1, total_time_stop = -1, trunc_time_stop = -1, purify_time_stop = -1, nnz_time_total = -1, frob_diff_time_stop = -1, eucl_diff_time_stop = -1, trace_diff_time_stop = -1, mixed_diff_time_stop = -1, stopping_criterion_time_stop = -1, inf_diff_time_stop = -1, eigenvectors_total_time = 0;
1112 int maxit_tmp = maxit;
1113 real estim_order = -1;
1114 real XmX2_trace = -1;
1115 real XmX2_fro_norm = -1;
1116 real XmX2_infty_norm = -1;
1117 real XmX2_mixed_norm = -1;
1118 real XmX2_eucl = -1;
1120 int already_stop_before = 0;
1122 info.Iterations.clear();
1123 info.Iterations.reserve(100);
1135#ifdef SAVE_MATRIX_IN_EACH_ITERATION
1139 save_matrix_now(str.str());
1145#ifdef PURI_OUTPUT_NNZ
1147 size_t nnzX_before_trunc;
1148 double nnzXp_before_trunc = get_nnz_X(nnzX_before_trunc);
1154 truncate_matrix(thresh, it);
1158#ifdef PURI_OUTPUT_NNZ
1160 if((
double)thresh >= 0)
1162 size_t nnzX_after_trunc;
1163 double nnzXp_after_trunc = get_nnz_X(nnzX_after_trunc);
1165 "Actual introduced error %e : "
1166 "nnz before %lu <-> %.2lf %% , "
1167 "nnz after %lu <-> %.2lf %%",
1168 (
double)thresh, nnzX_before_trunc, nnzXp_before_trunc, nnzX_after_trunc, nnzXp_after_trunc);
1172 size_t nnzX_after_trunc;
1173 double nnzXp_after_trunc = get_nnz_X(nnzX_after_trunc);
1175 "nnz before %lu <-> %.2lf %% , "
1176 "nnz after %lu <-> %.2lf %%",
1177 nnzX_before_trunc, nnzXp_before_trunc, nnzX_after_trunc, nnzXp_after_trunc);
1181 if((
double)thresh >= 0)
1188 Xsq = (
real)1.0 * X * X;
1192#ifdef PURI_OUTPUT_NNZ
1195 double nnzXsqp = get_nnz_Xsq(nnzXsq);
1203 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1207#if defined USE_CHUNKS_AND_TASKS && defined TEST_INFTY
1209 XmX2_infty_norm = MatrixType::infty_diff(X, Xsq);
1217#ifndef USE_CHUNKS_AND_TASKS
1219 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq,
mixed_acc);
1227 XmX2_trace = X.trace() - Xsq.trace();
1233#ifndef USE_CHUNKS_AND_TASKS
1235 XmX2_eucl = MatrixType::eucl_diff(X, Xsq,
eucl_acc);
1248#ifdef DEBUG_PURI_OUTPUT
1249 output_norms_and_traces(iter_info);
1252 if (compute_eigenvectors_in_this_SCF_cycle)
1255 compute_eigenvectors_without_diagonalization(it, iter_info);
1259 ostringstream str_out;
1260 str_out <<
"Iteration " << it;
1266 iter_info.
gap = 1 - homo_bounds.upp() - lumo_bounds.upp();
1271 iter_info.
NNZ_X = get_nnz_X();
1272 iter_info.
NNZ_X2 = get_nnz_Xsq();
1278 stopping_criterion_time_stop = 0;
1287 iter_info.
nnz_time = nnz_time_total;
1288 save_other_iter_info(iter_info, it);
1292 info.Iterations.push_back(iter_info);
1296 output_time_WriteAndReadAll();
1301 while (it <= maxit_tmp)
1316#ifdef SAVE_MATRIX_IN_EACH_ITERATION
1320 save_matrix_now(str.str());
1324#ifdef PURI_OUTPUT_NNZ
1326 size_t nnzX_before_trunc;
1327 double nnzXp_before_trunc = get_nnz_X(nnzX_before_trunc);
1333 truncate_matrix(thresh, it);
1337 #ifdef PURI_OUTPUT_NNZ
1339 if((
double)thresh >= 0)
1341 size_t nnzX_after_trunc;
1342 double nnzXp_after_trunc = get_nnz_X(nnzX_after_trunc);
1344 "Actual introduced error %e : "
1345 "nnz before %lu <-> %.2lf %% , "
1346 "nnz after %lu <-> %.2lf %%",
1347 (
double)thresh, nnzX_before_trunc, nnzXp_before_trunc, nnzX_after_trunc, nnzXp_after_trunc);
1351 size_t nnzX_after_trunc;
1352 double nnzXp_after_trunc = get_nnz_X(nnzX_after_trunc);
1354 "nnz before %lu <-> %.2lf %% , "
1355 "nnz after %lu <-> %.2lf %%",
1356 nnzX_before_trunc, nnzXp_before_trunc, nnzX_after_trunc, nnzXp_after_trunc);
1360 if((
double)thresh >= 0)
1368 Xsq = (
real)1.0 * X * X;
1371#ifdef PURI_OUTPUT_NNZ
1374 double nnzXsqp = get_nnz_Xsq(nnzXsq);
1385 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1389#if defined USE_CHUNKS_AND_TASKS && defined TEST_INFTY
1391 XmX2_infty_norm = MatrixType::infty_diff(X, Xsq);
1398#ifndef USE_CHUNKS_AND_TASKS
1400 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq,
mixed_acc);
1409#ifndef USE_CHUNKS_AND_TASKS
1411 XmX2_eucl = MatrixType::eucl_diff(X, Xsq,
eucl_acc);
1419 XmX2_trace = X.trace() - Xsq.trace();
1425 if (XmX2_trace < -1e10)
1427 throw std::runtime_error(
"Error in purification_process() : trace of X-X^2 is negative, seems as a"
1428 " misconvergence of the recursive expansion.");
1437#ifdef DEBUG_PURI_OUTPUT
1438 output_norms_and_traces(iter_info);
1454 iter_info.
nnz_time = nnz_time_total;
1459 if (compute_eigenvectors_in_this_SCF_cycle)
1462 compute_eigenvectors_without_diagonalization(it, iter_info);
1470 if (check_stopping_criterion_iter > 0 && it >= check_stopping_criterion_iter)
1473 stopping_criterion(iter_info, stop, estim_order);
1486 if ((already_stop_before == 0) && ((stop == 1) || (it == maxit)))
1499 assert(maxit_tmp <= (
int)VecPoly.size());
1500 maxit_tmp = it + additional_iterations;
1501 already_stop_before = 1;
1504 ostringstream str_out;
1505 str_out <<
"Iteration " << it;
1514 iter_info.
NNZ_X = get_nnz_X();
1515 iter_info.
NNZ_X2 = get_nnz_Xsq();
1524 if (use_new_stopping_criterion)
1526 iter_info.
order = estim_order;
1529 save_other_iter_info(iter_info, it);
1533 info.Iterations.push_back(iter_info);
1541 double nnzDp = get_nnz_X(nnzD);
1546 output_time_WriteAndReadAll();
1549 info.total_it = maxit_tmp;
1550 info.additional_iterations = additional_iterations;
1552 real acc_err_sub = total_subspace_error(maxit_tmp - additional_iterations);
1553#ifdef DEBUG_PURI_OUTPUT
1554 if (acc_err_sub != -1)
1559 info.accumulated_error_subspace = acc_err_sub;
1562 info.accumulated_time_calls_for_eigenvec_functions = eigenvectors_total_time;
1565 output_separate_total_times(info);
1570template<
typename MatrixType>
1573 #ifndef USE_CHUNKS_AND_TASKS
1581template<
typename MatrixType>
1605 if(total_Xmixed >= 0)
1608 if(total_Xfrob >= 0)
1612 if(total_Xstcrit >= 0)
1616 if(total_Xtrace >= 0)
1624template<
typename MatrixType>
1627 if (info.converged == 1)
1634 info.get_vec_mixed_norms(norms_mixed);
1636#if defined USE_CHUNKS_AND_TASKS && defined TEST_INFTY
1638 info.get_vec_infty_norms(norms_mixed);
1642 info.get_vec_frob_norms(norms_frob);
1643 info.get_vec_traces(traces);
1644 get_eigenvalue_estimates(norms_mixed, norms_frob, traces);
1657 info.homo_estim_low_F = homo_bounds_F_new.low();
1658 info.homo_estim_upp_F = homo_bounds_F_new.upp();
1659 info.lumo_estim_low_F = lumo_bounds_F_new.low();
1660 info.lumo_estim_upp_F = lumo_bounds_F_new.upp();
1667template<
typename MatrixType>
1671 info.eigValHOMO = 0;
1672 info.homo_eigenvector_is_computed =
false;
1676template<
typename MatrixType>
1679 eigVecUNOCC.clear();
1680 info.eigValLUMO = 0;
1681 info.lumo_eigenvector_is_computed =
false;
1685template<
typename MatrixType>
1688 if (info.converged != 1)
1691 discard_homo_eigenvector();
1692 discard_lumo_eigenvector();
1698 if (compute_eigenvectors_in_this_SCF_cycle)
1702 if (info.total_it < iter_for_homo)
1705 iter_for_homo, info.total_it);
1706 discard_homo_eigenvector();
1711 if ((info.total_it - iter_for_homo < ITER_DIFF) && info.homo_eigenvector_is_computed)
1714 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
1715 "thus result may not be accurate!", iter_for_homo, info.total_it);
1720 if (info.total_it < iter_for_lumo)
1723 iter_for_lumo, info.total_it);
1724 discard_lumo_eigenvector();
1729 if ((info.total_it - iter_for_lumo < ITER_DIFF) && info.lumo_eigenvector_is_computed)
1732 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, "
1733 "thus result may not be accurate!", iter_for_lumo, info.total_it);
1739 real low_lumo_F_bound = info.lumo_estim_low_F;
1740 real upp_lumo_F_bound = info.lumo_estim_upp_F;
1741 real low_homo_F_bound = info.homo_estim_low_F;
1742 real upp_homo_F_bound = info.homo_estim_upp_F;
1748 if (info.homo_eigenvector_is_computed)
1750 if ((low_homo_F_bound - flex_tolerance <= eigValHOMO) && (eigValHOMO <= upp_homo_F_bound + flex_tolerance))
1753 (
double)eigValHOMO, (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
1754 info.eigValHOMO = eigValHOMO;
1759 " discard computed HOMO eigenvector.",
1760 (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
1761 discard_homo_eigenvector();
1766 discard_homo_eigenvector();
1770 if (info.lumo_eigenvector_is_computed)
1772 if ((low_lumo_F_bound - flex_tolerance <= eigValLUMO) && (eigValLUMO <= upp_lumo_F_bound + flex_tolerance))
1775 (
double)eigValLUMO, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
1776 info.eigValLUMO = eigValLUMO;
1781 " discard computed LUMO eigenvector.",
1782 (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
1783 discard_lumo_eigenvector();
1788 discard_lumo_eigenvector();
1792 if (info.homo_eigenvector_is_computed && info.lumo_eigenvector_is_computed)
1794 real gap = eigValLUMO - eigValHOMO;
1803template<
typename MatrixType>
1806 if (spectrum_bounds.empty())
1811#ifdef DEBUG_PURI_OUTPUT
1815 real eigmin = spectrum_bounds.low();
1816 real eigmax = spectrum_bounds.upp();
1817 X.add_identity(-eigmax);
1818 X *= ((
real)1.0 / (eigmin - eigmax));
1822template<
typename MatrixType>
1825 if (spectrum_bounds.empty())
1830#ifdef DEBUG_PURI_OUTPUT
1834 real eigmin = spectrum_bounds.low();
1835 real eigmax = spectrum_bounds.upp();
1839 homo_bounds = homo_bounds_F;
1840 lumo_bounds = lumo_bounds_F;
1842 homo_bounds.intersect(spectrum_bounds);
1843 lumo_bounds.intersect(spectrum_bounds);
1845#ifdef DEBUG_PURI_OUTPUT
1846 if (homo_bounds.empty())
1850 if (lumo_bounds.empty())
1863 homo_bounds = (homo_bounds - eigmax) / (eigmin - eigmax);
1864 homo_bounds_X0 = homo_bounds;
1865 homo_bounds =
IntervalType(1 - homo_bounds.upp(), 1 - homo_bounds.low());
1869 lumo_bounds = (lumo_bounds - eigmax) / (eigmin - eigmax);
1870 lumo_bounds_X0 = lumo_bounds;
1881template<
typename MatrixType>
1884 real allowed_error = error_per_it;
1887 if (allowed_error == 0)
1892 assert((
int)VecGap.size() > it);
1904 tau = (allowed_error * VecGap[it]) / (1 + allowed_error);
1908 tau = allowed_error * 0.01;
1912#ifdef USE_CHUNKS_AND_TASKS
1913 threshold = X.thresh_frob(tau);
1915 threshold = X.thresh(tau, normPuriTrunc);
1925template<
typename MatrixType>
1928 int estim_num_iter = 0;
1930 estimate_number_of_iterations(estim_num_iter);
1931 info.estim_total_it = estim_num_iter;
1935 if (estim_num_iter < maxit)
1938 maxit = estim_num_iter;
1942 error_per_it = error_sub / estim_num_iter;
1953template<
typename MatrixType>
1957 computed_spectrum_bounds =
true;
1963template<
typename MatrixType>
1966 if (!computed_spectrum_bounds)
1968 compute_spectrum_bounds();
1970 eigmin = spectrum_bounds.low();
1971 eigmax = spectrum_bounds.upp();
1977template<
typename MatrixType>
1981 real eigminG, eigmaxG, eigmin, eigmax;
1984 X.gershgorin(eigminG, eigmaxG);
1997 eigminG -= smallNumberToExpandWith;
1998 eigmaxG += smallNumberToExpandWith;
1999#ifdef DEBUG_PURI_OUTPUT
2008#ifndef USE_CHUNKS_AND_TASKS
2015 Xshifted.add_identity((
real)(-1.0) * eigminG);
2020 eigmax = Xshifted.eucl(acc, maxIter) + eigminG + acc;
2033 Xshifted.add_identity((
real)(1.0) * eigminG);
2034 Xshifted.add_identity((
real)(-1.0) * eigmaxG);
2038 eigmin = -Xshifted.eucl(acc, maxIter) + eigmaxG - acc;
2051 if (eigmin == eigmax)
2060 computed_spectrum_bounds =
true;
2066template<
typename MatrixType>
2069 assert(check_stopping_criterion_iter > 0);
2071 int it = iter_info.
it;
2072 real XmX2_norm_it = -1, XmX2_norm_itm2 = -1;
2074 if (it < check_stopping_criterion_iter)
2079 if (use_new_stopping_criterion)
2086 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_eucl;
2092 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_fro_norm;
2098 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_mixed_norm;
2102 throw std::runtime_error(
"Error in stopping_criterion() : unknown matrix norm.");
2106 check_new_stopping_criterion(it, XmX2_norm_it, XmX2_norm_itm2, XmX2_trace, stop, estim_order);
2124 check_standard_stopping_criterion(XmX2_norm_it, stop);
2129template<
typename MatrixType>
2135 bool homoLumo_converged = (homo_bounds.upp() < error_eig &&
2136 lumo_bounds.upp() < error_eig);
2137 bool XmX2norm_converged = XmX2_norm < error_eig;
2138 if (homoLumo_converged || XmX2norm_converged)
2143#ifdef DEBUG_PURI_OUTPUT
2149template<
typename MatrixType>
2151 int& stop,
real& estim_order)
2165 return_constant_C(it, C);
2166 this->constant_C = C;
2176 (VecPoly[it - 1] != VecPoly[it])
2195#ifdef DEBUG_PURI_OUTPUT
2196 if ((VecPoly[it - 1] != VecPoly[it]) && (XmX2_norm_itm2 < 1))
2211template<
typename MatrixType>
2215 assert(it <= (
int)VecGap.size());
2219 for (
int i = 0; i <= it; ++i)
2221 if (VecGap[i] == -1)
2225 normE = info.Iterations[i].threshold_X;
2226 error += normE / (VecGap[i] - normE);
2233template<
typename MatrixType>
2236 return info.estim_total_it;
2240template<
typename MatrixType>
2243 if (info.converged == 1)
2245 return info.total_it;
2256template<
typename MatrixType>
2261 estimate_homo_lumo(XmX2_norm_mixed, XmX2_norm_frob, XmX2_trace);
2269template<
typename MatrixType>
2278 int total_it = info.total_it - info.additional_iterations;
2285 real STOP_NORM = gammaStopEstim - gammaStopEstim * gammaStopEstim;
2290 if (XmX2_norm_mixed.size() == XmX2_norm_frob.size())
2292 XmX2_norm_out = XmX2_norm_mixed;
2296 XmX2_norm_out = XmX2_norm_frob;
2299 for (
int it = total_it; it >= 0; it--)
2301 vi = XmX2_norm_frob[it];
2302 wi = XmX2_trace[it];
2303 mi = XmX2_norm_out[it];
2305 if (vi >= STOP_NORM)
2316 temp_value = vi * vi / wi;
2321 bounds_from_it[2] = bounds_from_it[0];
2322 bounds_from_it[3] = bounds_from_it[1];
2325 apply_inverse_poly_vector(it, bounds_from_it);
2328 final_bounds[0] = std::min(final_bounds[0], bounds_from_it[0]);
2329 final_bounds[1] = std::min(final_bounds[1], bounds_from_it[1]);
2331 final_bounds[2] = std::min(final_bounds[2], bounds_from_it[2]);
2332 final_bounds[3] = std::min(final_bounds[3], bounds_from_it[3]);
2336 real maxeig = spectrum_bounds.upp();
2337 real mineig = spectrum_bounds.low();
2338 lumo_bounds_F_new =
IntervalType(maxeig * (1 - final_bounds[1]) + mineig * final_bounds[1],
2339 maxeig * (1 - final_bounds[0]) + mineig * final_bounds[0]);
2340 homo_bounds_F_new =
IntervalType(mineig * (1 - final_bounds[2]) + maxeig * final_bounds[2],
2341 mineig * (1 - final_bounds[3]) + maxeig * final_bounds[3]);
2349template<
typename MatrixType>
2352 save_matrix_A_now(X, str);
2356template<
typename MatrixType>
2359#ifndef USE_CHUNKS_AND_TASKS
2360 vector<int> Itmp, I, Jtmp, J;
2361 vector<real> Vtmp, V;
2362 A.get_all_values(Itmp, Jtmp, Vtmp);
2366 for (
size_t i = 0; i < Itmp.size(); i++)
2368 nnz += (Vtmp[i] != 0);
2375 for (
size_t i = 0; i < Itmp.size(); i++)
2379 I.push_back(Itmp[i]);
2380 J.push_back(Jtmp[i]);
2381 V.push_back(Vtmp[i]);
2385 string name = str +
".mtx";
2388 throw std::runtime_error(
"Error in save_matrix_A_now : error in write_matrix_to_mtx.\n");
2400#ifdef USE_CHUNKS_AND_TASKS
2402template<
typename MatrixType>
2405 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
2409template<
typename MatrixType>
2412 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_last_iter_proj() is not implemented for CHTMatrix.");
2416template<
typename MatrixType>
2419 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization() is not implemented for CHTMatrix.");
2425template<
typename MatrixType>
2428 real start_shift = homo_bounds_F.upp();
2429 real end_shift = lumo_bounds_F.low();
2430 real dist = (end_shift - start_shift) / 15;
2431 real acc_eigv = eigensolver_accuracy;
2437 Fsq = (
real)1.0 * F * F;
2442 for (sigma = start_shift; sigma < end_shift; sigma += dist)
2445 M.add_identity(sigma * sigma);
2446 M += ((
real) - 2 * sigma) * F;
2447 M = ((
real) - 1.0) * M;
2449 vector<real> eigValTmp(1);
2450 vector<int> num_iter_solver(1, -1);
2459 eigenvectors_iterative_method_str, num_iter_solver,
2460 eigensolver_maxiter_for_F);
2468 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2470 printf(
"sigma = %lf , eigval = %lf , iters = %d\n", (
double)sigma, (
double)eigval, num_iter_solver[0]);
2476 M.add_identity(sigma * sigma);
2477 M += ((
real) - 2 * sigma) * F;
2478 M = ((
real) - 1.0) * M;
2480 vector<real> eigValTmp(1);
2481 vector<int> num_iter_solver(1, -1);
2490 eigenvectors_iterative_method_str, num_iter_solver,
2491 eigensolver_maxiter_for_F);
2498 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2500 printf(
"sigma = %lf , eigval = %lf , iters = %d\n", (
double)sigma, (
double)eigval, num_iter_solver[0]);
2505template<
typename MatrixType>
2508 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2514 bool compute_homo_now =
false;
2515 bool compute_lumo_now =
false;
2517 if (compute_eigenvectors_in_each_iteration)
2520 if (get_number_of_occ_eigenvectors_to_compute() >= 1)
2522 for (
size_t i = 0; i < good_iterations_homo.size(); ++i)
2524 if (good_iterations_homo[i] == it)
2526 compute_homo_now =
true;
2533 if (get_number_of_unocc_eigenvectors_to_compute() >= 1)
2535 for (
size_t i = 0; i < good_iterations_lumo.size(); ++i)
2537 if (good_iterations_lumo[i] == it)
2539 compute_lumo_now =
true;
2548 if (get_number_of_occ_eigenvectors_to_compute() >= 1)
2550 if (it == iter_for_homo)
2552 compute_homo_now =
true;
2555 if (get_number_of_unocc_eigenvectors_to_compute() >= 1)
2557 if (it == iter_for_lumo)
2559 compute_lumo_now =
true;
2564 if (compute_homo_now && !info.homo_eigenvector_is_computed)
2569 writeToTmpFile(Xsq);
2571 real sigma = SIGMA_HOMO_VEC[it];
2576 M.add_identity(sigma * sigma);
2577 M += ((
real) - 2 * sigma) * X;
2578 M = ((
real) - 1.0) * M;
2581 compute_eigenvector(M, eigVecOCC, eigValOCC, it,
true);
2593 readFromTmpFile(Xsq);
2596 if (compute_lumo_now && !info.lumo_eigenvector_is_computed)
2601 writeToTmpFile(Xsq);
2603 real sigma = SIGMA_LUMO_VEC[it];
2608 M.add_identity(sigma * sigma);
2609 M += ((
real) - 2 * sigma) * X;
2610 M = ((
real) - 1.0) * M;
2613 compute_eigenvector(M, eigVecUNOCC, eigValUNOCC, it,
false);
2624 readFromTmpFile(Xsq);
2634 assert(get_number_of_occ_eigenvectors_to_compute() + get_number_of_unocc_eigenvectors_to_compute() > 0);
2637 if(vec_matrices_Xi.empty())
2639 int num_allocated_matrices = info.estim_total_it;
2645 vec_matrices_Xi.resize(num_allocated_matrices+1, dummy);
2648 assert((
int)vec_matrices_Xi.size() >= it);
2651 if( it % get_jump_over_X_iter_proj_method() == 0 )
2655 vec_matrices_Xi[it] = X;
2667 if (compute_eigenvectors_in_each_iteration)
2670 info.homo_eigenvector_is_computed =
false;
2671 info.lumo_eigenvector_is_computed =
false;
2676template<
typename MatrixType>
2680 output_time_WriteAndReadAll();
2686 if (compute_eigenvectors_in_each_iteration)
2688 int total_it = info.total_it;
2693 for (
int it = 0; it <= total_it; ++it)
2695 projection_method_one_puri_iter(it);
2698 info.homo_eigenvector_is_computed =
false;
2699 info.lumo_eigenvector_is_computed =
false;
2709 int num_skip_iter_for_proj_method = std::min(get_go_back_X_iter_proj_method(), info.total_it);
2710 int step = get_jump_over_X_iter_proj_method();
2712 assert(num_skip_iter_for_proj_method >= 0);
2717 int current_iteration = std::max(info.total_it - num_skip_iter_for_proj_method, 0);
2720 while (current_iteration % step != 0) {
2721 current_iteration--;
2726 for(; current_iteration >= 0; current_iteration -= step)
2728 do_output(
LOG_CAT_INFO,
LOG_AREA_DENSFROMF,
"(EIGV,PROJ) Requested to compute %d occupied eigenvectors and %d unoccupied eigenvectors", get_number_of_occ_eigenvectors_to_compute(), get_number_of_unocc_eigenvectors_to_compute());
2730 if (! info.homo_eigenvector_is_computed)
2732 if (! info.lumo_eigenvector_is_computed)
2737 projection_method_one_puri_iter(current_iteration);
2739 if(!info.homo_eigenvector_is_computed)
2741 if(!info.lumo_eigenvector_is_computed)
2744 if(info.homo_eigenvector_is_computed && info.lumo_eigenvector_is_computed)
2749 if(!info.homo_eigenvector_is_computed || !info.lumo_eigenvector_is_computed)
2750 do_output(
LOG_CAT_INFO,
LOG_AREA_DENSFROMF,
"(EIGV,PROJ) Some eigenvectors are NOT COMPUTED, try to increase the number of allowed eigensolver iterations (input parameter scf.eigensolver_maxiter).");
2756template<
typename MatrixType>
2759 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2760 real DX_mult_time_homo_stop, DX_mult_time_lumo_stop;
2764 if (! info.homo_eigenvector_is_computed)
2769 X_homo = vec_matrices_Xi[current_iteration];
2770 readFromTmpFile(X_homo);
2782 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP, X_homo, 0, DXi);
2791 compute_eigenvector(Zh, eigVecOCC, eigValOCC, current_iteration,
true);
2795 info.Iterations[current_iteration].homo_eig_solver_time = homo_solver_time_stop;
2796 info.Iterations[current_iteration].DX_mult_homo_time = DX_mult_time_homo_stop;
2801 info.Iterations[current_iteration].orbital_homo_time = homo_total_time_stop;
2806 if (! info.lumo_eigenvector_is_computed)
2814 DXi = (
real)(-1) * DXi;
2817 compute_eigenvector(DXi, eigVecUNOCC, eigValUNOCC, current_iteration,
false);
2825 info.Iterations[current_iteration].DX_mult_lumo_time = 0;
2826 info.Iterations[current_iteration].lumo_eig_solver_time = lumo_solver_time_stop;
2827 info.Iterations[current_iteration].orbital_lumo_time = lumo_total_time_stop;
2834 if (! info.lumo_eigenvector_is_computed)
2838 X_lumo = vec_matrices_Xi[current_iteration];
2839 readFromTmpFile(X_lumo);
2849 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP, X_lumo, 0, DXi);
2855 DXi = (
real)(-1) * DXi;
2858 compute_eigenvector(DXi, eigVecUNOCC, eigValUNOCC, current_iteration,
false);
2861 info.Iterations[current_iteration].lumo_eig_solver_time = lumo_solver_time_stop;
2862 info.Iterations[current_iteration].DX_mult_lumo_time = DX_mult_time_lumo_stop;
2866 info.Iterations[current_iteration].orbital_lumo_time = lumo_total_time_stop;
2880template<
typename MatrixType>
2888 get_iterations_for_lumo_and_homo(iter_for_lumo, iter_for_homo);
2898template<
typename MatrixType>
2900 int& chosen_iter_homo)
2902 int maxiter = maxit;
2905 real homo0 = 1 - homo_bounds.upp();
2906 real lumo0 = lumo_bounds.upp();
2907 real homoi = homo0, lumoi = lumo0;
2909 real Dh_homo, Dh_lumo, Dgh_homo, Dgh_lumo,
2910 Dgh_homo_max = get_min_double(), Dgh_lumo_max = get_min_double();
2912 chosen_iter_lumo = -1;
2913 chosen_iter_homo = -1;
2915 good_iterations_homo.clear();
2916 good_iterations_lumo.clear();
2918 find_shifts_every_iter();
2920 for (
int i = 1; i <= maxiter; ++i)
2922 homoi = apply_poly(i, homoi);
2923 lumoi = apply_poly(i, lumoi);
2925 Dh_homo = compute_derivative(i, homo0, dummy);
2926 Dh_lumo = compute_derivative(i, lumo0, dummy);
2929 Dgh_homo = 2 * (homoi - SIGMA_HOMO_VEC[i]) * Dh_homo;
2930 Dgh_lumo = 2 * (lumoi - SIGMA_LUMO_VEC[i]) * Dh_lumo;
2932 if (homoi >= SIGMA_HOMO_VEC[i])
2934 good_iterations_homo.push_back(i);
2935 if (Dgh_homo >= Dgh_homo_max)
2937 Dgh_homo_max = Dgh_homo;
2938 chosen_iter_homo = i;
2942 if (lumoi <= SIGMA_LUMO_VEC[i])
2944 good_iterations_lumo.push_back(i);
2948 chosen_iter_lumo = i;
2953 if ((chosen_iter_homo == -1) || (chosen_iter_lumo == -1))
2955 throw "Error in get_iterations_for_lumo_and_homo() : something went wrong, cannot choose iteration to compute HOMO or LUMO eigenvector.";
2961template<
typename MatrixType>
2964 int maxiter = maxit;
2966 SIGMA_HOMO_VEC.resize(maxiter + 1);
2967 SIGMA_LUMO_VEC.resize(maxiter + 1);
2971 real homo = 1 - homo_bounds.upp();
2972 real lumo = lumo_bounds.upp();
2973 real homo_out = 1 - homo_bounds.low();
2974 real lumo_out = lumo_bounds.low();
2976 for (
int i = 1; i <= maxiter; ++i)
2978 homo = apply_poly(i, homo);
2979 lumo = apply_poly(i, lumo);
2981 homo_out = apply_poly(i, homo_out);
2982 lumo_out = apply_poly(i, lumo_out);
2984 SIGMA_HOMO_VEC[i] = (homo_out + lumo) / 2;
2985 SIGMA_LUMO_VEC[i] = (lumo_out + homo) / 2;
2991template<
typename MatrixType>
2994 int maxiter = this->maxit;
2996 this->ITER_ERROR_VEC.
clear();
2997 this->ITER_ERROR_VEC.resize(maxiter + 1);
2999 real error = error_per_it;
3000 if (error_per_it == 0)
3002 error = this->get_epsilon();
3005 for (
int i = 1; i <= maxiter; ++i)
3007 this->ITER_ERROR_VEC[i] = (error * this->VecGap[i]) / (1 + error);
3012template<
typename MatrixType>
3015 int maxiter = maxit;
3017 EIG_REL_GAP_HOMO_VEC.resize(maxiter + 1);
3018 EIG_REL_GAP_LUMO_VEC.
clear();
3019 EIG_REL_GAP_LUMO_VEC.resize(maxiter + 1);
3020 EIG_ABS_GAP_HOMO_VEC.clear();
3021 EIG_ABS_GAP_HOMO_VEC.resize(maxiter + 1);
3022 EIG_ABS_GAP_LUMO_VEC.clear();
3023 EIG_ABS_GAP_LUMO_VEC.resize(maxiter + 1);
3027 real homo0 = homo_bounds_X0.low();
3028 real lumo0 = lumo_bounds_X0.upp();
3032 real homo_map, lumo_map, one_map, zero_map, sigma;
3034 real homo = homo0, lumo = lumo0;
3036 for (
int i = 1; i <= maxiter; ++i)
3038 homo = apply_poly(i, homo);
3039 lumo = apply_poly(i, lumo);
3041 sigma = SIGMA_HOMO_VEC[i];
3043 homo_map = (homo - sigma) * (homo - sigma);
3044 lumo_map = (lumo - sigma) * (lumo - sigma);
3045 one_map = (one - sigma) * (one - sigma);
3046 zero_map = (
zero - sigma) * (zero - sigma);
3048 EIG_ABS_GAP_HOMO_VEC[i] =
min(lumo_map - homo_map, one_map - homo_map);
3049 EIG_REL_GAP_HOMO_VEC[i] = EIG_ABS_GAP_HOMO_VEC[i] /
3050 max(zero_map - homo_map, one_map - homo_map);
3053 homo = homo0, lumo = lumo0;
3055 for (
int i = 1; i <= maxiter; ++i)
3057 homo = apply_poly(i, homo);
3058 lumo = apply_poly(i, lumo);
3060 sigma = SIGMA_LUMO_VEC[i];
3062 homo_map = (homo - sigma) * (homo - sigma);
3063 lumo_map = (lumo - sigma) * (lumo - sigma);
3064 zero_map = (
zero - sigma) * (zero - sigma);
3065 one_map = (one - sigma) * (one - sigma);
3067 EIG_ABS_GAP_LUMO_VEC[i] =
min(homo_map - lumo_map, zero_map - lumo_map);
3068 EIG_REL_GAP_LUMO_VEC[i] = EIG_ABS_GAP_LUMO_VEC[i] /
3069 max(zero_map - lumo_map, one_map - lumo_map);
3077#ifdef USE_CHUNKS_AND_TASKS
3078template<
typename MatrixType>
3081 throw "compute_eigenvector() is not implemented with Chunks and Tasks.";
3086template<
typename MatrixType>
3089 real acc_eigv = eigensolver_accuracy;
3091 #ifdef DEBUG_PURI_OUTPUT
3095#ifdef SAVE_MATRIX_IN_EACH_ITERATION
3098 str <<
"M_" << is_homo <<
"_" << it;
3099 save_matrix_A_now(M, str.str());
3103 int number_of_eigenvalues_to_compute;
3105 number_of_eigenvalues_to_compute = get_number_of_occ_eigenvectors_to_compute();
3107 number_of_eigenvalues_to_compute = get_number_of_unocc_eigenvectors_to_compute();
3108 assert(number_of_eigenvalues_to_compute >= 0);
3111 if(number_of_eigenvalues_to_compute == 0)
return;
3124 eigVec = vector<VectorType>(number_of_eigenvalues_to_compute,
VectorType(
rows));
3126 vector<real> eigVal_of_M(number_of_eigenvalues_to_compute);
3127 eigVal.resize(number_of_eigenvalues_to_compute);
3128 vector<int> num_iter_solver(number_of_eigenvalues_to_compute, -1);
3130 if (use_prev_vector_as_initial_guess)
3136 eigVec[0]= eigVecHOMORef;
3141 eigVec[0] = eigVecLUMORef;
3153 bool eigensolver_converged =
false;
3157 number_of_eigenvalues_to_compute,
3158 eigenvectors_iterative_method_str, num_iter_solver,
3159 eigensolver_maxiter);
3161 eigensolver_converged =
true;
3166 eigensolver_converged =
false;
3169 eigVal_of_M.clear();
3172 catch (std::exception &e)
3175 eigensolver_converged =
false;
3178 eigVal_of_M.clear();
3189#ifdef SAVE_EIGEVECTORS_TO_FILES
3216 for (
size_t i = 0; i < number_of_eigenvalues_to_compute; i++) {
3217 save_selected_eigenvector_to_file(eigVec[i], i, is_homo, it);
3222 if (eigensolver_converged)
3225 for (
int i = 0; i < number_of_eigenvalues_to_compute; i++) {
3226 real tau = std::max( error_sub, get_epsilon() );
3227 if (eigVal_of_M[i] < tau) {
3229 " the computed eigenvector may not be accurate (tol = %e)", i, (
double)eigVal_of_M[i], it, tau);
3233 bool is_homo_computed_now = is_homo, is_lumo_computed_now = !is_homo;
3236 real eigValHOMO_or_LUMO;
3237 check_homo_lumo_eigenvalues(eigValHOMO_or_LUMO, eigVec[0], is_homo_computed_now, is_lumo_computed_now, it);
3238 if (is_homo_computed_now)
3240 really_good_iterations_homo.push_back(it);
3241 info.homo_eigenvector_is_computed =
true;
3242 info.homo_eigenvector_is_computed_in_iter = it;
3243 info.homo_eigensolver_iter = num_iter_solver[0];
3244 info.homo_eigensolver_time = eigv_elapsed_wall_sec;
3245 eigValHOMO = eigValHOMO_or_LUMO;
3246 eigVal[0] = eigValHOMO_or_LUMO;
3248 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
3250 else if (is_lumo_computed_now)
3252 really_good_iterations_lumo.push_back(it);
3253 info.lumo_eigenvector_is_computed =
true;
3254 info.lumo_eigenvector_is_computed_in_iter = it;
3255 info.lumo_eigensolver_iter = num_iter_solver[0];
3256 info.lumo_eigensolver_time = eigv_elapsed_wall_sec;
3257 eigValLUMO = eigValHOMO_or_LUMO;
3258 eigVal[0] = eigValHOMO_or_LUMO;
3260 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
3265 ": number of eigensolver iterations is %d", it, num_iter_solver[0]);
3269 eigVal_of_M.clear();
3274 if(is_lumo_computed_now || is_homo_computed_now)
3276 assert((
int)eigVal.size() >= number_of_eigenvalues_to_compute);
3278 for (
int i = 1; i < number_of_eigenvalues_to_compute; i++) {
3279 eigVal[i] = eigvec::compute_rayleigh_quotient<real>(F, eigVec[i]);
3292template<
typename MatrixType>
3296#ifndef USE_CHUNKS_AND_TASKS
3302template<
typename MatrixType>
3306#ifndef USE_CHUNKS_AND_TASKS
3312template<
typename MatrixType>
3316 std::vector<real> fullVector;
3320 ostringstream filename;
3322 filename <<
"occ_eigv_" << num;
3324 filename <<
"unocc_eigv_" << num;
3328 filename <<
"_it_" << it;
3333 filename <<
"_scf_step_" << scf_step;
3338 if (
write_vector(filename.str().c_str(), fullVector) == -1)
3340 throw std::runtime_error(
"Error in save_selected_eigenvector_to_file() : error in write_vector.");
3345template<
typename MatrixType>
3349 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3351 eigVal = approx_eig;
3355template<
typename MatrixType>
3358 assert(is_homo || is_lumo);
3360 bool is_homo_init = is_homo;
3361 bool is_lumo_init = is_lumo;
3373 real low_lumo_F_bound, low_homo_F_bound;
3374 real upp_lumo_F_bound, upp_homo_F_bound;
3379 low_lumo_F_bound = lumo_bounds_F.low();
3380 upp_lumo_F_bound = lumo_bounds_F.upp();
3381 low_homo_F_bound = homo_bounds_F.low();
3382 upp_homo_F_bound = homo_bounds_F.upp();
3387 low_lumo_F_bound = info.lumo_estim_low_F;
3388 upp_lumo_F_bound = info.lumo_estim_upp_F;
3389 low_homo_F_bound = info.homo_estim_low_F;
3390 upp_homo_F_bound = info.homo_estim_upp_F;
3394 throw std::runtime_error(
"Error in check_homo_lumo_eigenvalues() : unexpected eigenvectors_method value.");
3397#ifdef DEBUG_PURI_OUTPUT
3402 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3405 eigVal = approx_eig;
3410 if ((approx_eig <= upp_homo_F_bound + flex_tolerance) && (approx_eig >= low_homo_F_bound - flex_tolerance))
3416 "HOMO bounds are [ %lf , %lf ]", (
double)approx_eig, (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
3418 iter_for_homo = iter;
3420 if (is_lumo_init && (eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3422 iter_for_lumo = iter_for_lumo + 1;
3427 else if ((approx_eig <= upp_lumo_F_bound + flex_tolerance) && (approx_eig >= low_lumo_F_bound - flex_tolerance))
3433 "LUMO interval [ %lf , %lf ]", (
double)approx_eig, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
3435 iter_for_lumo = iter;
3437 if (is_homo_init && (eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3439 iter_for_homo = iter_for_homo + 1;
3450 (
double)approx_eig, (
double)low_homo_F_bound, (
double)upp_homo_F_bound, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
3455 if ((eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3457 iter_for_homo = iter_for_homo + 1;
3469template<
typename MatrixType>
3473 f.open(filename, std::ios::out);
3480 int it = info.total_it;
3483 for (
int i = 0; i <= it; ++i)
3485 f << VecPoly[i] <<
" ";
3487 f <<
"];" << std::endl;
3493 for (
int i = 0; i <= it; ++i)
3495 f << (double)info.Iterations[i].XmX2_eucl <<
" ";
3497 f <<
"];" << std::endl;
3498 f <<
" norm_letter = '2';" << std::endl;
3503 for (
int i = 0; i <= it; ++i)
3505 f << (double)info.Iterations[i].XmX2_fro_norm <<
" ";
3507 f <<
"];" << std::endl;
3508 f <<
" norm_letter = 'F';" << std::endl;
3513 for (
int i = 0; i <= it; ++i)
3515 f << (double)info.Iterations[i].XmX2_mixed_norm <<
" ";
3517 f <<
"];" << std::endl;
3518 f <<
" norm_letter = 'M';" << std::endl;
3522 throw "Wrong norm in PurificationGeneral::gen_matlab_file_norm_diff()";
3525 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3526 f <<
"it = " << it <<
";" << std::endl;
3527 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3528 f <<
"fighandle = figure; clf;" << std::endl;
3529 f <<
"MARKER = ['o', '+'];" << std::endl;
3530 f <<
"semilogy(0:stop_iteration, X_norm(1:stop_iteration+1), '-b', plot_props{:});" << std::endl;
3531 f <<
"hold on" << std::endl;
3532 f <<
"for i = 1:stop_iteration+1" << std::endl;
3533 f <<
" if POLY(i) == 1" << std::endl;
3534 f <<
" h1 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3535 f <<
" else" << std::endl;
3536 f <<
" h2 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3537 f <<
" end" << std::endl;
3538 f <<
"end" << std::endl;
3539 f <<
"if stop_iteration ~= it" << std::endl;
3540 f <<
"h3 = semilogy(stop_iteration+1:it, X_norm(stop_iteration+2:it+1), '-.vr', plot_props{:});" << std::endl;
3541 f <<
"semilogy(stop_iteration:stop_iteration+1, X_norm(stop_iteration+1:stop_iteration+2), '-.r', plot_props{:});" << std::endl;
3542 f <<
"legend([h1 h2 h3],{'$x^2$', '$2x-x^2$', 'After stop'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3543 f <<
"else" << std::endl;
3544 f <<
"legend([h1 h2],{'$x^2$', '$2x-x^2$'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3545 f <<
"end" << std::endl;
3546 f <<
"xlabel('Iteration SP2', 'Interpreter','latex');" << std::endl;
3547 f <<
"ylabel({['$\\|X_i-X_i^2\\|_{' norm_letter '}$']},'interpreter','latex');" << std::endl;
3548 f <<
"grid on" << std::endl;
3549 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3550 f <<
"set(gca,'FontSize',20);" << std::endl;
3551 f <<
"xlim([0 it]);" << std::endl;
3552 f <<
"ylim([-inf inf]);" << std::endl;
3553 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3554 f <<
"a = 16; S = logspace(-a, 1, a+2);" << std::endl;
3555 f <<
"set(gca,'YTick',S(1:2:end));" << std::endl;
3557 f <<
"hold off" << std::endl;
3559 f <<
"% print( fighandle, '-depsc2', 'norm_diff.eps');" << std::endl;
3565template<
typename MatrixType>
3569 f.open(filename, std::ios::out);
3576 int it = info.total_it;
3579 for (
int i = 0; i <= it; ++i)
3581 f << (double)info.Iterations[i].threshold_X <<
" ";
3583 f <<
"];" << std::endl;
3585 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3586 f <<
"it = " << it <<
";" << std::endl;
3587 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3588 f <<
"fighandle = figure; clf;" << std::endl;
3589 f <<
"semilogy(0:stop_iteration, Thresh(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3590 f <<
"hold on" << std::endl;
3591 f <<
"if stop_iteration ~= it" << std::endl;
3592 f <<
"semilogy(stop_iteration+1:it, Thresh(stop_iteration+2:it+1), '-^r', plot_props{:});" << std::endl;
3593 f <<
"semilogy(stop_iteration:stop_iteration+1, Thresh(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3594 f <<
"legend('before stop', 'after stop', 'Location','NorthWest');" << std::endl;
3595 f <<
"end" << std::endl;
3596 f <<
"grid on" << std::endl;
3597 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3598 f <<
"set(gca,'FontSize',20);" << std::endl;
3599 f <<
"xlim([0 it]);" << std::endl;
3600 f <<
"ylim([-inf inf]);" << std::endl;
3601 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3602 f <<
"hold off" << std::endl;
3603 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3604 f <<
"ylabel('Threshold value', 'interpreter','latex');" << std::endl;
3605 f <<
"% print( fighandle, '-depsc2', 'threshold.eps');" << std::endl;
3611template<
typename MatrixType>
3615 f.open(filename, std::ios::out);
3622 int it = info.total_it;
3625 for (
int i = 0; i <= it; ++i)
3627 f << (double)info.Iterations[i].NNZ_X <<
" ";
3629 f <<
"];" << std::endl;
3633 for (
int i = 0; i <= it; ++i)
3635 f << (double)info.Iterations[i].NNZ_X2 <<
" ";
3637 f <<
"];" << std::endl;
3641 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3642 f <<
"it = " << it <<
";" << std::endl;
3643 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3644 f <<
"fighandle = figure; clf;" << std::endl;
3645 f <<
"h2 = plot(0:stop_iteration, NNZ_X(1:stop_iteration+1), '-ob', plot_props{:});" << std::endl;
3646 f <<
"hold on" << std::endl;
3647 f <<
"h1 = plot(0:stop_iteration, NNZ_X2(1:stop_iteration+1), '-vm', plot_props{:});" << std::endl;
3648 f <<
"if stop_iteration ~= it" << std::endl;
3649 f <<
"plot(stop_iteration+1:it, NNZ_X(stop_iteration+2:it+1), '-vr', plot_props{:});" << std::endl;
3650 f <<
"plot(stop_iteration+1:it, NNZ_X2(stop_iteration+2:it+1), '-*r', plot_props{:});" << std::endl;
3651 f <<
"plot(stop_iteration:stop_iteration+1, NNZ_X(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3652 f <<
"plot(stop_iteration:stop_iteration+1, NNZ_X2(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3653 f <<
"end" << std::endl;
3654 f <<
"legend([h1, h2], {'$X^2$', '$X$'}, 'interpreter','latex');" << std::endl;
3655 f <<
"grid on" << std::endl;
3656 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3657 f <<
"set(gca,'FontSize',20);" << std::endl;
3658 f <<
"xlim([0 it]);" << std::endl;
3659 f <<
"ylim([0 inf]);" << std::endl;
3660 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3661 f <<
"hold off" << std::endl;
3662 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3663 f <<
"ylabel('NNZ [\\%]', 'interpreter','latex');" << std::endl;
3664 f <<
"% print( fighandle, '-depsc2', 'nnz.eps');" << std::endl;
3709template<
typename MatrixType>
3713 f.open(filename, std::ios::out);
3720 int it = info.total_it;
3723 for (
int i = 0; i <= it; ++i)
3725 f << (double)info.Iterations[i].homo_bound_low <<
" ";
3727 f <<
"];" << std::endl;
3731 for (
int i = 0; i <= it; ++i)
3733 f << (double)info.Iterations[i].homo_bound_upp <<
" ";
3735 f <<
"];" << std::endl;
3739 for (
int i = 0; i <= it; ++i)
3741 f << (double)info.Iterations[i].lumo_bound_low <<
" ";
3743 f <<
"];" << std::endl;
3747 for (
int i = 0; i <= it; ++i)
3749 f << (double)info.Iterations[i].lumo_bound_upp <<
" ";
3751 f <<
"];" << std::endl;
3755 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3756 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3757 f <<
"fighandle = figure; clf;" << std::endl;
3758 f <<
"semilogy(0:stop_iteration, homo_upp(1:stop_iteration+1), '-^b', plot_props{:});" << std::endl;
3759 f <<
"hold on" << std::endl;
3760 f <<
"semilogy(0:stop_iteration, homo_low(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3761 f <<
"semilogy(0:stop_iteration, lumo_low(1:stop_iteration+1), '-vr', plot_props{:});" << std::endl;
3762 f <<
"semilogy(0:stop_iteration, lumo_upp(1:stop_iteration+1), '-^r', plot_props{:});" << std::endl;
3763 f <<
"grid on" << std::endl;
3764 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3765 f <<
"set(gca,'FontSize',20);" << std::endl;
3766 f <<
"xlim([0 stop_iteration]);" << std::endl;
3767 f <<
"set(gca,'XTick',[0 5:5:stop_iteration]);" << std::endl;
3768 f <<
"ylim([-inf inf]);" << std::endl;
3769 f <<
"hold off" << std::endl;
3770 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3771 f <<
"ylabel('Homo and lumo bounds', 'interpreter','latex');" << std::endl;
3772 f <<
"% print( fighandle, '-depsc2', 'eigs.eps');" << std::endl;
3778template<
typename MatrixType>
3782 f.open(filename, std::ios::out);
3789 int it = info.total_it;
3791 f <<
"time_total = [";
3792 for (
int i = 0; i <= it; ++i)
3794 f << std::setprecision(16) << (double)info.Iterations[i].total_time <<
" ";
3796 f <<
"];" << std::endl;
3798 f <<
"time_square = [";
3799 for (
int i = 0; i <= it; ++i)
3801 f << std::setprecision(16) << (double)info.Iterations[i].Xsquare_time <<
" ";
3803 f <<
"];" << std::endl;
3805 f <<
"time_trunc = [";
3806 for (
int i = 0; i <= it; ++i)
3808 f << std::setprecision(16) << (double)info.Iterations[i].trunc_time <<
" ";
3810 f <<
"];" << std::endl;
3812 if (info.compute_eigenvectors_in_this_SCF_cycle)
3814 f <<
"time_eigenvectors_homo = [";
3815 for (
int i = 0; i <= it; ++i)
3817 f << std::setprecision(16) << (double)info.Iterations[i].orbital_homo_time <<
" ";
3819 f <<
"];" << std::endl;
3820 f <<
"time_eigenvectors_lumo = [";
3821 for (
int i = 0; i <= it; ++i)
3823 f << std::setprecision(16) << (double)info.Iterations[i].orbital_lumo_time <<
" ";
3825 f <<
"];" << std::endl;
3827 f <<
"time_solver_homo = [";
3828 for (
int i = 0; i <= it; ++i)
3830 f << std::setprecision(16) << (double)info.Iterations[i].homo_eig_solver_time <<
" ";
3832 f <<
"];" << std::endl;
3833 f <<
"time_solver_lumo = [";
3834 for (
int i = 0; i <= it; ++i)
3836 f << std::setprecision(16) << (double)info.Iterations[i].lumo_eig_solver_time <<
" ";
3838 f <<
"];" << std::endl;
3844 f <<
"X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_eigenvectors_homo-time_solver_homo;"
3845 " time_eigenvectors_lumo-time_solver_lumo; "
3846 " time_total - time_square - time_trunc - time_eigenvectors_homo - time_eigenvectors_lumo];" << std::endl;
3850 f <<
"time_DX_homo = [";
3851 for (
int i = 0; i <= it; ++i)
3853 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_homo_time <<
" ";
3855 f <<
"];" << std::endl;
3856 f <<
"time_DX_lumo = [";
3857 for (
int i = 0; i <= it; ++i)
3859 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_lumo_time <<
" ";
3861 f <<
"];" << std::endl;
3865 f <<
"X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_DX_homo; time_DX_lumo;"
3866 " time_eigenvectors_homo - time_DX_homo - time_solver_homo; time_eigenvectors_lumo - time_DX_lumo - time_solver_lumo;"
3867 " time_total - time_square - time_trunc];" << std::endl;
3872 f <<
"X = [time_square; time_trunc; time_total - time_square - time_trunc];" << std::endl;
3875 f <<
"it = " << it <<
";" << std::endl;
3876 f <<
"xtick = 0:it;" << std::endl;
3877 f <<
"fighandle = figure; clf;" << std::endl;
3878 f <<
"b=bar(xtick, X', 'stacked');" << std::endl;
3879 f <<
"axis tight;" << std::endl;
3880 f <<
"grid on" << std::endl;
3881 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3882 f <<
"set(gca,'FontSize',20);" << std::endl;
3883 f <<
"xlim([0 it]);" << std::endl;
3884 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3885 f <<
"ylim([-inf inf]);" << std::endl;
3886 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3887 f <<
"ylabel('Time [s]', 'interpreter','latex');" << std::endl;
3888 if (info.compute_eigenvectors_in_this_SCF_cycle)
3895 f <<
"legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3899 f <<
"legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'DX (lumo)', 'DX (homo)', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3904 f <<
"legend(flipud(b(:)), {'other', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3906 f <<
"% print( fighandle, '-depsc2', 'time.eps');" << std::endl;
generalVector VectorType
Definition GetDensFromFock.cc:62
Definition puri_info.h:52
real homo_bound_upp
Definition puri_info.h:88
real threshold_X
Definition puri_info.h:57
real lumo_bound_upp
Definition puri_info.h:90
double mixed_diff_time
Definition puri_info.h:65
double nnz_time
Definition puri_info.h:67
real XmX2_eucl
Definition puri_info.h:79
double trunc_time
Definition puri_info.h:59
double Xsquare_time
Definition puri_info.h:58
double lumo_eig_solver_time
Definition puri_info.h:74
real gap
Definition puri_info.h:82
real order
Definition puri_info.h:80
real NNZ_X2
Definition puri_info.h:84
real homo_bound_low
Definition puri_info.h:87
real constantC
Definition puri_info.h:96
int it
Definition puri_info.h:56
double frob_diff_time
Definition puri_info.h:66
real lumo_bound_low
Definition puri_info.h:89
double trace_diff_time
Definition puri_info.h:64
double purify_time
Definition puri_info.h:60
real XmX2_infty_norm
Definition puri_info.h:77
double inf_diff_time
Definition puri_info.h:68
real XmX2_fro_norm
Definition puri_info.h:76
double orbital_homo_time
Definition puri_info.h:69
double total_time
Definition puri_info.h:61
double homo_eig_solver_time
Definition puri_info.h:73
double orbital_lumo_time
Definition puri_info.h:70
real XmX2_mixed_norm
Definition puri_info.h:78
real NNZ_X
Definition puri_info.h:83
real XmX2_trace
Definition puri_info.h:75
double eucl_diff_time
Definition puri_info.h:63
double stopping_criterion_time
Definition puri_info.h:62
Definition puri_info.h:141
real get_total_Xtrunc_time()
Definition puri_info.cc:56
real get_total_mixed_diff_time()
Definition puri_info.cc:87
real get_total_nnz_time()
Definition puri_info.cc:64
real get_total_stopping_criterion_time()
Definition puri_info.cc:102
real get_total_purify_time()
Definition puri_info.cc:118
real get_total_trace_diff_time()
Definition puri_info.cc:110
real get_total_Xsquare_time()
Definition puri_info.cc:48
real get_total_frob_diff_time()
Definition puri_info.cc:94
real get_total_eucl_diff_time()
Definition puri_info.cc:80
real get_total_inf_diff_time()
Definition puri_info.cc:73
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition purification_general.h:117
void writeToTmpFile(MatrixType &A) const
Definition purification_general.h:3294
bool computed_spectrum_bounds
Definition purification_general.h:545
void unset_compute_eigenvectors_in_each_iteration()
Definition purification_general.h:246
int get_number_of_unocc_eigenvectors_to_compute() const
Definition purification_general.h:254
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition purification_general.h:482
void set_compute_eigenvectors_in_each_iteration()
Definition purification_general.h:245
void get_spectrum_bounds(real &eigmin, real &eigmax)
Get spectrum bounds.
Definition purification_general.h:1964
virtual void stopping_criterion(IterationInfo &iter_info, int &stop, real &estim_order)
Choose stopping criterion and check it.
Definition purification_general.h:2067
void set_go_back_X_iter_proj_method(int val)
Definition purification_general.h:261
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition purification_general.h:513
void compute_eigenvectors_without_diagonalization(int it, IterationInfo &iter_info)
Compute HOMO and LUMO eigenvalues and eigenvectors of the matrix F.
Definition purification_general.h:2506
std::vector< real > VectorTypeReal
Definition purification_general.h:123
std::vector< real > eigValUNOCC
Here we save eigenvalues corresponding to the unoccupied orbitals.
Definition purification_general.h:571
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)=0
VectorTypeReal EIG_REL_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for homo eigenvalue.
Definition purification_general.h:519
MatrixType F
Matrix F.
Definition purification_general.h:291
int get_est_number_of_puri_iterations()
Definition purification_general.h:2234
virtual real total_subspace_error(int it)
Definition purification_general.h:2213
int number_of_occ_eigenvectors
Definition purification_general.h:524
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition purification_general.h:503
virtual ~PurificationGeneral()
Create MATLAB .m file which plots a condition number of a problem of computing the density matrix in ...
Definition purification_general.h:287
virtual void PurificationStart()
Start recursive expansion.
Definition purification_general.h:953
virtual void truncate_matrix(real &thresh, int it)
Definition purification_general.h:1882
int get_go_back_X_iter_proj_method() const
Definition purification_general.h:263
VectorTypeReal EIG_ABS_GAP_LUMO_VEC
Definition purification_general.h:518
virtual real compute_derivative(const int it, real x, real &DDf)=0
virtual void check_new_stopping_criterion(const int it, const real XmX2_norm_it, const real XmX2_norm_itm2, const real XmX2_trace, int &stop, real &estim_order)
Check stopping criterion.
Definition purification_general.h:2150
static real get_min_double()
Get smallest number.
Definition purification_general.h:204
VectorTypeReal ITER_ERROR_VEC
(Eigenvectors) Maximum error introduced in each iteration.
Definition purification_general.h:516
std::vector< int > VectorTypeInt
Definition purification_general.h:122
int maxit
Maximum number of iterations.
Definition purification_general.h:481
VectorTypeInt really_good_iterations_lumo
Iterations where lumo eigenvector is actually computed.
Definition purification_general.h:584
int get_int_eig_method(string eigenvectors_method)
Definition purification_general.h:867
NormType normPuriStopCrit
Norm used in the stopping criterion Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm.
Definition purification_general.h:492
IntervalType lumo_bounds_X0
Initial lumo bounds for X.
Definition purification_general.h:535
IntervalType homo_bounds
(1-homo) bounds for Xi in iteration i
Definition purification_general.h:531
virtual void check_eigenvectors_at_the_end()
Definition purification_general.h:1686
string eigenvectors_method_str
Definition purification_general.h:557
void estimate_homo_lumo(const VectorTypeReal &XmX2_norm_mixed, const VectorTypeReal &XmX2_norm_frob, const VectorTypeReal &XmX2_trace)
Get homo and lumo bounds from traces and norms of Xi-Xi^2.
Definition purification_general.h:2270
IntervalType lumo_bounds_F
Initial homo bounds for F.
Definition purification_general.h:538
virtual void initialize(const MatrixType &F_, const IntervalType &lumo_bounds_, const IntervalType &homo_bounds_, int maxit_, real error_sub_, real error_eig_, bool use_new_stopping_criterion_, NormType norm_truncation, NormType norm_stop_crit, int nocc_)
Set imporatant parameters for the recursive expansion.
Definition purification_general.h:604
virtual void purify_bounds(const int it)=0
MatrixType MatrixTypeWrapper
Definition purification_general.h:125
VectorTypeInt really_good_iterations_homo
Iterations where homo eigenvector is actually computed.
Definition purification_general.h:583
int go_back_X_iter_proj_method
Definition purification_general.h:528
double get_nnz_Xsq()
Get nnz of X^2 in %.
Definition purification_general.h:389
void gen_matlab_file_time(const char *filename) const
Create MATLAB .m file which creates a bar plot presenting time spent on various parts of the iteratio...
Definition purification_general.h:3779
virtual void return_constant_C(const int it, real &Cval)=0
virtual void discard_lumo_eigenvector()
Definition purification_general.h:1677
real error_eig
Error in eigenvalues (used just in old stopping criterion).
Definition purification_general.h:497
void readFromTmpFile(MatrixType &A) const
Definition purification_general.h:3304
IntervalType lumo_bounds_F_new
Definition purification_general.h:541
static real get_epsilon()
Get machine epsilon.
Definition purification_general.h:196
void set_number_of_eigenvectors_to_compute(const int occ, const int unocc)
Definition purification_general.h:248
virtual void set_init_params()=0
MatrixType Xsq
Matrix X^2.
Definition purification_general.h:132
bool compute_eigenvectors_in_this_SCF_cycle
Definition purification_general.h:562
bool initialized_flag
Definition purification_general.h:475
virtual void purification_process()
Run recursive expansion.
Definition purification_general.h:1100
generalVector VectorType
Definition purification_general.h:124
void save_selected_eigenvector_to_file(const VectorType &v, int num, bool is_homo, int it=-1)
Definition purification_general.h:3314
void find_shifts_every_iter()
/brief Find shifts sigma which will be used for construction of the filtering polynomial for computin...
Definition purification_general.h:2962
ergo_real real
Definition purification_general.h:119
virtual void eigenvalue_bounds_estimation()
Estimate eigenvalues near homo-lumo gap.
Definition purification_general.h:1625
void set_eigenvectors_params(string eigenvectors_method_, string eigenvectors_iterative_method_, real eigensolver_accuracy_, int eigensolver_maxiter_, int scf_step_, bool try_eigv_on_next_iteration_if_fail_)
Set parameters for computing eigenvectors.
Definition purification_general.h:802
virtual void set_truncation_parameters()
Definition purification_general.h:1926
IntervalType spectrum_bounds
Outer bounds for the whole spectrum of F/Xi.
Definition purification_general.h:544
VectorTypeReal SIGMA_LUMO_VEC
(Eigenvectors) Approximation of shifts in each iteration.
Definition purification_general.h:517
real eigValHOMO
Definition purification_general.h:575
static real get_max_double()
Get largest number.
Definition purification_general.h:200
real eigensolver_accuracy
Accuracy of the eigenvalue problem solver.
Definition purification_general.h:553
VectorTypeReal EIG_ABS_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for lumo eigenvalue.
Definition purification_general.h:518
PurificationGeneral()
Definition purification_general.h:135
void gen_matlab_file_nnz(const char *filename) const
Create MATLAB .m file which plots the number of non-zero elements in matrices X_i and X_i^2 in each r...
Definition purification_general.h:3612
real eigValLUMO
Definition purification_general.h:574
real constant_C
Asymptotic constant C needed for the new stopping criterion.
Definition purification_general.h:501
int get_number_of_occ_eigenvectors_to_compute() const
Definition purification_general.h:252
NormType normPuriTrunc
Norm used for the truncation of matrices.
Definition purification_general.h:488
void gen_matlab_file_eigs(const char *filename) const
Create MATLAB .m file which plots the homo and lumo bounds in each recursive expansion iteration.
Definition purification_general.h:3710
void get_eigenvalue_estimates(const VectorTypeReal &XmX2_norm_mixed, const VectorTypeReal &XmX2_norm_frob, const VectorTypeReal &XmX2_trace)
Get homo and lumo bounds from traces and norms of Xi-Xi^2.
Definition purification_general.h:2257
double get_nnz_X(size_t &nnzX)
Get nnz of X in %.
Definition purification_general.h:377
IntervalType homo_bounds_F_new
Definition purification_general.h:540
VectorType eigVecHOMORef
Definition purification_general.h:566
PuriInfo info
Fill in during purification with useful information.
Definition purification_general.h:128
virtual void prepare_to_purification()
Prepare data for recursive expansion.
Definition purification_general.h:1065
void output_separate_total_times(PuriInfo &info) const
Definition purification_general.h:1582
virtual void purify_X(const int it)=0
void set_jump_over_X_iter_proj_method(int val)
Definition purification_general.h:257
void gen_matlab_file_norm_diff(const char *filename) const
Create MATLAB .m file which plots the idempotency error in each recursive expansion iteration.
Definition purification_general.h:3470
real error_per_it
Error allowed in each iteration due to truncation.
Definition purification_general.h:499
std::vector< real > eigValOCC
Here we save eigenvalues corresponding to the occupied orbitals.
Definition purification_general.h:570
VectorTypeInt VecPoly
Polynomials computed in the function estimated_number_of_iterations() VecPoly[i] = 1 if we use X=X^2 ...
Definition purification_general.h:508
IntervalType lumo_bounds
Lumo bounds for Xi in iteration i.
Definition purification_general.h:532
int eigenvectors_method
Chosen method for computing eigenvectors.
Definition purification_general.h:550
VectorTypeReal SIGMA_HOMO_VEC
Definition purification_general.h:517
intervalType IntervalType
Definition purification_general.h:120
virtual void get_iterations_for_lumo_and_homo(int &chosen_iter_lumo, int &chosen_iter_homo)
Find the best iterations for computing eigenvectors.
Definition purification_general.h:2899
real error_sub
Allowed error in invariant subspaces.
Definition purification_general.h:496
int get_exact_number_of_puri_iterations()
Definition purification_general.h:2241
void check_homo_lumo_eigenvalues(real &eigVal, VectorType &eigVec, bool &is_homo, bool &is_lumo, const int iter)
Definition purification_general.h:3356
int additional_iterations
Do a few more iterations after convergence.
Definition purification_general.h:479
int get_int_eig_iter_method(string eigenvectors_iterative_method)
Definition purification_general.h:886
int eigensolver_maxiter
Maximum number of iterations for eigensolver.
Definition purification_general.h:555
bool compute_eigenvectors_in_each_iteration
Compute homo and lumo eigenpairs in every iteration and save eigenvectors in txt files.
Definition purification_general.h:588
void projection_method_one_puri_iter(int current_iteration)
Definition purification_general.h:2757
virtual void prepare_to_purification_eigenvectors()
Prepare data related to the eigenvectors computations.
Definition purification_general.h:1011
MatrixType X
Matrix X.
Definition purification_general.h:130
VectorTypeInt good_iterations_homo
Iterations where homo eigenvector can be computed.
Definition purification_general.h:580
VectorTypeInt good_iterations_lumo
Iterations where homo eigenvector can be computed.
Definition purification_general.h:581
IntervalType homo_bounds_X0
Initial lumo bounds for X.
Definition purification_general.h:534
virtual void save_other_iter_info(IterationInfo &iter_info, int it)=0
virtual void estimate_number_of_iterations(int &estim_num_iter)=0
void set_spectrum_bounds(real eigmin, real eigmax)
Set spectrum bounds.
Definition purification_general.h:1954
int iter_for_lumo
Definition purification_general.h:578
void get_eigenvalue_of_F_from_eigv_of_Xi(real &eigVal, const VectorType &eigVec)
Definition purification_general.h:3346
void map_bounds_to_0_1()
Get eigenvalue bounds for X0.
Definition purification_general.h:1823
void compute_eigenvector(MatrixType const &M, std::vector< VectorType > &eigVec, std::vector< real > &eigVal, int it, bool is_homo)
Definition purification_general.h:3087
double get_nnz_Xsq(size_t &nnzXsq)
Get nnz of X^2 in %.
Definition purification_general.h:385
bool try_eigv_on_next_iteration_if_fail
Definition purification_general.h:563
void gen_matlab_file_threshold(const char *filename) const
Create MATLAB .m file which plots the actual introduced error after truncation of the matrix X_i in e...
Definition purification_general.h:3566
double get_nnz_X()
Get nnz of X in %.
Definition purification_general.h:381
void extract_computed_eigenpairs(std::vector< VectorType > &eigVecUNOCCref, std::vector< VectorType > &eigVecOCCref, std::vector< real > &eigValUNOCCref, std::vector< real > &eigValOCCref)
Definition purification_general.h:232
virtual bool is_initialized() const
Check is function initialize() is already called.
Definition purification_general.h:298
int jump_over_X_iter_proj_method
Definition purification_general.h:527
virtual void determine_iteration_for_eigenvectors()
Determine in which iterations will be computed homo and lumo eigenvectors.
Definition purification_general.h:2881
void save_matrix_now(string str)
Definition purification_general.h:2350
void compute_eigenvectors_without_diagonalization_last_iter_proj()
Definition purification_general.h:2677
bool puri_is_prepared_flag
Definition purification_general.h:476
void save_matrix_A_now(const MatrixType &A, string str)
Definition purification_general.h:2357
bool use_new_stopping_criterion
True for new stopping criterion.
Definition purification_general.h:478
virtual bool puri_is_prepared() const
Check is function prepare_to_purification() is already called.
Definition purification_general.h:302
virtual real apply_poly(const int it, real x)=0
IntervalType homo_bounds_F
Initial lumo bounds for F.
Definition purification_general.h:537
void compute_eigenvectors_without_diagonalization_on_F(const MatrixType &F, int eigensolver_maxiter_for_F)
Definition purification_general.h:2426
virtual void compute_X()
Get matrix X0 by mapping spectrum of F into [0,1] in reverse order.
Definition purification_general.h:1804
virtual void compute_spectrum_bounds()
Compute spectrum bounds.
Definition purification_general.h:1978
std::vector< MatrixType > vec_matrices_Xi
Save matrices Xi in each iteration (if used projection method for computing eigenvectors).
Definition purification_general.h:293
int get_jump_over_X_iter_proj_method() const
Definition purification_general.h:259
std::vector< VectorType > eigVecOCC
Here we save eigenvectors corresponding to the occupied orbitals.
Definition purification_general.h:568
virtual void discard_homo_eigenvector()
Definition purification_general.h:1668
void output_norms_and_traces(IterationInfo &iter_info) const
Definition purification_general.h:905
int nocc
Number of occupied orbitals.
Definition purification_general.h:484
int eigenvectors_iterative_method
Chosen eigensolver.
Definition purification_general.h:551
virtual void clear()
Clear all matrices in the class.
Definition purification_general.h:176
bool use_prev_vector_as_initial_guess
Definition purification_general.h:560
void output_time_WriteAndReadAll() const
Definition purification_general.h:1571
VectorTypeReal EIG_REL_GAP_LUMO_VEC
Definition purification_general.h:519
int number_of_unocc_eigenvectors
Definition purification_general.h:525
int scf_step
Definition purification_general.h:586
string eigenvectors_iterative_method_str
Definition purification_general.h:558
mat::normType NormType
Definition purification_general.h:121
int iter_for_homo
Definition purification_general.h:577
void set_eigenvectors_params_basic(string eigenvectors_method_, string eigenvectors_iterative_method_, real eigensolver_accuracy_, int eigensolver_maxiter_, int scf_step_, bool try_eigv_on_next_iteration_if_fail_, bool use_prev_vector_as_initial_guess_)
Set parameters for computing eigenvectors.
Definition purification_general.h:737
VectorType eigVecLUMORef
Definition purification_general.h:565
std::vector< VectorType > eigVecUNOCC
Here we save eigenvectors corresponding to the unoccupied orbitals.
Definition purification_general.h:569
virtual void check_standard_stopping_criterion(const real XmX2_norm, int &stop)
Check stopping criterion (obsolete).
Definition purification_general.h:2130
Time-measuring class.
Definition utilities.h:80
void print(int area, const char *routine)
Definition utilities.h:111
double get_elapsed_wall_seconds()
Definition utilities.h:107
virtual const char * what() const
Definition Failure.h:67
static std::string writeAndReadAll()
Definition FileWritable.cc:509
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
Definition VectorGeneral.h:48
void fullvector(std::vector< Treal > &fullVector) const
Definition VectorGeneral.h:88
File contataining definitions of constants used in the recursive expansion.
#define ORDER
Order of convergence of recursive expansion.
Definition constants.h:41
mat::SizesAndBlocks rows
Definition test.cc:51
mat::SizesAndBlocks cols
Definition test.cc:52
ergo_real real
Definition test.cc:46
int write_vector(const char *filename, const vector< real > &v)
Definition files_dense.cc:206
File containing declaration of functions for reading/writing dense matrices and vectors.
int write_matrix_to_mtx(const char *filename, const vector< int > &I, const vector< int > &J, const vector< real > &val, const int &N)
Definition files_sparse.cc:151
File containing declarations of functions for reading/writing sparse matrices from/to mtx (MatrixMark...
File containing declaration of functions for reading/writing sparse matrices from/to binary files.
Defined namespace eigvec containing functions for computing largest eigenvalues and corresponding eig...
Wrapper routines for different parts of the integral code, including conversion of matrices from/to t...
#define max(a, b)
Definition integrator.cc:87
int min(int a, int b)
Definition lin_trans.cc:66
Proxy structs used by the matrix API.
Header file with typedefs for matrix and vector types.
mat::Interval< ergo_real > intervalType
Definition matrix_typedefs.h:78
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
int computeEigenvectors(const MatrixType &A, Treal tol, std::vector< Treal > &eigVal, std::vector< VectorType > &eigVec, int number_of_eigenvalues_to_compute, std::string method, std::vector< int > &num_iter, int maxit=200, bool do_deflation=false)
Function for choosing method for computing eigenvectors.
Definition get_eigenvectors.h:232
normType
Definition matInclude.h:139
@ euclNorm
Definition matInclude.h:139
@ frobNorm
Definition matInclude.h:139
@ mixedNorm
Definition matInclude.h:139
@ zero
Definition matInclude.h:138
void output_current_memory_usage(int logArea, const char *contextString)
Definition output.cc:186
void enable_printf_output()
Definition output.cc:66
void do_output(int logCategory, int logArea, const char *format,...)
Definition output.cc:53
Functionality for writing output messages to a text file.
#define LOG_AREA_DENSFROMF
Definition output.h:61
#define LOG_CAT_INFO
Definition output.h:49
#define LOG_CAT_WARNING
Definition output.h:48
File containing classes IterationInfo and PuriInfo.
#define NUM_ADDITIONAL_ITERATIONS
Definition purification_general.h:76
int EIG_PROJECTION_INT
Definition purification_general.cc:55
real TOL_OVERLAPPING_BOUNDS
If the difference between inner bounds for homo and lumo is less then this tolerance,...
Definition purification_general.cc:45
int EIG_POWER_INT
Definition purification_general.cc:56
int EIG_LANCZOS_INT
Definition purification_general.cc:57
real THRESHOLD_EIG_TOLERANCE
Inner homo and lumo bounds may be too good, and it may happen that computed eigenvalue slightly outsi...
Definition purification_general.cc:49
ergo_real real
Definition purification_general.h:64
int EIG_SQUARE_INT
Definition purification_general.cc:54
real mixed_acc
Tolerance used for computation of mixed norm.
Definition purification_general.cc:43
real eucl_acc
Tolerance used for computation of spectral norm.
Definition purification_general.cc:42
int EIG_EMPTY
Definition purification_general.cc:53
intervalType IntervalType
Definition random_matrices.h:71
Definition of the main floating-point datatype used; the ergo_real type.
double ergo_real
Definition realtype.h:69
symmMatrix MatrixType
Definition recexp_diag_test.cc:66
Treal template_blas_sqrt(Treal x)
Treal template_blas_log(Treal x)
Treal template_blas_pow(Treal x, Treal y)
Treal template_blas_fabs(Treal x)
Constants for conversion between units for some common units like Angstrom, electron-volt (eV),...
#define UNIT_one_eV
Definition units.h:45
Basic OS access utilities.