39 #ifndef HEADER_PURIFICATION_GENERAL 40 #define HEADER_PURIFICATION_GENERAL 76 #define NUM_ADDITIONAL_ITERATIONS 0 80 #define PURI_OUTPUT_NNZ 110 template<
typename MatrixType>
143 int use_new_stopping_criterion_,
199 {
return mat::getMachineEpsilon<real>(); }
229 string eigenvectors_iterative_method_,
230 real eigensolver_accuracy_,
231 int eigensolver_maxiter_,
233 int use_prev_vector_as_initial_guess_,
234 int try_eigv_on_next_iteration_if_fail_,
284 int& stop,
real& estim_order);
307 { nnzX =
X.nnz();
return (
double)(((double)nnzX) / ((double)
X.get_ncols() *
X.get_nrows()) * 100); }
311 {
return (
double)(((double)
X.nnz()) / ((
double)
X.get_ncols() *
X.get_nrows()) * 100); }
315 { nnzXsq =
Xsq.nnz();
return (
double)(((double)nnzXsq) / ((double)
Xsq.get_ncols() *
Xsq.get_nrows()) * 100); }
319 {
return (
double)(((double)
Xsq.nnz()) / ((
double)
Xsq.get_ncols() *
Xsq.get_nrows()) * 100); }
331 const std::vector<real>& h_in,
332 const std::vector<real>& l_in,
351 int& chosen_iter_homo);
383 virtual void purify_X(
const int it) = 0;
520 template<
typename MatrixType>
527 int use_new_stopping_criterion_,
536 error_sub = error_sub_;
537 error_eig = error_eig_;
538 use_new_stopping_criterion = use_new_stopping_criterion_;
539 normPuriTrunc = norm_truncation_;
540 normPuriStopCrit = norm_stop_crit_;
543 initialized_flag =
true;
546 lumo_bounds_F = lumo_bounds_;
547 homo_bounds_F = homo_bounds_;
551 #ifdef ENABLE_PRINTF_OUTPUT 556 double nnzXp = get_nnz_X(nnzX);
558 " , nocc = %d , NNZ = %lu <-> %.5lf %%",
559 X.get_nrows(), nocc, nnzX, nnzXp);
563 switch (normPuriTrunc)
578 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
582 #ifdef USE_CHUNKS_AND_TASKS 591 switch (normPuriStopCrit)
606 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
610 #ifdef USE_CHUNKS_AND_TASKS 618 if (this->use_new_stopping_criterion == 1)
629 check_stopping_criterion_iter = -1;
630 compute_eigenvectors_in_this_SCF_cycle =
false;
635 info.stopping_criterion = use_new_stopping_criterion;
636 info.error_subspace = error_sub;
638 info.debug_output = 0;
639 #ifdef DEBUG_PURI_OUTPUT 640 info.debug_output = 1;
648 template<
typename MatrixType>
650 string eigenvectors_iterative_method_,
651 real eigensolver_accuracy_,
652 int eigensolver_maxiter_,
654 int use_prev_vector_as_initial_guess_,
655 int try_eigv_on_next_iteration_if_fail_,
661 number_of_eigenvalues = 1;
662 if (number_of_eigenvalues > 1)
664 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot compute more than 1 eigenpair.");
668 eigVecLUMO = eigVecLUMO_;
669 eigVecHOMO = eigVecHOMO_;
671 eigensolver_accuracy = eigensolver_accuracy_;
672 eigensolver_maxiter = eigensolver_maxiter_;
673 assert(eigensolver_maxiter >= 1);
674 scf_step = scf_step_;
675 eigenvectors_method_str = eigenvectors_method_;
676 eigenvectors_iterative_method_str = eigenvectors_iterative_method_;
677 eigenvectors_method = get_int_eig_method(eigenvectors_method_);
678 eigenvectors_iterative_method = get_int_eig_iter_method(eigenvectors_iterative_method_);
679 use_prev_vector_as_initial_guess = use_prev_vector_as_initial_guess_;
680 try_eigv_on_next_iteration_if_fail = try_eigv_on_next_iteration_if_fail_;
682 info.lumo_eigenvector_is_computed =
false;
683 info.homo_eigenvector_is_computed =
false;
689 if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && (eigenvectors_method ==
EIG_EMPTY))
693 "Eigenvectors will not be computed in this SCF cycle. Set eigenvectors to NULL.");
708 if (((eigVecLUMO != NULL) || (eigVecHOMO != NULL)) && use_prev_vector_as_initial_guess)
714 #ifndef USE_CHUNKS_AND_TASKS 715 if (compute_eigenvectors_in_each_iteration)
718 if ((eigVecLUMO != NULL) && use_prev_vector_as_initial_guess)
723 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for LUMO!");
726 eigVecLUMORef.resetSizesAndBlocks(
cols);
727 eigVecLUMORef = *eigVecLUMO;
730 if ((eigVecHOMO != NULL) && use_prev_vector_as_initial_guess)
735 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for HOMO!");
738 eigVecHOMORef.resetSizesAndBlocks(
cols);
739 eigVecHOMORef = *eigVecHOMO;
746 template<
typename MatrixType>
749 if (eigenvectors_method ==
"square")
753 if (eigenvectors_method ==
"projection")
757 if (eigenvectors_method ==
"")
761 throw std::runtime_error(
"Error in get_int_eig_method(): unknown method to compute eigenvectors");
765 template<
typename MatrixType>
768 if (eigenvectors_iterative_method ==
"power")
772 if (eigenvectors_iterative_method ==
"lanczos")
776 if (eigenvectors_iterative_method ==
"")
780 throw std::runtime_error(
"Error in get_int_eig_iter_method(): unknown iterative method to compute eigenvectors");
784 template<
typename MatrixType>
791 int it = iter_info.
it;
795 #ifndef USE_CHUNKS_AND_TASKS 798 assert(XmX2_eucl >= 0);
804 assert(XmX2_fro_norm >= 0);
805 assert(XmX2_mixed_norm >= 0);
812 assert(XmX2_fro_norm >= 0);
824 template<
typename MatrixType>
829 prepare_to_purification();
831 #ifdef DEBUG_PURI_OUTPUT 834 purification_process();
836 eigenvalue_bounds_estimation();
841 if (info.converged == 1)
843 if (compute_eigenvectors_in_this_SCF_cycle && (eigenvectors_method ==
EIG_PROJECTION_INT))
845 compute_eigenvectors_without_diagonalization_last_iter_proj();
854 check_eigenvectors_at_the_end();
858 info.total_time = total_time_stop;
862 template<
typename MatrixType>
865 if (!is_initialized())
867 throw std::runtime_error(
"Error in prepare_to_purification() : function is called for an uninitialized class.");
870 if (!computed_spectrum_bounds)
873 compute_spectrum_bounds();
876 info.time_spectrum_bounds = total_time_spectrum_bounds_stop;
879 info.set_spectrum_bounds(spectrum_bounds.low(), spectrum_bounds.upp());
883 set_truncation_parameters();
885 if ((eigVecLUMO != NULL) || (eigVecHOMO != NULL))
890 compute_eigenvectors_in_this_SCF_cycle =
true;
891 info.compute_eigenvectors_in_this_SCF_cycle = compute_eigenvectors_in_this_SCF_cycle;
897 determine_iteration_for_eigenvectors();
908 puri_is_prepared_flag =
true;
912 template<
typename MatrixType>
915 if (!puri_is_prepared())
917 throw std::runtime_error(
"Error in purification_process() : " 918 "first expect a call for function prepare_to_purification().");
924 real Xsquare_time_stop = -1, total_time_stop = -1, trunc_time_stop = -1, purify_time_stop = -1;
925 real frob_diff_time_stop = -1, eucl_diff_time_stop = -1, trace_diff_time_stop = -1, mixed_diff_time_stop = -1, stopping_criterion_time_stop = -1;
926 int maxit_tmp = maxit;
927 real estim_order = -1;
928 real XmX2_trace = -1;
929 real XmX2_fro_norm = -1;
930 real XmX2_mixed_norm = -1;
933 int already_stop_before = 0;
935 info.Iterations.clear();
936 info.Iterations.reserve(100);
949 #ifdef SAVE_MATRIX_IN_EACH_ITERATION 953 save_matrix_now(str.str());
959 #ifdef PURI_OUTPUT_NNZ 960 double nnzX = get_nnz_X();
965 truncate_matrix(thresh, it);
969 #ifdef PURI_OUTPUT_NNZ 970 if((
double)thresh >= 0)
972 double nnzXafter = get_nnz_X();
977 double nnzXafter = get_nnz_X();
981 if((
double)thresh >= 0)
988 Xsq = (
real)1.0 * X * X;
992 #ifdef PURI_OUTPUT_NNZ 994 double nnzXsqp = get_nnz_Xsq(nnzXsq);
1001 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1007 #ifndef USE_CHUNKS_AND_TASKS 1009 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq,
mixed_acc);
1017 XmX2_trace = X.trace() - Xsq.trace();
1023 #ifndef USE_CHUNKS_AND_TASKS 1025 XmX2_eucl = MatrixType::eucl_diff(X, Xsq,
eucl_acc);
1037 #ifdef DEBUG_PURI_OUTPUT 1038 output_norms_and_traces(iter_info);
1042 if (compute_eigenvectors_in_this_SCF_cycle)
1044 compute_eigenvectors_without_diagonalization(it, iter_info);
1047 ostringstream str_out;
1048 str_out <<
"Iteration " << it;
1054 iter_info.
gap = 1 - homo_bounds.upp() - lumo_bounds.upp();
1059 iter_info.
NNZ_X = get_nnz_X();
1060 iter_info.
NNZ_X2 = get_nnz_Xsq();
1066 stopping_criterion_time_stop = 0;
1073 save_other_iter_info(iter_info, it);
1077 info.Iterations.push_back(iter_info);
1088 while (it <= maxit_tmp)
1102 #ifdef SAVE_MATRIX_IN_EACH_ITERATION 1106 save_matrix_now(str.str());
1110 #ifdef PURI_OUTPUT_NNZ 1111 double nnzX = get_nnz_X();
1116 truncate_matrix(thresh, it);
1120 #ifdef PURI_OUTPUT_NNZ 1121 double nnzXafter = get_nnz_X();
1131 Xsq = (
real)1.0 * X * X;
1134 #ifdef PURI_OUTPUT_NNZ 1136 double nnzXsqp = get_nnz_Xsq(nnzXsq);
1146 XmX2_fro_norm = MatrixType::frob_diff(X, Xsq);
1152 #ifndef USE_CHUNKS_AND_TASKS 1154 XmX2_mixed_norm = MatrixType::mixed_diff(X, Xsq,
mixed_acc);
1163 #ifndef USE_CHUNKS_AND_TASKS 1165 XmX2_eucl = MatrixType::eucl_diff(X, Xsq,
eucl_acc);
1173 XmX2_trace = X.trace() - Xsq.trace();
1179 if (XmX2_trace < -1e10)
1181 throw std::runtime_error(
"Error in purification_process() : trace of X-X^2 is negative, seems as a" 1182 " misconvergence of the recursive expansion.");
1190 #ifdef DEBUG_PURI_OUTPUT 1191 output_norms_and_traces(iter_info);
1210 if (compute_eigenvectors_in_this_SCF_cycle)
1212 compute_eigenvectors_without_diagonalization(it, iter_info);
1218 if (it >= check_stopping_criterion_iter)
1221 stopping_criterion(iter_info, stop, estim_order);
1234 if ((already_stop_before == 0) && ((stop == 1) || (it == maxit)))
1247 assert(maxit_tmp <= (
int)VecPoly.size());
1248 maxit_tmp = it + additional_iterations;
1249 already_stop_before = 1;
1252 ostringstream str_out;
1253 str_out <<
"Iteration " << it;
1262 iter_info.
NNZ_X = get_nnz_X();
1263 iter_info.
NNZ_X2 = get_nnz_Xsq();
1272 if (use_new_stopping_criterion)
1274 iter_info.
order = estim_order;
1277 save_other_iter_info(iter_info, it);
1281 info.Iterations.push_back(iter_info);
1288 double nnzD = get_nnz_X();
1299 info.total_it = maxit_tmp;
1300 info.additional_iterations = additional_iterations;
1302 real acc_err_sub = total_subspace_error(maxit_tmp - additional_iterations);
1303 #ifdef DEBUG_PURI_OUTPUT 1304 if (acc_err_sub != -1)
1309 info.accumulated_error_subspace = acc_err_sub;
1315 template<
typename MatrixType>
1318 if (info.converged == 1)
1325 info.get_vec_mixed_norms(norms_mixed);
1327 info.get_vec_frob_norms(norms_frob);
1328 info.get_vec_traces(traces);
1329 get_eigenvalue_estimates(norms_mixed, norms_frob, traces);
1342 info.homo_estim_low_F = homo_bounds_F_new.low();
1343 info.homo_estim_upp_F = homo_bounds_F_new.upp();
1344 info.lumo_estim_low_F = lumo_bounds_F_new.low();
1345 info.lumo_estim_upp_F = lumo_bounds_F_new.upp();
1352 template<
typename MatrixType>
1355 if ((eigVecHOMO == NULL) || eigVecHOMO->is_empty())
1359 info.eigValHOMO = 0;
1360 info.homo_eigenvector_is_computed =
false;
1364 eigVecHOMO->clear_structure();
1368 template<
typename MatrixType>
1371 if ((eigVecLUMO == NULL) || eigVecLUMO->is_empty())
1375 info.eigValLUMO = 0;
1376 info.lumo_eigenvector_is_computed =
false;
1377 eigVecLUMO->clear_structure();
1381 template<
typename MatrixType>
1384 if (info.converged != 1)
1387 discard_homo_eigenvector();
1388 discard_lumo_eigenvector();
1394 if (compute_eigenvectors_in_this_SCF_cycle)
1398 if (!eigVecHOMO->is_empty() && (info.total_it < iter_for_homo))
1401 iter_for_homo, info.total_it);
1402 discard_homo_eigenvector();
1407 if (!eigVecHOMO->is_empty() && (info.total_it - iter_for_homo < ITER_DIFF) && info.homo_eigenvector_is_computed)
1410 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, " 1411 "thus result may not be accurate!", iter_for_homo, info.total_it);
1416 if (!eigVecLUMO->is_empty() && (info.total_it < iter_for_lumo))
1419 iter_for_lumo, info.total_it);
1420 discard_lumo_eigenvector();
1425 if (!eigVecLUMO->is_empty() && (info.total_it - iter_for_lumo < ITER_DIFF) && info.lumo_eigenvector_is_computed)
1428 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, " 1429 "thus result may not be accurate!", iter_for_lumo, info.total_it);
1435 real low_lumo_F_bound = info.lumo_estim_low_F;
1436 real upp_lumo_F_bound = info.lumo_estim_upp_F;
1437 real low_homo_F_bound = info.homo_estim_low_F;
1438 real upp_homo_F_bound = info.homo_estim_upp_F;
1444 if (info.homo_eigenvector_is_computed)
1446 if ((low_homo_F_bound - flex_tolerance <= eigValHOMO) && (eigValHOMO <= upp_homo_F_bound + flex_tolerance))
1449 (
double)eigValHOMO, (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
1450 info.eigValHOMO = eigValHOMO;
1455 " discard computed HOMO eigenvector.",
1456 (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
1457 discard_homo_eigenvector();
1462 discard_homo_eigenvector();
1466 if (info.lumo_eigenvector_is_computed)
1468 if ((low_lumo_F_bound - flex_tolerance <= eigValLUMO) && (eigValLUMO <= upp_lumo_F_bound + flex_tolerance))
1471 (
double)eigValLUMO, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
1472 info.eigValLUMO = eigValLUMO;
1477 " discard computed LUMO eigenvector.",
1478 (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
1479 discard_lumo_eigenvector();
1484 discard_lumo_eigenvector();
1488 if (info.homo_eigenvector_is_computed && info.lumo_eigenvector_is_computed)
1490 real gap = eigValLUMO - eigValHOMO;
1499 template<
typename MatrixType>
1502 if (spectrum_bounds.empty())
1507 #ifdef DEBUG_PURI_OUTPUT 1511 real eigmin = spectrum_bounds.low();
1512 real eigmax = spectrum_bounds.upp();
1513 X.add_identity(-eigmax);
1514 X *= ((
real)1.0 / (eigmin - eigmax));
1518 template<
typename MatrixType>
1521 if (spectrum_bounds.empty())
1526 #ifdef DEBUG_PURI_OUTPUT 1530 real eigmin = spectrum_bounds.low();
1531 real eigmax = spectrum_bounds.upp();
1535 homo_bounds = homo_bounds_F;
1536 lumo_bounds = lumo_bounds_F;
1538 homo_bounds.intersect(spectrum_bounds);
1539 lumo_bounds.intersect(spectrum_bounds);
1541 #ifdef DEBUG_PURI_OUTPUT 1542 if (homo_bounds.empty())
1546 if (lumo_bounds.empty())
1559 homo_bounds = (homo_bounds - eigmax) / (eigmin - eigmax);
1560 homo_bounds_X0 = homo_bounds;
1561 homo_bounds =
IntervalType(1 - homo_bounds.upp(), 1 - homo_bounds.low());
1565 lumo_bounds = (lumo_bounds - eigmax) / (eigmin - eigmax);
1566 lumo_bounds_X0 = lumo_bounds;
1577 template<
typename MatrixType>
1580 int estim_num_iter = 0;
1582 estimate_number_of_iterations(estim_num_iter);
1583 info.estim_total_it = estim_num_iter;
1587 if (estim_num_iter < maxit)
1590 maxit = estim_num_iter;
1594 error_per_it = error_sub / estim_num_iter;
1605 template<
typename MatrixType>
1609 computed_spectrum_bounds =
true;
1615 template<
typename MatrixType>
1618 if (!computed_spectrum_bounds)
1620 compute_spectrum_bounds();
1622 eigmin = spectrum_bounds.low();
1623 eigmax = spectrum_bounds.upp();
1629 template<
typename MatrixType>
1633 real eigminG, eigmaxG, eigmin, eigmax;
1636 X.gershgorin(eigminG, eigmaxG);
1649 eigminG -= smallNumberToExpandWith;
1650 eigmaxG += smallNumberToExpandWith;
1651 #ifdef DEBUG_PURI_OUTPUT 1659 #if 1 // 0 - without Lanczos, 1 - with Lanczos 1660 #ifndef USE_CHUNKS_AND_TASKS 1667 Xshifted.add_identity((
real)(-1.0) * eigminG);
1672 eigmax = Xshifted.eucl(acc, maxIter) + eigminG + acc;
1685 Xshifted.add_identity((
real)(1.0) * eigminG);
1686 Xshifted.add_identity((
real)(-1.0) * eigmaxG);
1690 eigmin = -Xshifted.eucl(acc, maxIter) + eigmaxG - acc;
1699 #endif // USE_CHUNKS_AND_TASKS 1703 if (eigmin == eigmax)
1712 computed_spectrum_bounds =
true;
1718 template<
typename MatrixType>
1721 int it = iter_info.
it;
1722 real XmX2_norm_it = -1, XmX2_norm_itm2 = -1;
1724 if (it < check_stopping_criterion_iter)
1737 if(use_new_stopping_criterion && check_stopping_criterion_iter == -1 && it < 2)
1739 if (use_new_stopping_criterion)
1746 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_eucl;
1752 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_fro_norm;
1758 XmX2_norm_itm2 = info.Iterations[it - 2].XmX2_mixed_norm;
1762 throw std::runtime_error(
"Error in stopping_criterion() : unknown matrix norm.");
1766 check_new_stopping_criterion(it, XmX2_norm_it, XmX2_norm_itm2, XmX2_trace, stop, estim_order);
1784 check_standard_stopping_criterion(XmX2_norm_it, stop);
1789 template<
typename MatrixType>
1795 bool homoLumo_converged = (homo_bounds.upp() < error_eig &&
1796 lumo_bounds.upp() < error_eig);
1797 bool XmX2norm_converged = XmX2_norm < error_eig;
1798 if (homoLumo_converged || XmX2norm_converged)
1803 #ifdef DEBUG_PURI_OUTPUT 1809 template<
typename MatrixType>
1811 int& stop,
real& estim_order)
1825 return_constant_C(it, C);
1826 this->constant_C = C;
1835 if ((VecPoly[it - 1] != VecPoly[it]) &&
1836 (XmX2_norm_itm2 < 1) &&
1853 #ifdef DEBUG_PURI_OUTPUT 1854 if ((VecPoly[it - 1] != VecPoly[it]) && (XmX2_norm_itm2 < 1))
1869 template<
typename MatrixType>
1873 assert(it <= (
int)VecGap.size());
1877 for (
int i = 0; i <= it; ++i)
1879 if (VecGap[i] == -1)
1883 normE = info.Iterations[i].threshold_X;
1884 error += normE / (VecGap[i] - normE);
1891 template<
typename MatrixType>
1894 return info.estim_total_it;
1898 template<
typename MatrixType>
1901 if (info.converged == 1)
1903 return info.total_it;
1913 template<
typename MatrixType>
1920 out_occ.resize(nmax);
1922 out_unocc.resize(nmax);
1924 out_unocc[0] = value_unocc;
1925 out_occ[0] = value_occ;
1928 for (
int i = 1; i < nmax; ++i)
1930 occ = out_occ[i - 1];
1931 unocc = out_unocc[i - 1];
1932 apply_poly_to_eigs(i, occ, unocc);
1934 out_unocc[i] = unocc;
1939 template<
typename MatrixType>
1944 estimate_homo_lumo(XmX2_norm_mixed, XmX2_norm_frob, XmX2_trace);
1948 real eigmax = spectrum_bounds.upp();
1949 real eigmin = spectrum_bounds.low();
1950 real homo_in_0 = 1 - (homo_bounds_F_new.low() - eigmax) / (eigmin - eigmax);
1951 real lumo_in_0 = (lumo_bounds_F_new.upp() - eigmax) / (eigmin - eigmax);
1953 int n = get_exact_number_of_puri_iterations();
1956 propagate_values_in_each_iter(lumo_in_0, homo_in_0, l_in, h_in, n);
1959 get_frob_norm_est(XmX2_norm_frob, h_in, l_in, YmY2_norm_frob_est);
1963 template<
typename MatrixType>
1969 int n = get_exact_number_of_puri_iterations();
1973 YmY2_norm_frob_est.clear();
1974 YmY2_norm_frob_est.resize(n);
1977 for (
int i = 0; i < n; ++i)
1979 th = h_in[i] - h_in[i] * h_in[i];
1980 tl = l_in[i] - l_in[i] * l_in[i];
1981 YmY2_norm_frob_est[i] =
template_blas_sqrt(XmX2_norm_frob[i] * XmX2_norm_frob[i] - th * th - tl * tl);
1991 template<
typename MatrixType>
2000 int total_it = info.total_it - info.additional_iterations;
2007 real STOP_NORM = gammaStopEstim - gammaStopEstim * gammaStopEstim;
2012 if (XmX2_norm_mixed.size() == XmX2_norm_frob.size())
2014 XmX2_norm_out = XmX2_norm_mixed;
2018 XmX2_norm_out = XmX2_norm_frob;
2021 for (
int it = total_it; it >= 0; it--)
2023 vi = XmX2_norm_frob[it];
2024 wi = XmX2_trace[it];
2025 mi = XmX2_norm_out[it];
2027 if (vi >= STOP_NORM)
2038 temp_value = vi * vi / wi;
2043 bounds_from_it[2] = bounds_from_it[0];
2044 bounds_from_it[3] = bounds_from_it[1];
2047 apply_inverse_poly_vector(it, bounds_from_it);
2050 final_bounds[0] =
std::min(final_bounds[0], bounds_from_it[0]);
2051 final_bounds[1] =
std::min(final_bounds[1], bounds_from_it[1]);
2053 final_bounds[2] =
std::min(final_bounds[2], bounds_from_it[2]);
2054 final_bounds[3] =
std::min(final_bounds[3], bounds_from_it[3]);
2058 real maxeig = spectrum_bounds.upp();
2059 real mineig = spectrum_bounds.low();
2060 lumo_bounds_F_new =
IntervalType(maxeig * (1 - final_bounds[1]) + mineig * final_bounds[1],
2061 maxeig * (1 - final_bounds[0]) + mineig * final_bounds[0]);
2062 homo_bounds_F_new =
IntervalType(mineig * (1 - final_bounds[2]) + maxeig * final_bounds[2],
2063 mineig * (1 - final_bounds[3]) + maxeig * final_bounds[3]);
2069 template<
typename MatrixType>
2072 #ifdef USE_CHUNKS_AND_TASKS 2073 vector<int> Itmp, I, Jtmp, J;
2074 vector<real> Vtmp, V;
2075 X.get_all_values(Itmp, Jtmp, Vtmp);
2079 for (
size_t i = 0; i < Itmp.size(); i++)
2081 nnz += (Vtmp[i] != 0);
2088 for (
size_t i = 0; i < Itmp.size(); i++)
2092 I.push_back(Itmp[i]);
2093 J.push_back(Jtmp[i]);
2094 V.push_back(Vtmp[i]);
2098 string name =
"X_" + str +
".mtx";
2101 throw std::runtime_error(
"Error in save_matrix_now : error in write_matrix_to_mtx.\n");
2113 #ifdef USE_CHUNKS_AND_TASKS 2115 template<
typename MatrixType>
2118 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_on_F() is not implemented for CHTMatrix.");
2122 template<
typename MatrixType>
2125 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization_last_iter_proj() is not implemented for CHTMatrix.");
2129 template<
typename MatrixType>
2132 throw std::runtime_error(
"compute_eigenvectors_without_diagonalization() is not implemented for CHTMatrix.");
2138 template<
typename MatrixType>
2141 real start_shift = homo_bounds_F.upp();
2142 real end_shift = lumo_bounds_F.low();
2143 real dist = (end_shift - start_shift) / 15;
2144 real acc_eigv = eigensolver_accuracy;
2150 Fsq = (
real)1.0 * F * F;
2155 for (sigma = start_shift; sigma < end_shift; sigma += dist)
2158 M.add_identity(sigma * sigma);
2159 M += ((
real) - 2 * sigma) * F;
2160 M = ((
real) - 1.0) * M;
2162 vector<real> eigValTmp(1);
2163 vector<int> num_iter_solver(1, -1);
2170 eigenvectors_iterative_method_str, num_iter_solver,
2171 eigensolver_maxiter_for_F);
2173 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2175 printf(
"sigma = %lf , eigval = %lf , iters = %d\n", (
double)sigma, (
double)eigval, num_iter_solver[0]);
2181 M.add_identity(sigma * sigma);
2182 M += ((
real) - 2 * sigma) * F;
2183 M = ((
real) - 1.0) * M;
2185 vector<real> eigValTmp(1);
2186 vector<int> num_iter_solver(1, -1);
2193 eigenvectors_iterative_method_str, num_iter_solver,
2194 eigensolver_maxiter_for_F);
2196 eigval = eigvec::compute_rayleigh_quotient<real>(F, eigVec[0]);
2198 printf(
"sigma = %lf , eigval = %lf , iters = %d\n", (
double)sigma, (
double)eigval, num_iter_solver[0]);
2203 template<
typename MatrixType>
2206 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2209 bool compute_homo_now =
false;
2210 bool compute_lumo_now =
false;
2213 if (compute_eigenvectors_in_each_iteration)
2218 if (eigVecHOMO != NULL)
2220 for (
size_t i = 0; i < good_iterations_homo.size(); ++i)
2222 if (good_iterations_homo[i] == it)
2224 compute_homo_now =
true;
2231 if (eigVecLUMO != NULL)
2233 for (
size_t i = 0; i < good_iterations_lumo.size(); ++i)
2235 if (good_iterations_lumo[i] == it)
2237 compute_lumo_now =
true;
2246 if (eigVecHOMO != NULL)
2248 compute_homo_now =
true;
2250 if (eigVecLUMO != NULL)
2252 compute_lumo_now =
true;
2259 if (eigVecHOMO != NULL)
2261 if (it == iter_for_homo)
2263 compute_homo_now =
true;
2266 if (eigVecLUMO != NULL)
2268 if (it == iter_for_lumo)
2270 compute_lumo_now =
true;
2275 if (compute_homo_now && !info.homo_eigenvector_is_computed)
2282 writeToTmpFile(Xsq);
2284 real sigma = SIGMA_HOMO_VEC[it];
2289 M.add_identity(sigma * sigma);
2290 M += ((
real) - 2 * sigma) * X;
2291 M = ((
real) - 1.0) * M;
2294 compute_eigenvector(M, eigVecHOMO, it,
true);
2306 readFromTmpFile(Xsq);
2310 if (compute_eigenvectors_in_each_iteration)
2314 name <<
"homo_" << it;
2316 save_matrix_now(name.str());
2333 if (compute_lumo_now && !info.lumo_eigenvector_is_computed)
2340 writeToTmpFile(Xsq);
2342 real sigma = SIGMA_LUMO_VEC[it];
2347 M.add_identity(sigma * sigma);
2348 M += ((
real) - 2 * sigma) * X;
2349 M = ((
real) - 1.0) * M;
2352 compute_eigenvector(M, eigVecLUMO, it,
false);
2363 readFromTmpFile(Xsq);
2367 if (compute_eigenvectors_in_each_iteration)
2371 name <<
"lumo_" << it;
2373 save_matrix_now(name.str());
2389 if (compute_eigenvectors_in_each_iteration)
2392 info.homo_eigenvector_is_computed =
false;
2393 info.lumo_eigenvector_is_computed =
false;
2398 template<
typename MatrixType>
2401 real homo_total_time_stop, lumo_total_time_stop, homo_solver_time_stop, lumo_solver_time_stop;
2402 real DX_mult_time_homo_stop, DX_mult_time_lumo_stop;
2412 assert(eigVecHOMO != NULL);
2413 assert(eigVecLUMO != NULL);
2415 if (compute_eigenvectors_in_each_iteration)
2417 int total_it = info.total_it;
2423 for (
int it = 0; it <= total_it; ++it)
2428 name <<
"X_lumo_" << it <<
".mtx";
2443 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP, Xi, 0, DXi);
2452 compute_eigenvector(Zh, eigVecHOMO, it,
true);
2463 DXi = (
real)(-1) * DXi;
2466 compute_eigenvector(DXi, eigVecLUMO, it,
false);
2473 info.Iterations[it].DX_mult_homo_time = DX_mult_time_homo_stop;
2474 info.Iterations[it].DX_mult_lumo_time = 0;
2476 info.Iterations[it].homo_eig_solver_time = homo_solver_time_stop;
2477 info.Iterations[it].lumo_eig_solver_time = lumo_solver_time_stop;
2479 info.Iterations[it].orbital_homo_time = homo_total_time_stop;
2480 info.Iterations[it].orbital_lumo_time = lumo_total_time_stop;
2489 if (info.total_it >= iter_for_homo)
2493 readFromTmpFile(X_homo);
2505 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP, X_homo, 0, DXi);
2514 compute_eigenvector(Zh, eigVecHOMO, iter_for_homo,
true);
2518 info.Iterations[iter_for_homo].homo_eig_solver_time = homo_solver_time_stop;
2519 info.Iterations[iter_for_homo].DX_mult_homo_time = DX_mult_time_homo_stop;
2524 info.Iterations[iter_for_homo].orbital_homo_time = homo_total_time_stop;
2528 if (iter_for_homo == iter_for_lumo)
2536 DXi = (
real)(-1) * DXi;
2539 compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo,
false);
2547 info.Iterations[iter_for_lumo].DX_mult_lumo_time = 0;
2548 info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop;
2549 info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
2557 if ((info.total_it >= iter_for_lumo) && (iter_for_homo != iter_for_lumo))
2560 readFromTmpFile(X_lumo);
2570 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP, X_lumo, 0, DXi);
2576 DXi = (
real)(-1) * DXi;
2579 compute_eigenvector(DXi, eigVecLUMO, iter_for_lumo,
false);
2582 info.Iterations[iter_for_lumo].lumo_eig_solver_time = lumo_solver_time_stop;
2583 info.Iterations[iter_for_lumo].DX_mult_lumo_time = DX_mult_time_lumo_stop;
2587 info.Iterations[iter_for_lumo].orbital_lumo_time = lumo_total_time_stop;
2594 #endif // USE_CHUNKS_AND_TASKS 2598 template<
typename MatrixType>
2603 get_iterations_for_lumo_and_homo(iter_for_lumo, iter_for_homo);
2617 throw std::runtime_error(
"Error in determine_iteration_for_eigenvectors: unknown method for computing eigenvectors.");
2627 template<
typename MatrixType>
2629 int& chosen_iter_homo)
2631 int maxiter = maxit;
2634 real homo0 = 1 - homo_bounds.upp();
2635 real lumo0 = lumo_bounds.upp();
2636 real homoi = homo0, lumoi = lumo0;
2638 real Dh_homo, Dh_lumo, Dgh_homo, Dgh_lumo,
2639 Dgh_homo_max = get_min_double(), Dgh_lumo_max = get_min_double();
2641 chosen_iter_lumo = -1;
2642 chosen_iter_homo = -1;
2644 good_iterations_homo.clear();
2645 good_iterations_lumo.clear();
2647 find_shifts_every_iter();
2649 for (
int i = 1; i <= maxiter; ++i)
2651 homoi = apply_poly(i, homoi);
2652 lumoi = apply_poly(i, lumoi);
2654 Dh_homo = compute_derivative(i, homo0, dummy);
2655 Dh_lumo = compute_derivative(i, lumo0, dummy);
2658 Dgh_homo = 2 * (homoi - SIGMA_HOMO_VEC[i]) * Dh_homo;
2659 Dgh_lumo = 2 * (lumoi - SIGMA_LUMO_VEC[i]) * Dh_lumo;
2661 if (homoi >= SIGMA_HOMO_VEC[i])
2663 good_iterations_homo.push_back(i);
2664 if (Dgh_homo >= Dgh_homo_max)
2666 Dgh_homo_max = Dgh_homo;
2667 chosen_iter_homo = i;
2671 if (lumoi <= SIGMA_LUMO_VEC[i])
2673 good_iterations_lumo.push_back(i);
2677 chosen_iter_lumo = i;
2682 if ((chosen_iter_homo == -1) || (chosen_iter_lumo == -1))
2684 throw "Error in get_iterations_for_lumo_and_homo() : something went wrong, cannot choose iteration to compute HOMO or LUMO eigenvector.";
2689 template<
typename MatrixType>
2692 int maxiter = this->maxit;
2694 this->ITER_ERROR_VEC.clear();
2695 this->ITER_ERROR_VEC.resize(maxiter + 1);
2697 real error = error_per_it;
2698 if (error_per_it == 0)
2700 error = this->get_epsilon();
2703 for (
int i = 1; i <= maxiter; ++i)
2705 this->ITER_ERROR_VEC[i] = (error * this->VecGap[i]) / (1 + error);
2713 template<
typename MatrixType>
2716 int maxiter = maxit;
2718 SIGMA_HOMO_VEC.resize(maxiter + 1);
2719 SIGMA_LUMO_VEC.resize(maxiter + 1);
2723 real homo = 1 - homo_bounds.upp();
2724 real lumo = lumo_bounds.upp();
2725 real homo_out = 1 - homo_bounds.low();
2726 real lumo_out = lumo_bounds.low();
2728 for (
int i = 1; i <= maxiter; ++i)
2730 homo = apply_poly(i, homo);
2731 lumo = apply_poly(i, lumo);
2733 homo_out = apply_poly(i, homo_out);
2734 lumo_out = apply_poly(i, lumo_out);
2736 SIGMA_HOMO_VEC[i] = (homo_out + lumo) / 2;
2737 SIGMA_LUMO_VEC[i] = (lumo_out + homo) / 2;
2742 template<
typename MatrixType>
2745 int maxiter = maxit;
2747 EIG_REL_GAP_HOMO_VEC.resize(maxiter + 1);
2748 EIG_REL_GAP_LUMO_VEC.clear();
2749 EIG_REL_GAP_LUMO_VEC.resize(maxiter + 1);
2750 EIG_ABS_GAP_HOMO_VEC.clear();
2751 EIG_ABS_GAP_HOMO_VEC.resize(maxiter + 1);
2752 EIG_ABS_GAP_LUMO_VEC.clear();
2753 EIG_ABS_GAP_LUMO_VEC.resize(maxiter + 1);
2757 real homo0 = homo_bounds_X0.low();
2758 real lumo0 = lumo_bounds_X0.upp();
2762 real homo_map, lumo_map, one_map, zero_map, sigma;
2764 real homo = homo0, lumo = lumo0;
2766 for (
int i = 1; i <= maxiter; ++i)
2768 homo = apply_poly(i, homo);
2769 lumo = apply_poly(i, lumo);
2771 sigma = SIGMA_HOMO_VEC[i];
2773 homo_map = (homo - sigma) * (homo - sigma);
2774 lumo_map = (lumo - sigma) * (lumo - sigma);
2775 one_map = (one - sigma) * (one - sigma);
2776 zero_map = (
zero - sigma) * (
zero - sigma);
2778 EIG_ABS_GAP_HOMO_VEC[i] =
min(lumo_map - homo_map, one_map - homo_map);
2779 EIG_REL_GAP_HOMO_VEC[i] = EIG_ABS_GAP_HOMO_VEC[i] /
2780 max(zero_map - homo_map, one_map - homo_map);
2783 homo = homo0, lumo = lumo0;
2785 for (
int i = 1; i <= maxiter; ++i)
2787 homo = apply_poly(i, homo);
2788 lumo = apply_poly(i, lumo);
2790 sigma = SIGMA_LUMO_VEC[i];
2792 homo_map = (homo - sigma) * (homo - sigma);
2793 lumo_map = (lumo - sigma) * (lumo - sigma);
2794 zero_map = (
zero - sigma) * (
zero - sigma);
2795 one_map = (one - sigma) * (one - sigma);
2797 EIG_ABS_GAP_LUMO_VEC[i] =
min(homo_map - lumo_map, zero_map - lumo_map);
2798 EIG_REL_GAP_LUMO_VEC[i] = EIG_ABS_GAP_LUMO_VEC[i] /
2799 max(zero_map - lumo_map, one_map - lumo_map);
2805 #ifdef USE_CHUNKS_AND_TASKS 2806 template<
typename MatrixType>
2809 throw "compute_eigenvector() is not implemented for CHTMatrix.";
2814 template<
typename MatrixType>
2817 assert(eigVecHOMOorLUMO != NULL);
2818 real acc_eigv = eigensolver_accuracy;
2820 #ifdef DEBUG_PURI_OUTPUT 2825 if (compute_eigenvectors_in_each_iteration && use_prev_vector_as_initial_guess)
2830 *eigVecHOMO = eigVecHOMORef;
2834 *eigVecLUMO = eigVecLUMORef;
2839 vector<real> eigValTmp(number_of_eigenvalues);
2849 vector<VectorType> eigVec(number_of_eigenvalues,
VectorType(
rows));
2850 if (use_prev_vector_as_initial_guess)
2852 use_prev_vector_as_initial_guess = 0;
2855 if (!eigVecHOMO->is_empty())
2857 eigVec[0] = *eigVecHOMO;
2859 use_prev_vector_as_initial_guess = 1;
2864 if (!eigVecLUMO->is_empty())
2866 eigVec[0] = *eigVecLUMO;
2868 use_prev_vector_as_initial_guess = 1;
2878 vector<int> num_iter_solver(number_of_eigenvalues, -1);
2885 eigenvectors_iterative_method_str, num_iter_solver,
2886 eigensolver_maxiter);
2890 if (num_iter_solver.empty())
2892 throw std::runtime_error(
"Error in compute_eigenvector() : (num_iter_solver.empty())");
2897 bool is_homo_tmp =
false, is_lumo_tmp =
false;
2899 if (num_iter_solver[0] == eigensolver_maxiter)
2905 is_homo_tmp = is_homo, is_lumo_tmp = !is_homo;
2907 get_interval_with_lambda(eigVal, eigVec[0], is_homo_tmp, is_lumo_tmp, it);
2910 really_good_iterations_homo.push_back(it);
2911 *eigVecHOMO = eigVec[0];
2912 info.homo_eigenvector_is_computed =
true;
2913 info.homo_eigenvector_is_computed_in_iter = it;
2914 info.homo_eigensolver_iter = num_iter_solver[0];
2915 info.homo_eigensolver_time = eigv_elapsed_wall_sec;
2916 eigValHOMO = eigVal;
2918 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
2920 else if (is_lumo_tmp)
2922 really_good_iterations_lumo.push_back(it);
2923 *eigVecLUMO = eigVec[0];
2924 info.lumo_eigenvector_is_computed =
true;
2925 info.lumo_eigenvector_is_computed_in_iter = it;
2926 info.lumo_eigensolver_iter = num_iter_solver[0];
2927 info.lumo_eigensolver_time = eigv_elapsed_wall_sec;
2928 eigValLUMO = eigVal;
2930 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
2935 ": number of eigensolver iterations is %d", it, num_iter_solver[0]);
2940 if (compute_eigenvectors_in_each_iteration)
2942 save_eigenvectors_to_file(is_homo_tmp, is_lumo_tmp, it);
2947 #endif // USE_CHUNKS_AND_TASKS 2950 template<
typename MatrixType>
2954 #ifndef USE_CHUNKS_AND_TASKS 2960 template<
typename MatrixType>
2964 #ifndef USE_CHUNKS_AND_TASKS 2970 template<
typename MatrixType>
2974 if (is_lumo || is_homo)
2976 std::vector<real> fullVector;
2980 eigVecHOMO->fullvector(fullVector);
2985 eigVecLUMO->fullvector(fullVector);
2993 name << eig_name <<
"_" << it <<
"_scf_step_" << scf_step <<
".txt";
2997 name << eig_name <<
"_" << it <<
".txt";
2999 if (
write_vector(name.str().c_str(), fullVector) == -1)
3001 throw std::runtime_error(
"Error in save_eigenvectors_to_file() : error in write_vector.");
3009 template<
typename MatrixType>
3013 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3015 eigVal = approx_eig;
3019 template<
typename MatrixType>
3022 assert(is_homo || is_lumo);
3024 bool is_homo_init = is_homo;
3025 bool is_lumo_init = is_lumo;
3037 real low_lumo_F_bound, low_homo_F_bound;
3038 real upp_lumo_F_bound, upp_homo_F_bound;
3043 low_lumo_F_bound = lumo_bounds_F.low();
3044 upp_lumo_F_bound = lumo_bounds_F.upp();
3045 low_homo_F_bound = homo_bounds_F.low();
3046 upp_homo_F_bound = homo_bounds_F.upp();
3051 low_lumo_F_bound = info.lumo_estim_low_F;
3052 upp_lumo_F_bound = info.lumo_estim_upp_F;
3053 low_homo_F_bound = info.homo_estim_low_F;
3054 upp_homo_F_bound = info.homo_estim_upp_F;
3058 throw std::runtime_error(
"Error in get_interval_with_lambda() : unexpected eigenvectors_method value.");
3061 #ifdef DEBUG_PURI_OUTPUT 3066 real approx_eig = eigvec::compute_rayleigh_quotient<real>(F, eigVec);
3069 eigVal = approx_eig;
3074 if ((approx_eig <= upp_homo_F_bound + flex_tolerance) && (approx_eig >= low_homo_F_bound - flex_tolerance))
3080 "HOMO bounds are [ %lf , %lf ]", (
double)approx_eig, (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
3082 iter_for_homo = iter;
3084 if (is_lumo_init && (eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3086 iter_for_lumo = iter_for_lumo + 1;
3091 else if ((approx_eig <= upp_lumo_F_bound + flex_tolerance) && (approx_eig >= low_lumo_F_bound - flex_tolerance))
3097 "LUMO interval [ %lf , %lf ]", (
double)approx_eig, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
3099 iter_for_lumo = iter;
3101 if (is_homo_init && (eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3103 iter_for_homo = iter_for_homo + 1;
3114 (
double)approx_eig, (
double)low_homo_F_bound, (
double)upp_homo_F_bound, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
3119 if ((eigenvectors_method ==
EIG_SQUARE_INT) && try_eigv_on_next_iteration_if_fail)
3121 iter_for_homo = iter_for_homo + 1;
3133 template<
typename MatrixType>
3137 f.open(filename, std::ios::out);
3144 int it = info.total_it;
3147 for (
int i = 0; i <= it; ++i)
3149 f << VecPoly[i] <<
" ";
3151 f <<
"];" << std::endl;
3157 for (
int i = 0; i <= it; ++i)
3159 f << (double)info.Iterations[i].XmX2_eucl <<
" ";
3161 f <<
"];" << std::endl;
3162 f <<
" norm_letter = '2';" << std::endl;
3167 for (
int i = 0; i <= it; ++i)
3169 f << (double)info.Iterations[i].XmX2_fro_norm <<
" ";
3171 f <<
"];" << std::endl;
3172 f <<
" norm_letter = 'F';" << std::endl;
3177 for (
int i = 0; i <= it; ++i)
3179 f << (double)info.Iterations[i].XmX2_mixed_norm <<
" ";
3181 f <<
"];" << std::endl;
3182 f <<
" norm_letter = 'M';" << std::endl;
3186 throw "Wrong norm in PurificationGeneral::gen_matlab_file_norm_diff()";
3189 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3190 f <<
"it = " << it <<
";" << std::endl;
3191 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3192 f <<
"fighandle = figure; clf;" << std::endl;
3193 f <<
"MARKER = ['o', '+'];" << std::endl;
3194 f <<
"semilogy(0:stop_iteration, X_norm(1:stop_iteration+1), '-b', plot_props{:});" << std::endl;
3195 f <<
"hold on" << std::endl;
3196 f <<
"for i = 1:stop_iteration+1" << std::endl;
3197 f <<
" if POLY(i) == 1" << std::endl;
3198 f <<
" h1 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3199 f <<
" else" << std::endl;
3200 f <<
" h2 = semilogy(i-1, X_norm(i), [MARKER((POLY(i) == 1) + 1) 'b'], plot_props{:});" << std::endl;
3201 f <<
" end" << std::endl;
3202 f <<
"end" << std::endl;
3203 f <<
"if stop_iteration ~= it" << std::endl;
3204 f <<
"h3 = semilogy(stop_iteration+1:it, X_norm(stop_iteration+2:it+1), '-.vr', plot_props{:});" << std::endl;
3205 f <<
"semilogy(stop_iteration:stop_iteration+1, X_norm(stop_iteration+1:stop_iteration+2), '-.r', plot_props{:});" << std::endl;
3206 f <<
"legend([h1 h2 h3],{'$x^2$', '$2x-x^2$', 'After stop'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3207 f <<
"else" << std::endl;
3208 f <<
"legend([h1 h2],{'$x^2$', '$2x-x^2$'}, 'Interpreter','latex', 'Location','SouthWest');" << std::endl;
3209 f <<
"end" << std::endl;
3210 f <<
"xlabel('Iteration SP2', 'Interpreter','latex');" << std::endl;
3211 f <<
"ylabel({['$\\|X_i-X_i^2\\|_{' norm_letter '}$']},'interpreter','latex');" << std::endl;
3212 f <<
"grid on" << std::endl;
3213 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3214 f <<
"set(gca,'FontSize',20);" << std::endl;
3215 f <<
"xlim([0 it]);" << std::endl;
3216 f <<
"ylim([-inf inf]);" << std::endl;
3217 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3218 f <<
"a = 16; S = logspace(-a, 1, a+2);" << std::endl;
3219 f <<
"set(gca,'YTick',S(1:2:end));" << std::endl;
3221 f <<
"hold off" << std::endl;
3223 f <<
"% print( fighandle, '-depsc2', 'norm_diff.eps');" << std::endl;
3229 template<
typename MatrixType>
3233 f.open(filename, std::ios::out);
3240 int it = info.total_it;
3243 for (
int i = 0; i <= it; ++i)
3245 f << (double)info.Iterations[i].threshold_X <<
" ";
3247 f <<
"];" << std::endl;
3249 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3250 f <<
"it = " << it <<
";" << std::endl;
3251 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3252 f <<
"fighandle = figure; clf;" << std::endl;
3253 f <<
"semilogy(0:stop_iteration, Thresh(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3254 f <<
"hold on" << std::endl;
3255 f <<
"if stop_iteration ~= it" << std::endl;
3256 f <<
"semilogy(stop_iteration+1:it, Thresh(stop_iteration+2:it+1), '-^r', plot_props{:});" << std::endl;
3257 f <<
"semilogy(stop_iteration:stop_iteration+1, Thresh(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3258 f <<
"legend('before stop', 'after stop', 'Location','NorthWest');" << std::endl;
3259 f <<
"end" << std::endl;
3260 f <<
"grid on" << std::endl;
3261 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3262 f <<
"set(gca,'FontSize',20);" << std::endl;
3263 f <<
"xlim([0 it]);" << std::endl;
3264 f <<
"ylim([-inf inf]);" << std::endl;
3265 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3266 f <<
"hold off" << std::endl;
3267 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3268 f <<
"ylabel('Threshold value', 'interpreter','latex');" << std::endl;
3269 f <<
"% print( fighandle, '-depsc2', 'threshold.eps');" << std::endl;
3275 template<
typename MatrixType>
3279 f.open(filename, std::ios::out);
3286 int it = info.total_it;
3289 for (
int i = 0; i <= it; ++i)
3291 f << (double)info.Iterations[i].NNZ_X <<
" ";
3293 f <<
"];" << std::endl;
3297 for (
int i = 0; i <= it; ++i)
3299 f << (double)info.Iterations[i].NNZ_X2 <<
" ";
3301 f <<
"];" << std::endl;
3305 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3306 f <<
"it = " << it <<
";" << std::endl;
3307 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3308 f <<
"fighandle = figure; clf;" << std::endl;
3309 f <<
"h2 = plot(0:stop_iteration, NNZ_X(1:stop_iteration+1), '-ob', plot_props{:});" << std::endl;
3310 f <<
"hold on" << std::endl;
3311 f <<
"h1 = plot(0:stop_iteration, NNZ_X2(1:stop_iteration+1), '-vm', plot_props{:});" << std::endl;
3312 f <<
"if stop_iteration ~= it" << std::endl;
3313 f <<
"plot(stop_iteration+1:it, NNZ_X(stop_iteration+2:it+1), '-vr', plot_props{:});" << std::endl;
3314 f <<
"plot(stop_iteration+1:it, NNZ_X2(stop_iteration+2:it+1), '-*r', plot_props{:});" << std::endl;
3315 f <<
"plot(stop_iteration:stop_iteration+1, NNZ_X(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3316 f <<
"plot(stop_iteration:stop_iteration+1, NNZ_X2(stop_iteration+1:stop_iteration+2), '-r', plot_props{:});" << std::endl;
3317 f <<
"end" << std::endl;
3318 f <<
"legend([h1, h2], {'$X^2$', '$X$'}, 'interpreter','latex');" << std::endl;
3319 f <<
"grid on" << std::endl;
3320 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3321 f <<
"set(gca,'FontSize',20);" << std::endl;
3322 f <<
"xlim([0 it]);" << std::endl;
3323 f <<
"ylim([0 inf]);" << std::endl;
3324 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3325 f <<
"hold off" << std::endl;
3326 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3327 f <<
"ylabel('NNZ [\\%]', 'interpreter','latex');" << std::endl;
3328 f <<
"% print( fighandle, '-depsc2', 'nnz.eps');" << std::endl;
3334 template<
typename MatrixType>
3338 f.open(filename, std::ios::out);
3345 int it = info.total_it;
3348 for (
int i = 0; i <= it; ++i)
3350 f << (double)(1 - info.Iterations[i].homo_bound_upp - info.Iterations[i].lumo_bound_upp) <<
" ";
3352 f <<
"];" << std::endl;
3354 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3355 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3356 f <<
"fighandle = figure; clf;" << std::endl;
3357 f <<
"plot(0:stop_iteration, 1./in(1:stop_iteration+1), '-*r', plot_props{:});" << std::endl;
3358 f <<
"grid on" << std::endl;
3359 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3360 f <<
"set(gca,'FontSize',20);" << std::endl;
3361 f <<
"xlim([0 it]);" << std::endl;
3362 f <<
"ylim([-inf inf]);" << std::endl;
3363 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3364 f <<
"hold off" << std::endl;
3365 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3366 f <<
"ylabel('$\\kappa$', 'interpreter','latex');" << std::endl;
3367 f <<
"% print( fighandle, '-depsc2', 'cond.eps');" << std::endl;
3373 template<
typename MatrixType>
3377 f.open(filename, std::ios::out);
3384 int it = info.total_it;
3387 for (
int i = 0; i <= it; ++i)
3389 f << (double)info.Iterations[i].homo_bound_low <<
" ";
3391 f <<
"];" << std::endl;
3395 for (
int i = 0; i <= it; ++i)
3397 f << (double)info.Iterations[i].homo_bound_upp <<
" ";
3399 f <<
"];" << std::endl;
3403 for (
int i = 0; i <= it; ++i)
3405 f << (double)info.Iterations[i].lumo_bound_low <<
" ";
3407 f <<
"];" << std::endl;
3411 for (
int i = 0; i <= it; ++i)
3413 f << (double)info.Iterations[i].lumo_bound_upp <<
" ";
3415 f <<
"];" << std::endl;
3419 f <<
"stop_iteration = " << it - info.additional_iterations <<
";" << std::endl;
3420 f <<
"plot_props = {'LineWidth', 2, 'MarkerSize', 8};" << std::endl;
3421 f <<
"fighandle = figure; clf;" << std::endl;
3422 f <<
"semilogy(0:stop_iteration, homo_upp(1:stop_iteration+1), '-^b', plot_props{:});" << std::endl;
3423 f <<
"hold on" << std::endl;
3424 f <<
"semilogy(0:stop_iteration, homo_low(1:stop_iteration+1), '-vb', plot_props{:});" << std::endl;
3425 f <<
"semilogy(0:stop_iteration, lumo_low(1:stop_iteration+1), '-vr', plot_props{:});" << std::endl;
3426 f <<
"semilogy(0:stop_iteration, lumo_upp(1:stop_iteration+1), '-^r', plot_props{:});" << std::endl;
3427 f <<
"grid on" << std::endl;
3428 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3429 f <<
"set(gca,'FontSize',20);" << std::endl;
3430 f <<
"xlim([0 stop_iteration]);" << std::endl;
3431 f <<
"set(gca,'XTick',[0 5:5:stop_iteration]);" << std::endl;
3432 f <<
"ylim([-inf inf]);" << std::endl;
3433 f <<
"hold off" << std::endl;
3434 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3435 f <<
"ylabel('Homo and lumo bounds', 'interpreter','latex');" << std::endl;
3436 f <<
"% print( fighandle, '-depsc2', 'eigs.eps');" << std::endl;
3442 template<
typename MatrixType>
3446 f.open(filename, std::ios::out);
3453 int it = info.total_it;
3455 f <<
"time_total = [";
3456 for (
int i = 0; i <= it; ++i)
3458 f << std::setprecision(16) << (double)info.Iterations[i].total_time <<
" ";
3460 f <<
"];" << std::endl;
3462 f <<
"time_square = [";
3463 for (
int i = 0; i <= it; ++i)
3465 f << std::setprecision(16) << (double)info.Iterations[i].Xsquare_time <<
" ";
3467 f <<
"];" << std::endl;
3469 f <<
"time_trunc = [";
3470 for (
int i = 0; i <= it; ++i)
3472 f << std::setprecision(16) << (double)info.Iterations[i].trunc_time <<
" ";
3474 f <<
"];" << std::endl;
3476 if (info.compute_eigenvectors_in_this_SCF_cycle)
3478 f <<
"time_eigenvectors_homo = [";
3479 for (
int i = 0; i <= it; ++i)
3481 f << std::setprecision(16) << (double)info.Iterations[i].orbital_homo_time <<
" ";
3483 f <<
"];" << std::endl;
3484 f <<
"time_eigenvectors_lumo = [";
3485 for (
int i = 0; i <= it; ++i)
3487 f << std::setprecision(16) << (double)info.Iterations[i].orbital_lumo_time <<
" ";
3489 f <<
"];" << std::endl;
3491 f <<
"time_solver_homo = [";
3492 for (
int i = 0; i <= it; ++i)
3494 f << std::setprecision(16) << (double)info.Iterations[i].homo_eig_solver_time <<
" ";
3496 f <<
"];" << std::endl;
3497 f <<
"time_solver_lumo = [";
3498 for (
int i = 0; i <= it; ++i)
3500 f << std::setprecision(16) << (double)info.Iterations[i].lumo_eig_solver_time <<
" ";
3502 f <<
"];" << std::endl;
3508 f <<
"X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_eigenvectors_homo-time_solver_homo;" 3509 " time_eigenvectors_lumo-time_solver_lumo; " 3510 " time_total - time_square - time_trunc - time_eigenvectors_homo - time_eigenvectors_lumo];" << std::endl;
3514 f <<
"time_DX_homo = [";
3515 for (
int i = 0; i <= it; ++i)
3517 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_homo_time <<
" ";
3519 f <<
"];" << std::endl;
3520 f <<
"time_DX_lumo = [";
3521 for (
int i = 0; i <= it; ++i)
3523 f << std::setprecision(16) << (double)info.Iterations[i].DX_mult_lumo_time <<
" ";
3525 f <<
"];" << std::endl;
3529 f <<
"X = [time_square; time_trunc; time_solver_homo; time_solver_lumo; time_DX_homo; time_DX_lumo;" 3530 " time_eigenvectors_homo - time_DX_homo - time_solver_homo; time_eigenvectors_lumo - time_DX_lumo - time_solver_lumo;" 3531 " time_total - time_square - time_trunc];" << std::endl;
3536 f <<
"X = [time_square; time_trunc; time_total - time_square - time_trunc];" << std::endl;
3539 f <<
"it = " << it <<
";" << std::endl;
3540 f <<
"xtick = 0:it;" << std::endl;
3541 f <<
"fighandle = figure; clf;" << std::endl;
3542 f <<
"b=bar(xtick, X', 'stacked');" << std::endl;
3543 f <<
"axis tight;" << std::endl;
3544 f <<
"grid on" << std::endl;
3545 f <<
"set(gca, 'XMinorGrid','off', 'YMinorGrid','off', 'GridAlpha',0.6, 'GridLineStyle', '--')" << std::endl;
3546 f <<
"set(gca,'FontSize',20);" << std::endl;
3547 f <<
"xlim([0 it]);" << std::endl;
3548 f <<
"set(gca,'XTick',[0 5:5:it]);" << std::endl;
3549 f <<
"ylim([-inf inf]);" << std::endl;
3550 f <<
"xlabel('Iteration SP2', 'interpreter','latex');" << std::endl;
3551 f <<
"ylabel('Time [s]', 'interpreter','latex');" << std::endl;
3552 if (info.compute_eigenvectors_in_this_SCF_cycle)
3559 f <<
"legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3563 f <<
"legend(flipud(b(:)), {'other', 'lumo other', 'homo other', 'DX (lumo)', 'DX (homo)', 'lumo solver', 'homo solver', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3568 f <<
"legend(flipud(b(:)), {'other', 'truncation', '$X^2$'}, 'interpreter','latex');" << std::endl;
3570 f <<
"% print( fighandle, '-depsc2', 'time.eps');" << std::endl;
3581 template<
typename MatrixType>
3585 f.open(filename, std::ios::out);
3592 f <<
"import numpy as np" << std::endl;
3593 f <<
"import pylab as pl" << std::endl;
3594 f <<
"import matplotlib.font_manager as font_manager" << std::endl;
3595 f <<
"from matplotlib import rc" << std::endl;
3596 f <<
"rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})" << std::endl;
3597 f <<
"rc('text', usetex=True)" << std::endl;
3600 int it = info.total_it;
3601 f <<
"NNZ_X = np.array([";
3603 for (
int i = 0; i <= it; ++i)
3605 f << (double)info.Iterations[i].NNZ_X <<
", ";
3607 f <<
"])" << std::endl;
3609 f <<
"NNZ_X2 = np.array([";
3611 for (
int i = 0; i <= it; ++i)
3613 f << (double)info.Iterations[i].NNZ_X2 <<
", ";
3615 f <<
"])" << std::endl;
3617 f <<
"stop_iteration = " << it - info.additional_iterations << std::endl;
3618 f <<
"it = " << it << std::endl;
3619 f <<
"font_prop = font_manager.FontProperties(size=20)" << std::endl;
3620 f <<
"prop = {'markersize':8, 'fillstyle':'none', 'linewidth':2, 'markeredgewidth':2}" << std::endl;
3621 f <<
"fig1 = pl.figure(figsize = (8, 6), num='nnz')" << std::endl;
3622 f <<
"p1, = pl.plot(range(0, stop_iteration+1), NNZ_X2[0:stop_iteration+1], '-vm', **prop);" << std::endl;
3623 f <<
"p2, = pl.plot(range(0, stop_iteration+1), NNZ_X[0:stop_iteration+1], '-ob', **prop);" << std::endl;
3624 f <<
"if stop_iteration != it:" << std::endl;
3625 f <<
" pl.plot(range(stop_iteration+1,it+1), NNZ_X[stop_iteration+1:it+1], '-vr', **prop);" << std::endl;
3626 f <<
" pl.plot(range(stop_iteration+1,it+1), NNZ_X2[stop_iteration+1:it+1], '-*r', **prop);" << std::endl;
3627 f <<
" pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
3628 f <<
" pl.plot(range(stop_iteration,stop_iteration+2), NNZ_X2[stop_iteration:stop_iteration+2], '-r', **prop);" << std::endl;
3630 f <<
"pl.xlim(0, it);" << std::endl;
3631 f <<
"pl.ylim(ymin=0);" << std::endl;
3632 f <<
"pl.xlabel('Iteration SP2', fontproperties=font_prop);" << std::endl;
3633 f <<
"pl.ylabel('NNZ [\\%]', fontproperties=font_prop); " << std::endl;
3634 f <<
"pl.legend((p1, p2), ('$X^2$', '$X$'), loc='best', numpoints=1)" << std::endl;
3635 f <<
"pl.xticks(range(5, it, 5))" << std::endl;
3636 f <<
"pl.tick_params(labelsize=20)" << std::endl;
3637 f <<
"pl.grid(which='major', alpha=0.5, color='k', linestyle='--', linewidth=1) " << std::endl;
3638 f <<
"pl.savefig('nnz.eps', format='eps', dpi=1000)" << std::endl;
3639 f <<
"pl.show()" << std::endl;
3645 #endif //HEADER_PURIFICATION_GENERAL int EIG_LANCZOS_INT
Definition: purification_general.cc:57
real error_sub
Allowed error in invariant subspaces.
Definition: purification_general.h:419
real error_per_it
Error allowed in each iteration due to truncation.
Definition: purification_general.h:422
Treal upp() const
Definition: Interval.h:145
void set_compute_eigenvectors_in_each_iteration()
Definition: purification_general.h:239
VectorTypeInt really_good_iterations_homo
Iterations where homo eigenvector is actually computed.
Definition: purification_general.h:500
real total_time
Definition: puri_info.h:61
real lumo_bound_upp
Definition: puri_info.h:87
int eigensolver_maxiter
Definition: purification_general.h:474
VectorTypeReal EIG_ABS_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for lumo eigenvalue.
Definition: purification_general.h:441
#define LOG_CAT_WARNING
Definition: output.h:48
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)=0
Treal template_blas_pow(Treal x, Treal y)
PurificationGeneral()
Definition: purification_general.h:122
virtual bool is_initialized() const
Check is function initialize() is already called.
Definition: purification_general.h:152
double ergo_real
Definition: realtype.h:69
VectorTypeReal ITER_ERROR_VEC
(Eigenvectors) Maximum error introduced in each iteration.
Definition: purification_general.h:439
ergo_real real
Definition: purification_general.h:114
IntervalType spectrum_bounds
Outer bounds for the whole spectrum of F/Xi.
Definition: purification_general.h:462
mat::SizesAndBlocks rows
Definition: test.cc:51
bool computed_spectrum_bounds
Definition: purification_general.h:463
Functionality for writing output messages to a text file.
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition: purification_general.h:436
int use_new_stopping_criterion
True for new stopping criterion.
Definition: purification_general.h:401
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:237
mat::SizesAndBlocks cols
Definition: test.cc:52
virtual void set_init_params()=0
real homo_bound_upp
Definition: puri_info.h:85
real lumo_bound_low
Definition: puri_info.h:86
void enable_printf_output()
Definition: output.cc:66
void output_norms_and_traces(IterationInfo &iter_info) const
Definition: purification_general.h:785
intervalType IntervalType
Definition: random_matrices.h:67
Defined namespace eigvec containing functions for computing largest eigenvalues and corresponding eig...
VectorTypeReal EIG_ABS_GAP_LUMO_VEC
Definition: purification_general.h:441
Wrapper routines for different parts of the integral code, including conversion of matrices from/to t...
Treal low() const
Definition: Interval.h:144
int write_vector(const char *filename, const vector< real > &v)
Definition: files_dense.cc:206
virtual void stopping_criterion(IterationInfo &iter_info, int &stop, real &estim_order)
Choose stopping criterion and check it.
Definition: purification_general.h:1719
int use_prev_vector_as_initial_guess
Definition: purification_general.h:479
ergo_real real
Definition: purification_general.h:64
Basic OS access utilities.
virtual void purify_bounds(const int it)=0
IntervalType lumo_bounds_X0
Initial lumo bounds for X.
Definition: purification_general.h:453
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:111
real eigensolver_accuracy
Definition: purification_general.h:472
MatrixType X_homo
Definition: purification_general.h:257
real TOL_OVERLAPPING_BOUNDS
If the difference between inner bounds for homo and lumo is less then this tolerance, i.e.
Definition: purification_general.cc:45
void get_eigenvalue_of_F_from_eigv_of_Xi(real &eigVal, const VectorType &eigVec)
Definition: purification_general.h:3010
static real get_epsilon()
Get machine epsilon.
Definition: purification_general.h:198
void find_truncation_thresh_every_iter()
Definition: purification_general.h:2690
VectorTypeInt good_iterations_homo
Iterations where homo eigenvector can be computed.
Definition: purification_general.h:497
mat::Interval< ergo_real > intervalType
Definition: matrix_typedefs.h:78
Proxy structs used by the matrix API.
void get_interval_with_lambda(real &eigVal, VectorType &eigVec, bool &is_homo, bool &is_lumo, const int iter)
Definition: purification_general.h:3020
virtual void check_standard_stopping_criterion(const real XmX2_norm, int &stop)
Check stopping criterion (obsolete).
Definition: purification_general.h:1790
VectorTypeInt good_iterations_lumo
Iterations where homo eigenvector can be computed.
Definition: purification_general.h:498
string eigenvectors_iterative_method_str
Definition: purification_general.h:477
File contataining definitions of constants used in the recursive expansion.
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
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:2628
IntervalType homo_bounds
(1-homo) bounds for Xi in iteration i
Definition: purification_general.h:449
IntervalType lumo_bounds
Lumo bounds for Xi in iteration i.
Definition: purification_general.h:450
int get_int_eig_iter_method(string eigenvectors_iterative_method)
Definition: purification_general.h:766
File containing declaration of functions for reading/writing dense matrices and vectors.
real constant_C
Asymptotic constant C needed for the new stopping criterion.
Definition: purification_general.h:424
real homo_bound_low
Definition: puri_info.h:84
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:405
real XmX2_mixed_norm
Definition: puri_info.h:75
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
generalVector VectorType
Definition: purification_general.h:119
MatrixType X_lumo
Save matrix Xi in certain iterations for computing eigenvectors (projection method).
Definition: purification_general.h:257
#define LOG_AREA_DENSFROMF
Definition: output.h:61
virtual ~PurificationGeneral()
Definition: purification_general.h:171
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
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:3230
virtual void initialize(const MatrixType &F_, const IntervalType &lumo_bounds_, const IntervalType &homo_bounds_, int maxit_, real error_sub_, real error_eig_, int use_new_stopping_criterion_, NormType norm_truncation, NormType norm_stop_crit, int nocc_)
Set imporatant parameters for the recursive expansion.
Definition: purification_general.h:521
MatrixType MatrixTypeWrapper
Definition: purification_general.h:120
std::vector< int > VectorTypeInt
Definition: purification_general.h:117
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
File containing declaration of functions for reading/writing sparse matrices from/to binary files...
int get_int_eig_method(string eigenvectors_method)
Definition: purification_general.h:747
virtual void PurificationStart()
Start recursive expansion.
Definition: purification_general.h:825
VectorTypeReal EIG_REL_GAP_LUMO_VEC
Definition: purification_general.h:442
real trunc_time
Definition: puri_info.h:59
int EIG_PROJECTION_INT
Definition: purification_general.cc:55
void save_matrix_now(string str)
Definition: purification_general.h:2070
generalVector VectorType
Definition: GetDensFromFock.cc:62
#define ORDER
Order of convergence of recursive expansion.
Definition: constants.h:41
Treal template_blas_fabs(Treal x)
bool compute_eigenvectors_in_this_SCF_cycle
Definition: purification_general.h:481
static real get_max_double()
Get largest number.
Definition: purification_general.h:202
ergo_real real
Definition: test.cc:46
int additional_iterations
Do a few more iterations after convergence.
Definition: purification_general.h:402
void compute_spectrum_bounds()
Compute spectrum bounds.
Definition: purification_general.h:1630
File containing declarations of functions for reading/writing sparse matrices from/to mtx (MatrixMark...
void set_spectrum_bounds(real eigmin, real eigmax)
Set spectrum bounds.
Definition: purification_general.h:1606
int iter_for_homo
Definition: purification_general.h:494
real frob_diff_time
Definition: puri_info.h:66
virtual void clear()
Clear all matrices in the class.
Definition: purification_general.h:177
VectorTypeReal SIGMA_HOMO_VEC
Definition: purification_general.h:440
void unset_compute_eigenvectors_in_each_iteration()
Definition: purification_general.h:240
#define UNIT_one_eV
Definition: units.h:45
#define max(a, b)
Definition: integrator.cc:87
int nocc
Number of occupied orbitals.
Definition: purification_general.h:407
VectorType * eigVecHOMO
Definition: purification_general.h:485
intervalType IntervalType
Definition: purification_general.h:115
real Xsquare_time
Definition: puri_info.h:58
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:3374
void discard_homo_eigenvector()
Definition: purification_general.h:1353
Treal template_blas_log(Treal x)
Definition of the main floating-point datatype used; the ergo_real type.
#define LOG_CAT_INFO
Definition: output.h:49
Definition: matInclude.h:139
NormType normPuriStopCrit
Norm used in the stopping criterion Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm.
Definition: purification_general.h:415
virtual bool puri_is_prepared() const
Check is function prepare_to_purification() is already called.
Definition: purification_general.h:156
double get_nnz_X(size_t &nnzX)
Get nnz of X in %.
Definition: purification_general.h:306
Time-measuring class.
Definition: utilities.h:80
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
real NNZ_X
Definition: puri_info.h:80
virtual void prepare_to_purification()
Prepare data for recursive expansion.
Definition: purification_general.h:863
void gen_python_file_nnz(const char *filename) const
Create PYTHON .py file which plots number of non-zero elements in matrices X_i and X_i^2 in each recu...
Definition: purification_general.h:3582
int min(int a, int b)
Definition: lin_trans.cc:66
int EIG_EMPTY
Definition: purification_general.cc:53
void compute_eigenvector(MatrixType const &M, VectorType *eigVecHOMOorLUMO, int it, bool is_homo)
Definition: purification_general.h:2815
int iter_for_lumo
Definition: purification_general.h:495
real orbital_homo_time
Definition: puri_info.h:67
VectorTypeReal SIGMA_LUMO_VEC
(Eigenvectors) Approximation of shifts in each iteration.
Definition: purification_general.h:440
static real get_min_double()
Get smallest number.
Definition: purification_general.h:206
void print(int area, const char *routine)
Definition: utilities.h:111
virtual real apply_poly(const int it, real x)=0
IntervalType lumo_bounds_F
Initial homo bounds for F.
Definition: purification_general.h:456
bool compute_eigenvectors_in_each_iteration
Compute homo and lumo eigenpairs in every iteration and save eigenvectors in txt files.
Definition: purification_general.h:505
real purify_time
Definition: puri_info.h:60
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:3134
void compute_eigenvectors_without_diagonalization_last_iter_proj()
Definition: purification_general.h:2399
MatrixType Xsq
Matrix X^2.
Definition: purification_general.h:250
void find_eig_gaps_every_iter()
Definition: purification_general.h:2743
virtual real total_subspace_error(int it)
Definition: purification_general.h:1871
Definition: matInclude.h:139
void map_bounds_to_0_1()
Get eigenvalue bounds for X0.
Definition: purification_general.h:1519
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:3276
double get_nnz_X()
Get nnz of X in %.
Definition: purification_general.h:310
virtual void purification_process()
Run recursive expansion.
Definition: purification_general.h:913
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
int eigenvectors_method
Chosen method for computing eigenvectors.
Definition: purification_general.h:467
int EIG_POWER_INT
Definition: purification_general.cc:56
int it
Definition: puri_info.h:56
void writeToTmpFile(MatrixType &A) const
Definition: purification_general.h:2952
VectorTypeReal EIG_REL_GAP_HOMO_VEC
(Eigenvectors) Absolute and relative gap in filter for homo eigenvalue.
Definition: purification_general.h:442
VectorType * eigVecLUMO
Definition: purification_general.h:484
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:2204
virtual void purify_X(const int it)=0
Constants for conversion between units for some common units like Angstrom, electron-volt (eV)...
real eigValHOMO
Definition: purification_general.h:492
real constantC
Definition: puri_info.h:93
VectorType eigVecHOMORef
Definition: purification_general.h:488
Definition: matInclude.h:138
int EIG_SQUARE_INT
Definition: purification_general.cc:54
double get_elapsed_wall_seconds()
Definition: utilities.h:107
MatrixType X
Matrix X.
Definition: purification_general.h:248
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
int maxit
Maximum number of iterations.
Definition: purification_general.h:404
real eucl_acc
Tolerance used for computation of spectral norm.
Definition: purification_general.cc:42
virtual void truncate_matrix(real &thresh, int it)=0
bool puri_is_prepared_flag
Definition: purification_general.h:399
Header file with typedefs for matrix and vector types.
real mixed_acc
Tolerance used for computation of mixed norm.
Definition: purification_general.cc:43
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:1810
Definition: puri_info.h:136
void check_eigenvectors_at_the_end()
Definition: purification_general.h:1382
real mixed_diff_time
Definition: puri_info.h:65
NormType normPuriTrunc
Norm used for the truncation of matrices.
Definition: purification_general.h:411
IntervalType lumo_bounds_F_new
Definition: purification_general.h:459
bool try_eigv_on_next_iteration_if_fail
Definition: purification_general.h:482
virtual void estimate_number_of_iterations(int &estim_num_iter)=0
IntervalType homo_bounds_F
Initial lumo bounds for F.
Definition: purification_general.h:455
double get_nnz_Xsq()
Get nnz of X^2 in %.
Definition: purification_general.h:318
real homo_eig_solver_time
Definition: puri_info.h:71
Definition: puri_info.h:52
int scf_step
Definition: purification_general.h:503
void compute_X()
Get matrix X0 by mapping spectrum of F into [0,1] in reverse order.
Definition: purification_general.h:1500
real lumo_eig_solver_time
Definition: puri_info.h:72
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:1992
real eigValLUMO
Definition: purification_general.h:491
VectorTypeInt really_good_iterations_lumo
Iterations where lumo eigenvector is actually computed.
Definition: purification_general.h:501
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:246
virtual void apply_poly_to_eigs(const int it, real &homo, real &lumo)=0
void get_frob_norm_est(const VectorTypeReal &XmX2_norm_frob, const std::vector< real > &h_in, const std::vector< real > &l_in, VectorTypeReal &YmY2_norm_frob_est)
Definition: purification_general.h:1964
void gen_matlab_file_cond_num(const char *filename) const
Create MATLAB .m file which plots a condition number of a problem of computing the density matrix in ...
Definition: purification_general.h:3335
int get_exact_number_of_puri_iterations()
Definition: purification_general.h:1899
real XmX2_eucl
Definition: puri_info.h:76
real XmX2_fro_norm
Definition: puri_info.h:74
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:3443
Definition: matInclude.h:139
real trace_diff_time
Definition: puri_info.h:64
VectorTypeInt VecPoly
Polynomials computed in the function estimated_number_of_iterations() VecPoly[i] = 1 if we use X=X^2 ...
Definition: purification_general.h:431
double get_nnz_Xsq(size_t &nnzXsq)
Get nnz of X^2 in %.
Definition: purification_general.h:314
mat::normType NormType
Definition: purification_general.h:116
virtual real compute_derivative(const int it, real x, real &DDf)=0
real order
Definition: puri_info.h:77
VectorType eigVecLUMORef
Definition: purification_general.h:487
real NNZ_X2
Definition: puri_info.h:81
void get_spectrum_bounds(real &eigmin, real &eigmax)
Get spectrum bounds.
Definition: purification_general.h:1616
real eucl_diff_time
Definition: puri_info.h:63
real XmX2_trace
Definition: puri_info.h:73
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
virtual void return_constant_C(const int it, real &Cval)=0
int number_of_eigenvalues
Definition: purification_general.h:447
void get_eigenvalue_estimates(const VectorTypeReal &XmX2_norm_mixed, const VectorTypeReal &XmX2_norm_frob, const VectorTypeReal &XmX2_trace)
Definition: purification_general.h:1940
void set_truncation_parameters()
Definition: purification_general.h:1578
static std::string writeAndReadAll()
Definition: FileWritable.cc:509
void set_eigenvectors_params(string eigenvectors_method_, string eigenvectors_iterative_method_, real eigensolver_accuracy_, int eigensolver_maxiter_, int scf_step_, int use_prev_vector_as_initial_guess_, int try_eigv_on_next_iteration_if_fail_, VectorType *eigVecLUMO_, VectorType *eigVecHOMO_)
Set parameters for computing eigenvectors.
Definition: purification_general.h:649
void output_current_memory_usage(int logArea, const char *contextString)
Definition: output.cc:186
virtual void save_other_iter_info(IterationInfo &iter_info, int it)=0
real gap
Definition: puri_info.h:79
MatrixType F
Matrix F.
Definition: purification_general.h:255
int get_est_number_of_puri_iterations()
Definition: purification_general.h:1892
File containing classes IterationInfo and PuriInfo.
void discard_lumo_eigenvector()
Definition: purification_general.h:1369
real orbital_lumo_time
Definition: puri_info.h:68
#define NUM_ADDITIONAL_ITERATIONS
Definition: purification_general.h:76
void compute_eigenvectors_without_diagonalization_on_F(const MatrixType &F, int eigensolver_maxiter_for_F)
Definition: purification_general.h:2139
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:2714
string eigenvectors_method_str
Definition: purification_general.h:476
bool initialized_flag
Definition: purification_general.h:398
Treal template_blas_sqrt(Treal x)
IntervalType homo_bounds_F_new
Definition: purification_general.h:458
real error_eig
Error in eigenvalues (used just in old stopping criterion).
Definition: purification_general.h:420
int eigenvectors_iterative_method
Chosen eigensolver.
Definition: purification_general.h:468
void determine_iteration_for_eigenvectors()
Definition: purification_general.h:2599
normType
Definition: matInclude.h:139
virtual void eigenvalue_bounds_estimation()
Estimate eigenvalues near homo-lumo gap.
Definition: purification_general.h:1316
real stopping_criterion_time
Definition: puri_info.h:62
void save_eigenvectors_to_file(bool is_homo, bool is_lumo, int it)
Definition: purification_general.h:2972
void readFromTmpFile(MatrixType &A) const
Definition: purification_general.h:2962
real threshold_X
Definition: puri_info.h:57
IntervalType homo_bounds_X0
Initial lumo bounds for X.
Definition: purification_general.h:452
void propagate_values_in_each_iter(real value_unocc, real value_occ, VectorTypeReal &out_unocc, VectorTypeReal &out_occ, int nmax)
Definition: purification_general.h:1914