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);
101 real Xtrace = this->
X.trace();
102 real Xsqtrace = this->
Xsq.trace();
123 template<
typename MatrixType>
126 assert((
int)this->
VecPoly.size() > it);
127 assert(this->
VecPoly[it] != -1);
132 template<
typename MatrixType>
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;
222 template<
typename MatrixType>
231 template<
typename MatrixType>
237 if (allowed_error == 0)
242 assert((
int)this->
VecGap.size() > it);
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);
277 template<
typename MatrixType>
288 int max_size = maxit_tmp + 1 + this->additional_iterations + 2;
291 this->
VecPoly.resize(max_size, -1);
294 this->
VecGap.resize(max_size, -1);
307 estim_num_iter = this->
maxit;
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]))
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);
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 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 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
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition: purification_general.h:436
Treal low() const
Definition: Interval.h:144
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 apply_poly(const int it, real x)
Definition: purification_sp2.h:453
static real get_epsilon()
Get machine epsilon.
Definition: purification_general.h:198
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
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 check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:405
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)
int additional_iterations
Do a few more iterations after convergence.
Definition: purification_general.h:402
int nocc
Number of occupied orbitals.
Definition: purification_general.h:407
void truncate_matrix(real &threshold, int it)
Definition: purification_sp2.h:232
#define LOG_CAT_INFO
Definition: output.h:49
NormType normPuriStopCrit
Norm used in the stopping criterion Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm.
Definition: purification_general.h:415
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
MatrixType Xsq
Matrix X^2.
Definition: purification_general.h:250
int method
Definition: puri_info.h:171
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
generalVector VectorType
Definition: purification_sp2.h:62
MatrixType X
Matrix X.
Definition: purification_general.h:248
int maxit
Maximum number of iterations.
Definition: purification_general.h:404
Recursive density matrix expansion (or density matrix purification).
NormType normPuriTrunc
Norm used for the truncation of matrices.
Definition: purification_general.h:411
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
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
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
bool empty() const
Definition: Interval.h:51
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