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
115 template<
typename MatrixType>
159 bool use_new_stopping_criterion_,
197 {
return mat::getMachineEpsilon<real>(); }
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();
603 template<
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;
736 template<
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;
801 template<
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);
814 template<
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];
866 template<
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");
885 template<
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");
904 template<
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);
952 template<
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;
1010 template<
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();
1064 template<
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;
1099 template<
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);
1570 template<
typename MatrixType>
1573 #ifndef USE_CHUNKS_AND_TASKS
1581 template<
typename MatrixType>
1605 if(total_Xmixed >= 0)
1608 if(total_Xfrob >= 0)
1612 if(total_Xstcrit >= 0)
1616 if(total_Xtrace >= 0)
1624 template<
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();
1667 template<
typename MatrixType>
1671 info.eigValHOMO = 0;
1672 info.homo_eigenvector_is_computed =
false;
1676 template<
typename MatrixType>
1679 eigVecUNOCC.clear();
1680 info.eigValLUMO = 0;
1681 info.lumo_eigenvector_is_computed =
false;
1685 template<
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;
1803 template<
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));
1822 template<
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;
1881 template<
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);
1925 template<
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;
1953 template<
typename MatrixType>
1957 computed_spectrum_bounds =
true;
1963 template<
typename MatrixType>
1966 if (!computed_spectrum_bounds)
1968 compute_spectrum_bounds();
1970 eigmin = spectrum_bounds.low();
1971 eigmax = spectrum_bounds.upp();
1977 template<
typename MatrixType>
1981 real eigminG, eigmaxG, eigmin, eigmax;
1984 X.gershgorin(eigminG, eigmaxG);
1997 eigminG -= smallNumberToExpandWith;
1998 eigmaxG += smallNumberToExpandWith;
1999 #ifdef DEBUG_PURI_OUTPUT
2007 #if 1 // 0 - without Lanczos, 1 - with Lanczos
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;
2047 #endif // USE_CHUNKS_AND_TASKS
2051 if (eigmin == eigmax)
2060 computed_spectrum_bounds =
true;
2066 template<
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);
2129 template<
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
2149 template<
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))
2211 template<
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);
2233 template<
typename MatrixType>
2236 return info.estim_total_it;
2240 template<
typename MatrixType>
2243 if (info.converged == 1)
2245 return info.total_it;
2256 template<
typename MatrixType>
2261 estimate_homo_lumo(XmX2_norm_mixed, XmX2_norm_frob, XmX2_trace);
2269 template<
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]);
2349 template<
typename MatrixType>
2352 save_matrix_A_now(X, str);
2356 template<
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
2402 template<
typename MatrixType>
2405 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
2409 template<
typename MatrixType>
2412 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_last_iter_proj() is not implemented for CHTMatrix.");
2416 template<
typename MatrixType>
2419 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization() is not implemented for CHTMatrix.");
2425 template<
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]);
2505 template<
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;
2676 template<
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).");
2756 template<
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;
2876 #endif // USE_CHUNKS_AND_TASKS
2880 template<
typename MatrixType>
2888 get_iterations_for_lumo_and_homo(iter_for_lumo, iter_for_homo);
2898 template<
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.";
2961 template<
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;
2991 template<
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);
3012 template<
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
3078 template<
typename MatrixType>
3081 throw "compute_eigenvector() is not implemented with Chunks and Tasks.";
3086 template<
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++) {
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]);
3289 #endif // USE_CHUNKS_AND_TASKS
3292 template<
typename MatrixType>
3296 #ifndef USE_CHUNKS_AND_TASKS
3302 template<
typename MatrixType>
3306 #ifndef USE_CHUNKS_AND_TASKS
3312 template<
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.");
3345 template<
typename MatrixType>
3349 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3351 eigVal = approx_eig;
3355 template<
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;
3469 template<
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;
3565 template<
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;
3611 template<
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;
3709 template<
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;
3778 template<
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;
3981 #endif //HEADER_PURIFICATION_GENERAL