39 #ifndef HEADER_PURIFICATION_SP2ACC 40 #define HEADER_PURIFICATION_SP2ACC 50 template<
typename MatrixType>
76 virtual void get_poly(
const int it,
int& poly,
real& alpha);
102 template<
typename MatrixType>
108 template<
typename MatrixType>
111 assert((
int)this->VecPoly.size() > it);
114 if (this->VecPoly[it] == -1)
116 real Xtrace = this->X.trace();
117 real Xsqtrace = this->Xsq.trace();
119 real delta = deltaTurnOffAcc;
122 if ((this->check_stopping_criterion_iter == -1) && (this->lumo_bounds.low() < delta) && (this->homo_bounds.low() < delta))
129 this->lumo_bounds =
IntervalType(0, this->lumo_bounds.upp());
135 this->check_stopping_criterion_iter = it + 1;
139 this->check_stopping_criterion_iter = it + 2;
153 this->VecPoly[it] = 1;
154 VecAlpha[it] = 2 / (2 - this->lumo_bounds.low());
158 this->VecPoly[it] = 0;
159 VecAlpha[it] = 2 / (2 - this->homo_bounds.low());
168 template<
typename MatrixType>
171 assert((
int)this->VecPoly.size() > it);
172 assert(this->VecPoly[it] != -1);
175 assert(this->VecAlpha[it] != -1);
177 poly = this->VecPoly[it];
178 alpha = VecAlpha[it];
182 template<
typename MatrixType>
193 get_poly(it, poly, alpha_tmp);
218 this->X *= ((
real)2.0 * (1 - alpha_tmp) * alpha_tmp);
219 this->X.add_identity((
real)(1 - alpha_tmp) * (1 - alpha_tmp));
220 this->Xsq *= ((
real)alpha_tmp * alpha_tmp);
221 this->Xsq += this->X;
234 this->X *= ((
real) - 2.0 * alpha_tmp);
235 this->Xsq *= ((
real) - alpha_tmp * alpha_tmp);
236 this->Xsq -= this->X;
240 this->Xsq *= ((
real) - 1.0);
241 this->X *= (
real)2.0;
242 this->Xsq += this->X;
248 this->Xsq.transfer(this->X);
252 template<
typename MatrixType>
259 real homo_low, homo_upp, lumo_upp, lumo_low;
263 get_poly(it, poly, alpha_tmp);
268 homo_low = 2 * alpha_tmp * this->homo_bounds.low() - alpha_tmp * alpha_tmp * this->homo_bounds.low() * this->homo_bounds.low();
269 homo_upp = 2 * alpha_tmp * this->homo_bounds.upp() - alpha_tmp * alpha_tmp * this->homo_bounds.upp() * this->homo_bounds.upp();
270 lumo_low = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.low());
271 lumo_low *= lumo_low;
272 lumo_upp = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.upp());
273 lumo_upp *= lumo_upp;
281 lumo_low = 2 * alpha_tmp * this->lumo_bounds.low() - alpha_tmp * alpha_tmp * this->lumo_bounds.low() * this->lumo_bounds.low();
282 lumo_upp = 2 * alpha_tmp * this->lumo_bounds.upp() - alpha_tmp * alpha_tmp * this->lumo_bounds.upp() * this->lumo_bounds.upp();
283 homo_low = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.low());
284 homo_low *= homo_low;
285 homo_upp = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.upp());
286 homo_upp *= homo_upp;
293 this->homo_bounds.intersect(zero_one);
294 this->lumo_bounds.intersect(zero_one);
297 if (this->homo_bounds.empty())
301 if (this->lumo_bounds.empty())
315 template<
typename MatrixType>
320 real alpha1 = VecAlpha[it - 1];
321 real alpha2 = VecAlpha[it];
334 if (((alpha1 == 1) && (alpha2 == 1)))
351 real homo_low = this->info.Iterations[it - 2].homo_bound_low;
352 real homo_upp = this->info.Iterations[it - 2].homo_bound_upp;
353 real lumo_low = this->info.Iterations[it - 2].lumo_bound_low;
354 real lumo_upp = this->info.Iterations[it - 2].lumo_bound_upp;
358 if ((homo_upp > 0.5) || (lumo_upp > 0.5))
373 " and alpha2 = %g", a, alpha1, alpha2);
379 C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
391 template<
typename MatrixType>
394 real allowed_error = this->error_per_it;
397 if (allowed_error == 0)
402 assert((
int)this->VecGap.size() > it);
408 if (this->VecGap[it] > 0)
413 tau = (allowed_error * this->VecGap[it]) / (1 + allowed_error);
417 tau = allowed_error * 0.01;
422 #ifdef USE_CHUNKS_AND_TASKS 423 threshold = (this->X).thresh_frob(tau);
425 threshold = (this->X).thresh(tau, this->normPuriTrunc);
438 template<
typename MatrixType>
442 int maxit_tmp = this->maxit + this->additional_iterations;
443 real x, y, x_out, y_out;
445 real epsilon = this->get_epsilon();
447 int max_size = maxit_tmp + 1 + this->additional_iterations + 2;
449 this->VecPoly.clear();
450 this->VecPoly.resize(max_size, -1);
452 this->VecGap.clear();
453 this->VecGap.resize(max_size, -1);
456 VecAlpha.resize(max_size, -1);
459 x = this->lumo_bounds.upp();
460 y = this->homo_bounds.upp();
468 estim_num_iter = this->maxit;
476 x_out = this->lumo_bounds.low();
477 y_out = this->homo_bounds.low();
481 real delta = deltaTurnOffAcc;
483 this->VecPoly[0] = -1;
484 this->VecGap[0] = 1 - x - y;
488 while (it <= maxit_tmp)
491 if ((x > y) || (it % 2 && (x == y)))
493 alpha_tmp = 2 / (2 - x_out);
495 x = (1 - alpha_tmp + alpha_tmp * x);
497 y = 2 * alpha_tmp * y - alpha_tmp * alpha_tmp * y * y;
499 x_out = (1 - alpha_tmp + alpha_tmp * x_out);
501 y_out = 2 * alpha_tmp * y_out - alpha_tmp * alpha_tmp * y_out * y_out;
503 this->VecPoly[it] = 1;
507 alpha_tmp = 2 / (2 - y_out);
509 x = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
510 y = (1 - alpha_tmp + alpha_tmp * y);
513 x_out = 2 * alpha_tmp * x_out - alpha_tmp * alpha_tmp * x_out * x_out;
514 y_out = (1 - alpha_tmp + alpha_tmp * y_out);
517 this->VecPoly[it] = 0;
520 VecAlpha[it] = alpha_tmp;
521 this->VecGap[it] = 1 - x - y;
524 if ((x_out < delta) && (y_out < delta) && (this->check_stopping_criterion_iter == -1))
529 this->check_stopping_criterion_iter = it + 1;
533 this->check_stopping_criterion_iter = it + 2;
542 if ((estim_num_iter == -1) && (it >= this->check_stopping_criterion_iter) &&
543 (x - x * x < epsilon) && (y - y * y < epsilon) &&
544 (this->VecPoly[it] != this->VecPoly[it - 1]))
547 maxit_tmp = it + this->additional_iterations;
559 if ((estim_num_iter == -1) && (it == maxit_tmp + 1))
564 estim_num_iter = this->maxit;
567 this->VecPoly.resize(maxit_tmp + 1);
568 this->VecGap.resize(maxit_tmp + 1);
569 this->VecAlpha.resize(maxit_tmp + 1);
575 template<
typename MatrixType>
578 assert((
int)this->VecPoly.size() > it);
579 assert((
int)this->VecGap.size() > it);
580 assert((
int)this->VecAlpha.size() > it);
582 iter_info.
poly = this->VecPoly[it];
583 iter_info.
gap = this->VecGap[it];
584 iter_info.
alpha = VecAlpha[it];
590 template<
typename MatrixType>
597 for (
int i = it; i >= 1; i--)
601 get_poly(i, poly, alpha_tmp);
606 bounds_from_it[0] = (bounds_from_it[0] - 1 + alpha_tmp) / alpha_tmp - tau;
608 bounds_from_it[1] = (bounds_from_it[1] - 1 + alpha_tmp) / alpha_tmp + tau;
610 bounds_from_it[2] = bounds_from_it[2] / (1 +
template_blas_sqrt(1 - bounds_from_it[2]));
611 bounds_from_it[2] = bounds_from_it[2] / alpha_tmp - tau;
612 bounds_from_it[3] = bounds_from_it[3] / (1 +
template_blas_sqrt(1 - bounds_from_it[3]));
613 bounds_from_it[3] = bounds_from_it[3] / alpha_tmp + tau;
617 bounds_from_it[0] = bounds_from_it[0] / (1 +
template_blas_sqrt(1 - bounds_from_it[0]));
618 bounds_from_it[0] = bounds_from_it[0] / alpha_tmp - tau;
619 bounds_from_it[1] = bounds_from_it[1] / (1 +
template_blas_sqrt(1 - bounds_from_it[1]));
620 bounds_from_it[1] = bounds_from_it[1] / alpha_tmp + tau;
623 bounds_from_it[2] = (bounds_from_it[2] - 1 + alpha_tmp) / alpha_tmp - tau;
625 bounds_from_it[3] = (bounds_from_it[3] - 1 + alpha_tmp) / alpha_tmp + tau;
631 template<
typename MatrixType>
644 get_poly(it, poly, alpha_tmp);
648 fx = (1 - alpha_tmp + alpha_tmp * x) * (1 - alpha_tmp + alpha_tmp * x);
652 fx = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
659 template<
typename MatrixType>
670 get_poly(it, poly, alpha_tmp);
674 homo = 2 * alpha_tmp * homo - alpha_tmp * alpha_tmp * homo * homo;
675 lumo = (1 - alpha_tmp + alpha_tmp * lumo) * (1 - alpha_tmp + alpha_tmp * lumo);
679 homo = (1 - alpha_tmp + alpha_tmp * homo) * (1 - alpha_tmp + alpha_tmp * homo);
680 lumo = 2 * alpha_tmp * lumo - alpha_tmp * alpha_tmp * lumo * lumo;
685 template<
typename MatrixType>
700 for (
int i = 1; i <= it; i++)
704 get_poly(i, poly, alpha);
708 a = ((1 - alpha) + alpha * temp) * ((1 - alpha) + alpha * temp);
709 b = 2 * alpha * ((1 - alpha) + alpha * temp);
713 a = 2 * alpha * temp - alpha * alpha * temp * temp;
714 b = 2 * alpha - 2 * alpha * alpha * temp;
723 #endif //HEADER_PURIFICATION_SP2ACC Treal upp() const
Definition: Interval.h:145
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
VectorTypeReal VecAlpha
Definition: purification_sp2acc.h:96
virtual void purify_X(const int it)
Definition: purification_sp2acc.h:183
virtual void apply_poly_to_eigs(const int it, real &homo, real &lumo)
Definition: purification_sp2acc.h:660
intervalType IntervalType
Definition: random_matrices.h:67
virtual void set_poly(const int it)
Definition: purification_sp2acc.h:109
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:111
generalVector VectorType
Definition: purification_sp2acc.h:62
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
real alpha
Definition: puri_info.h:91
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2acc.h:591
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:405
int poly
Definition: puri_info.h:78
#define LOG_AREA_DENSFROMF
Definition: output.h:61
virtual void set_init_params()
Definition: purification_sp2acc.h:66
std::vector< int > VectorTypeInt
Definition: purification_general.h:117
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2acc.h:56
virtual void return_constant_C(const int it, real &Cval)
Definition: purification_sp2acc.h:316
Treal template_blas_fabs(Treal x)
ergo_real real
Definition: test.cc:46
virtual void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2acc.h:576
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2acc.h:59
#define max(a, b)
Definition: integrator.cc:87
Purification_sp2acc is a class which provides an interface for SP2ACC recursive expansion.
Definition: purification_sp2acc.h:51
#define LOG_CAT_INFO
Definition: output.h:49
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
static const real deltaTurnOffAcc
Definition: purification_sp2acc.h:99
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2acc.h:55
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2acc.h:57
int method
Definition: puri_info.h:171
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
Recursive density matrix expansion (or density matrix purification).
Definition: puri_info.h:52
virtual void truncate_matrix(real &threshold, int it)
Definition: purification_sp2acc.h:392
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:246
Definition: matInclude.h:139
virtual real apply_poly(const int it, real x)
Definition: purification_sp2acc.h:633
virtual real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2acc.h:687
Purification_sp2acc()
Definition: purification_sp2acc.h:64
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
real gap
Definition: puri_info.h:79
virtual void estimate_number_of_iterations(int &numit)
Definition: purification_sp2acc.h:439
virtual void get_poly(const int it, int &poly, real &alpha)
Definition: purification_sp2acc.h:169
virtual void purify_bounds(const int it)
Definition: purification_sp2acc.h:253
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2acc.h:60
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:139