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>
138 const IntervalType& lumo_bounds_,
139 const IntervalType& homo_bounds_,
143 int use_new_stopping_criterion_,
144 NormType norm_truncation,
145 NormType norm_stop_crit,
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_,
235 VectorType *eigVecLUMO_,
236 VectorType *eigVecHOMO_
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); }
326 const VectorTypeReal& XmX2_norm_frob,
327 const VectorTypeReal& XmX2_trace);
331 const std::vector<real>& h_in,
332 const std::vector<real>& l_in,
333 VectorTypeReal& YmY2_norm_frob_est);
338 const VectorTypeReal& XmX2_norm_frob,
339 const VectorTypeReal& XmX2_trace);
343 VectorTypeReal& out_unocc,
344 VectorTypeReal& out_occ,
351 int& chosen_iter_homo);
383 virtual void purify_X(
const int it) = 0;
390 virtual real
apply_poly(
const int it, real x) = 0;
520 template<
typename MatrixType>
527 int use_new_stopping_criterion_,
551 #ifdef ENABLE_PRINTF_OUTPUT 558 " , nocc = %d , NNZ = %lu <-> %.5lf %%",
559 X.get_nrows(),
nocc, nnzX, nnzXp);
578 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
582 #ifdef USE_CHUNKS_AND_TASKS 606 throw std::runtime_error(
"Unknown norm in PurificationGeneral");
610 #ifdef USE_CHUNKS_AND_TASKS 639 #ifdef DEBUG_PURI_OUTPUT 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_,
664 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot compute more than 1 eigenpair.");
693 "Eigenvectors will not be computed in this SCF cycle. Set eigenvectors to NULL.");
714 #ifndef USE_CHUNKS_AND_TASKS 723 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for LUMO!");
735 throw std::runtime_error(
"Error in set_eigenvectors_params() : cannot save initial guess for HOMO!");
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>
831 #ifdef DEBUG_PURI_OUTPUT 862 template<
typename MatrixType>
867 throw std::runtime_error(
"Error in prepare_to_purification() : function is called for an uninitialized class.");
912 template<
typename MatrixType>
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;
949 #ifdef SAVE_MATRIX_IN_EACH_ITERATION 959 #ifdef PURI_OUTPUT_NNZ 969 #ifdef PURI_OUTPUT_NNZ 970 if((
double)thresh >= 0)
981 if((
double)thresh >= 0)
992 #ifdef PURI_OUTPUT_NNZ 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 1037 #ifdef DEBUG_PURI_OUTPUT 1047 ostringstream str_out;
1048 str_out <<
"Iteration " << it;
1050 total_time_stop = total_time.get_elapsed_wall_seconds();
1066 stopping_criterion_time_stop = 0;
1088 while (it <= maxit_tmp)
1102 #ifdef SAVE_MATRIX_IN_EACH_ITERATION 1110 #ifdef PURI_OUTPUT_NNZ 1120 #ifdef PURI_OUTPUT_NNZ 1134 #ifdef PURI_OUTPUT_NNZ 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 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 1234 if ((already_stop_before == 0) && ((stop == 1) || (it ==
maxit)))
1247 assert(maxit_tmp <= (
int)
VecPoly.size());
1249 already_stop_before = 1;
1252 ostringstream str_out;
1253 str_out <<
"Iteration " << it;
1274 iter_info.
order = estim_order;
1303 #ifdef DEBUG_PURI_OUTPUT 1304 if (acc_err_sub != -1)
1315 template<
typename MatrixType>
1352 template<
typename MatrixType>
1368 template<
typename MatrixType>
1381 template<
typename MatrixType>
1410 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, " 1428 "Eigenvalues of the matrix X in this iteration probably already reached the level of numerical errors, " 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);
1455 " discard computed HOMO eigenvector.",
1456 (
double)low_homo_F_bound, (
double)upp_homo_F_bound);
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);
1477 " discard computed LUMO eigenvector.",
1478 (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
1499 template<
typename MatrixType>
1507 #ifdef DEBUG_PURI_OUTPUT 1513 X.add_identity(-eigmax);
1514 X *= ((
real)1.0 / (eigmin - eigmax));
1518 template<
typename MatrixType>
1526 #ifdef DEBUG_PURI_OUTPUT 1541 #ifdef DEBUG_PURI_OUTPUT 1565 lumo_bounds = (lumo_bounds - eigmax) / (eigmin - eigmax);
1577 template<
typename MatrixType>
1580 int estim_num_iter = 0;
1587 if (estim_num_iter <
maxit)
1590 maxit = estim_num_iter;
1605 template<
typename MatrixType>
1615 template<
typename MatrixType>
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)
1718 template<
typename MatrixType>
1721 int it = iter_info.
it;
1722 real XmX2_norm_it = -1, XmX2_norm_itm2 = -1;
1762 throw std::runtime_error(
"Error in stopping_criterion() : unknown matrix norm.");
1789 template<
typename MatrixType>
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)
1836 (XmX2_norm_itm2 < 1) &&
1853 #ifdef DEBUG_PURI_OUTPUT 1869 template<
typename MatrixType>
1873 assert(it <= (
int)
VecGap.size());
1877 for (
int i = 0; i <= it; ++i)
1884 error += normE / (
VecGap[i] - normE);
1891 template<
typename MatrixType>
1898 template<
typename MatrixType>
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];
1934 out_unocc[i] = unocc;
1939 template<
typename MatrixType>
1963 template<
typename MatrixType>
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>
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];
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]);
2061 maxeig * (1 - final_bounds[0]) + mineig * final_bounds[0]);
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>
2143 real dist = (end_shift - start_shift) / 15;
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);
2167 vector<VectorType> eigVec(1,
VectorType(rows));
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);
2190 vector<VectorType> eigVec(1,
VectorType(rows));
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;
2224 compute_homo_now =
true;
2237 compute_lumo_now =
true;
2248 compute_homo_now =
true;
2252 compute_lumo_now =
true;
2263 compute_homo_now =
true;
2270 compute_lumo_now =
true;
2289 M.add_identity(sigma * sigma);
2290 M += ((
real) - 2 * sigma) *
X;
2291 M = ((
real) - 1.0) * M;
2314 name <<
"homo_" << it;
2347 M.add_identity(sigma * sigma);
2348 M += ((
real) - 2 * sigma) *
X;
2349 M = ((
real) - 1.0) * M;
2371 name <<
"lumo_" << it;
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;
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);
2463 DXi = (
real)(-1) * DXi;
2476 info.
Iterations[it].homo_eig_solver_time = homo_solver_time_stop;
2477 info.
Iterations[it].lumo_eig_solver_time = lumo_solver_time_stop;
2505 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP,
X_homo, 0, DXi);
2536 DXi = (
real)(-1) * DXi;
2570 MatrixType::ssmmUpperTriangleOnly((
real)1.0, TMP,
X_lumo, 0, DXi);
2576 DXi = (
real)(-1) * DXi;
2594 #endif // USE_CHUNKS_AND_TASKS 2598 template<
typename MatrixType>
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;
2636 real homoi = homo0, lumoi = lumo0;
2638 real Dh_homo, Dh_lumo, Dgh_homo, Dgh_lumo,
2641 chosen_iter_lumo = -1;
2642 chosen_iter_homo = -1;
2649 for (
int i = 1; i <= maxiter; ++i)
2664 if (Dgh_homo >= Dgh_homo_max)
2666 Dgh_homo_max = Dgh_homo;
2667 chosen_iter_homo = 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;
2703 for (
int i = 1; i <= maxiter; ++i)
2713 template<
typename MatrixType>
2716 int maxiter =
maxit;
2728 for (
int i = 1; i <= maxiter; ++i)
2742 template<
typename MatrixType>
2745 int maxiter =
maxit;
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)
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);
2780 max(zero_map - homo_map, one_map - homo_map);
2783 homo = homo0, lumo = lumo0;
2785 for (
int i = 1; i <= maxiter; ++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);
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);
2820 #ifdef DEBUG_PURI_OUTPUT 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;
2905 is_homo_tmp = is_homo, is_lumo_tmp = !is_homo;
2918 ": %d iterations, %lf wall sec", it, num_iter_solver[0], eigv_elapsed_wall_sec);
2920 else if (is_lumo_tmp)
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]);
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;
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;
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);
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);
3114 (
double)approx_eig, (
double)low_homo_F_bound, (
double)upp_homo_F_bound, (
double)low_lumo_F_bound, (
double)upp_lumo_F_bound);
3133 template<
typename MatrixType>
3137 f.open(filename, std::ios::out);
3147 for (
int i = 0; i <= it; ++i)
3151 f <<
"];" << std::endl;
3157 for (
int i = 0; i <= it; ++i)
3161 f <<
"];" << std::endl;
3162 f <<
" norm_letter = '2';" << std::endl;
3167 for (
int i = 0; i <= it; ++i)
3171 f <<
"];" << std::endl;
3172 f <<
" norm_letter = 'F';" << std::endl;
3177 for (
int i = 0; i <= it; ++i)
3181 f <<
"];" << std::endl;
3182 f <<
" norm_letter = 'M';" << std::endl;
3186 throw "Wrong norm in PurificationGeneral::gen_matlab_file_norm_diff()";
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);
3243 for (
int i = 0; i <= it; ++i)
3247 f <<
"];" << 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);
3289 for (
int i = 0; i <= it; ++i)
3293 f <<
"];" << std::endl;
3297 for (
int i = 0; i <= it; ++i)
3301 f <<
"];" << 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);
3348 for (
int i = 0; i <= it; ++i)
3352 f <<
"];" << 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);
3387 for (
int i = 0; i <= it; ++i)
3391 f <<
"];" << std::endl;
3395 for (
int i = 0; i <= it; ++i)
3399 f <<
"];" << std::endl;
3403 for (
int i = 0; i <= it; ++i)
3407 f <<
"];" << std::endl;
3411 for (
int i = 0; i <= it; ++i)
3415 f <<
"];" << 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);
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;
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;
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;
3601 f <<
"NNZ_X = np.array([";
3603 for (
int i = 0; i <= it; ++i)
3607 f <<
"])" << std::endl;
3609 f <<
"NNZ_X2 = np.array([";
3611 for (
int i = 0; i <= it; ++i)
3615 f <<
"])" << 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 accumulated_error_subspace
Definition: puri_info.h:184
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
real lumo_estim_low_F
Definition: puri_info.h:189
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)
void set_spectrum_bounds(const real lower_spectrum_bound_, const real upper_spectrum_bound_)
Definition: puri_info.cc:92
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.
real homo_estim_upp_F
Definition: puri_info.h:186
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition: purification_general.h:436
void get_vec_traces(std::vector< real > &traces)
Definition: puri_info.cc:76
int use_new_stopping_criterion
True for new stopping criterion.
Definition: purification_general.h:401
int estim_total_it
Definition: puri_info.h:178
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
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
int debug_output
Definition: puri_info.h:204
ergo_real real
Definition: purification_general.h:64
Basic OS access utilities.
virtual void purify_bounds(const int it)=0
real total_time
Definition: puri_info.h:175
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
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:53
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
bool homo_eigenvector_is_computed
Definition: puri_info.h:192
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
int lumo_eigenvector_is_computed_in_iter
Definition: puri_info.h:195
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.
int homo_eigensolver_iter
Definition: puri_info.h:196
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
int stopping_criterion
Definition: puri_info.h:173
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
real eigValLUMO
Definition: puri_info.h:201
void fullvector(std::vector< Treal > &fullVector) const
Definition: VectorGeneral.h:82
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
double homo_eigensolver_time
Definition: puri_info.h:198
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
std::vector< IterationInfo > Iterations
Definition: puri_info.h:203
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
void clear_structure()
Definition: VectorGeneral.h:60
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
#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
void get_vec_mixed_norms(std::vector< real > &norms)
Definition: puri_info.cc:66
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
int total_it
Definition: puri_info.h:176
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
real error_subspace
Definition: puri_info.h:183
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
bool is_empty() const
Definition: VectorGeneral.h:56
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
void resetSizesAndBlocks(SizesAndBlocks const &newRows)
Definition: VectorGeneral.h:51
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
real eigValHOMO
Definition: puri_info.h:200
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
int homo_eigenvector_is_computed_in_iter
Definition: puri_info.h:194
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
int additional_iterations
Definition: puri_info.h:179
real XmX2_eucl
Definition: puri_info.h:76
real XmX2_fro_norm
Definition: puri_info.h:74
int converged
Definition: puri_info.h:181
bool compute_eigenvectors_in_this_SCF_cycle
Definition: puri_info.h:191
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 lumo_estim_upp_F
Definition: puri_info.h:188
real eucl_diff_time
Definition: puri_info.h:63
bool empty() const
Definition: Interval.h:51
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_vec_frob_norms(std::vector< real > &norms)
Definition: puri_info.cc:56
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
double lumo_eigensolver_time
Definition: puri_info.h:199
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.
real homo_estim_low_F
Definition: puri_info.h:187
void discard_lumo_eigenvector()
Definition: purification_general.h:1369
real orbital_lumo_time
Definition: puri_info.h:68
real time_spectrum_bounds
Definition: puri_info.h:177
#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
int lumo_eigensolver_iter
Definition: puri_info.h:197
void determine_iteration_for_eigenvectors()
Definition: purification_general.h:2599
bool lumo_eigenvector_is_computed
Definition: puri_info.h:193
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