39 #ifndef HEADER_PURIFICATION_SP2 40 #define HEADER_PURIFICATION_SP2 50 template<
typename MatrixType>
74 void get_poly(
const int it,
int& poly);
93 template<
typename MatrixType>
96 assert((
int)this->VecPoly.size() > it);
99 if (this->VecPoly[it] == -1)
101 real Xtrace = this->X.trace();
102 real Xsqtrace = this->Xsq.trace();
113 this->VecPoly[it] = 1;
117 this->VecPoly[it] = 0;
123 template<
typename MatrixType>
126 assert((
int)this->VecPoly.size() > it);
127 assert(this->VecPoly[it] != -1);
128 poly = this->VecPoly[it];
132 template<
typename MatrixType>
159 this->Xsq *= ((
real) - 1.0);
160 this->X *= (
real)2.0;
161 this->Xsq += this->X;
164 this->Xsq.transfer(this->X);
168 template<
typename MatrixType>
174 real homo_low, homo_upp, lumo_upp, lumo_low;
182 homo_low = 2 * this->homo_bounds.low() - this->homo_bounds.low() * this->homo_bounds.low();
183 homo_upp = 2 * this->homo_bounds.upp() - this->homo_bounds.upp() * this->homo_bounds.upp();
184 lumo_low = this->lumo_bounds.low() * this->lumo_bounds.low();
185 lumo_upp = this->lumo_bounds.upp() * this->lumo_bounds.upp();
192 lumo_low = 2 * this->lumo_bounds.low() - this->lumo_bounds.low() * this->lumo_bounds.low();
193 lumo_upp = 2 * this->lumo_bounds.upp() - this->lumo_bounds.upp() * this->lumo_bounds.upp();
194 homo_low = this->homo_bounds.low() * this->homo_bounds.low();
195 homo_upp = this->homo_bounds.upp() * this->homo_bounds.upp();
201 this->homo_bounds.intersect(zero_one);
202 this->lumo_bounds.intersect(zero_one);
205 if (this->homo_bounds.empty())
209 if (this->lumo_bounds.empty())
222 template<
typename MatrixType>
231 template<
typename MatrixType>
234 real allowed_error = this->error_per_it;
237 if (allowed_error == 0)
242 assert((
int)this->VecGap.size() > it);
249 if (this->VecGap[it] > 0)
254 tau = (allowed_error * this->VecGap[it]) / (1 + allowed_error);
258 tau = allowed_error * 0.01;
262 #ifdef USE_CHUNKS_AND_TASKS 263 threshold = (this->X).thresh_frob(tau);
265 threshold = (this->X).thresh(tau, this->normPuriTrunc);
277 template<
typename MatrixType>
281 int maxit_tmp = this->maxit + this->additional_iterations;
283 real epsilon = this->get_epsilon();
286 this->check_stopping_criterion_iter = 2;
288 int max_size = maxit_tmp + 1 + this->additional_iterations + 2;
290 this->VecPoly.clear();
291 this->VecPoly.resize(max_size, -1);
293 this->VecGap.clear();
294 this->VecGap.resize(max_size, -1);
297 x = this->lumo_bounds.upp();
298 y = this->homo_bounds.upp();
307 estim_num_iter = this->maxit;
316 this->VecPoly[0] = -1;
317 this->VecGap[0] = 1 - x - y;
321 while (it <= maxit_tmp)
324 if ((x > y) || (it % 2 && (x == y)))
328 this->VecPoly[it] = 1;
334 this->VecPoly[it] = 0;
337 this->VecGap[it] = 1 - x - y;
340 if ((estim_num_iter == -1) &&
341 (x - x * x < epsilon) && (y - y * y < epsilon) &&
342 (this->VecPoly[it] != this->VecPoly[it - 1]))
345 maxit_tmp = it + this->additional_iterations;
359 if ((estim_num_iter == -1) && (it == maxit_tmp + 1))
364 estim_num_iter = this->maxit;
368 this->VecPoly.resize(maxit_tmp + 1);
369 this->VecGap.resize(maxit_tmp + 1);
373 for (
int i = 0; i < (int)this->VecPoly.size(); ++i)
384 template<
typename MatrixType>
387 assert((
int)this->VecPoly.size() > it);
388 assert((
int)this->VecGap.size() > it);
390 iter_info.
poly = this->VecPoly[it];
391 iter_info.
gap = this->VecGap[it];
398 template<
typename MatrixType>
402 for (
int i = it; i >= 1; i--)
411 bounds_from_it[2] = bounds_from_it[2] / (1 +
template_blas_sqrt(1 - bounds_from_it[2]));
412 bounds_from_it[3] = bounds_from_it[3] / (1 +
template_blas_sqrt(1 - bounds_from_it[3]));
416 bounds_from_it[0] = bounds_from_it[0] / (1 +
template_blas_sqrt(1 - bounds_from_it[0]));
417 bounds_from_it[1] = bounds_from_it[1] / (1 +
template_blas_sqrt(1 - bounds_from_it[1]));
451 template<
typename MatrixType>
478 template<
typename MatrixType>
492 homo = 2 * homo - homo * homo;
498 lumo = 2 * lumo - lumo * lumo;
503 template<
typename MatrixType>
517 for (
int i = 1; i <= it; i++)
523 DDf = 2 * Df * Df + (2 * a) * DDf;
529 DDf = -2 * Df * Df + (2 - 2 * a) * DDf;
541 #endif //HEADER_PURIFICATION_SP2 void purify_X(const int it)
Definition: purification_sp2.h:133
ergo_real real
Definition: purification_general.h:114
#define C_SP2
Constant used on the stopping criterion for the SP2 method.
Definition: constants.h:42
intervalType IntervalType
Definition: random_matrices.h:67
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:111
real apply_poly(const int it, real x)
Definition: purification_sp2.h:453
void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2.h:385
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
int poly
Definition: puri_info.h:78
Purification_sp2()
Definition: purification_sp2.h:64
#define LOG_AREA_DENSFROMF
Definition: output.h:61
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2.h:59
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2.h:57
void get_poly(const int it, int &poly)
Definition: purification_sp2.h:124
std::vector< int > VectorTypeInt
Definition: purification_general.h:117
void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2.h:399
real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2.h:505
Treal template_blas_fabs(Treal x)
ergo_real real
Definition: test.cc:46
void truncate_matrix(real &threshold, int it)
Definition: purification_sp2.h:232
#define LOG_CAT_INFO
Definition: output.h:49
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
void return_constant_C(const int it, real &Cval)
Definition: purification_sp2.h:223
Purification_sp2acc is a class which provides an interface for SP2 recursive expansion.
Definition: purification_sp2.h:51
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2.h:60
void set_poly(const int it)
Definition: purification_sp2.h:94
int method
Definition: puri_info.h:171
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
generalVector VectorType
Definition: purification_sp2.h:62
Recursive density matrix expansion (or density matrix purification).
Definition: puri_info.h:52
void set_init_params()
Definition: purification_sp2.h:66
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2.h:56
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:246
Definition: matInclude.h:139
void estimate_number_of_iterations(int &numit)
Definition: purification_sp2.h:278
void apply_poly_to_eigs(const int it, real &homo, real &lumo)
Definition: purification_sp2.h:479
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
real gap
Definition: puri_info.h:79
void purify_bounds(const int it)
Definition: purification_sp2.h:169
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2.h:55
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:139