39#ifndef HEADER_PURIFICATION_SP2
40#define HEADER_PURIFICATION_SP2
50template<
typename MatrixType>
74 void get_poly(
const int it,
int& poly);
90template<
typename MatrixType>
93 assert((
int)this->VecPoly.size() > it);
96 if (this->VecPoly[it] == -1)
98 real Xtrace = this->X.trace();
99 real Xsqtrace = this->Xsq.trace();
110 this->VecPoly[it] = 1;
114 this->VecPoly[it] = 0;
120template<
typename MatrixType>
123 assert((
int)this->VecPoly.size() > it);
124 assert(this->VecPoly[it] != -1);
125 poly = this->VecPoly[it];
129template<
typename MatrixType>
156 this->Xsq *= ((
real) - 1.0);
157 this->X *= (
real)2.0;
158 this->Xsq += this->X;
161 this->Xsq.transfer(this->X);
165template<
typename MatrixType>
171 real homo_low, homo_upp, lumo_upp, lumo_low;
179 homo_low = 2 * this->homo_bounds.low() - this->homo_bounds.low() * this->homo_bounds.low();
180 homo_upp = 2 * this->homo_bounds.upp() - this->homo_bounds.upp() * this->homo_bounds.upp();
181 lumo_low = this->lumo_bounds.low() * this->lumo_bounds.low();
182 lumo_upp = this->lumo_bounds.upp() * this->lumo_bounds.upp();
189 lumo_low = 2 * this->lumo_bounds.low() - this->lumo_bounds.low() * this->lumo_bounds.low();
190 lumo_upp = 2 * this->lumo_bounds.upp() - this->lumo_bounds.upp() * this->lumo_bounds.upp();
191 homo_low = this->homo_bounds.low() * this->homo_bounds.low();
192 homo_upp = this->homo_bounds.upp() * this->homo_bounds.upp();
198 this->homo_bounds.intersect(zero_one);
199 this->lumo_bounds.intersect(zero_one);
202 if (this->homo_bounds.empty())
206 if (this->lumo_bounds.empty())
219template<
typename MatrixType>
229template<
typename MatrixType>
233 int maxit_total = this->maxit + this->additional_iterations;
234 int maxit_tmp = maxit_total;
236 real epsilon = this->get_epsilon();
239 this->check_stopping_criterion_iter = 2;
241 int max_size = maxit_total + 1 + 2;
242 this->VecPoly.clear();
243 this->VecPoly.resize(max_size, -1);
245 this->VecGap.clear();
246 this->VecGap.resize(max_size, -1);
249 x = this->lumo_bounds.upp();
250 y = this->homo_bounds.upp();
259 estim_num_iter = this->maxit;
268 this->VecPoly[0] = -1;
269 this->VecGap[0] = 1 - x - y;
273 while (it <= maxit_tmp)
276 if ((x > y) || (it % 2 && (x == y)))
280 this->VecPoly[it] = 1;
286 this->VecPoly[it] = 0;
289 this->VecGap[it] = 1 - x - y;
292 if ((estim_num_iter == -1) &&
293 (x - x * x < epsilon) && (y - y * y < epsilon)
296 (this->VecPoly[it] != this->VecPoly[it - 1]))
299 maxit_tmp = it + this->additional_iterations;
316 if ( ((estim_num_iter == -1) && (it == maxit_tmp + 1) )
317 || (estim_num_iter > maxit_total))
320 estim_num_iter = this->maxit;
321 maxit_tmp = maxit_total;
325 this->VecPoly.resize(maxit_tmp + 1);
326 this->VecGap.resize(maxit_tmp + 1);
330 for (
int i = 0; i < (int)this->VecPoly.size(); ++i)
341template<
typename MatrixType>
344 assert((
int)this->VecPoly.size() > it);
345 assert((
int)this->VecGap.size() > it);
347 iter_info.
poly = this->VecPoly[it];
348 iter_info.
gap = this->VecGap[it];
355template<
typename MatrixType>
359 for (
int i = it; i >= 1; i--)
368 bounds_from_it[2] = bounds_from_it[2] / (1 +
template_blas_sqrt(1 - bounds_from_it[2]));
369 bounds_from_it[3] = bounds_from_it[3] / (1 +
template_blas_sqrt(1 - bounds_from_it[3]));
373 bounds_from_it[0] = bounds_from_it[0] / (1 +
template_blas_sqrt(1 - bounds_from_it[0]));
374 bounds_from_it[1] = bounds_from_it[1] / (1 +
template_blas_sqrt(1 - bounds_from_it[1]));
408template<
typename MatrixType>
437template<
typename MatrixType>
451 for (
int i = 1; i <= it; i++)
457 DDf = 2 * Df * Df + (2 * a) * DDf;
463 DDf = -2 * Df * Df + (2 - 2 * a) * DDf;
Definition: puri_info.h:52
real gap
Definition: puri_info.h:82
int poly
Definition: puri_info.h:81
int method
Definition: puri_info.h:189
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:117
std::vector< real > VectorTypeReal
Definition: purification_general.h:123
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:503
std::vector< int > VectorTypeInt
Definition: purification_general.h:122
ergo_real real
Definition: purification_general.h:119
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:128
Purification_sp2acc is a class which provides an interface for SP2 recursive expansion.
Definition: purification_sp2.h:52
void set_init_params()
Definition: purification_sp2.h:66
void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2.h:356
real apply_poly(const int it, real x)
Definition: purification_sp2.h:410
real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2.h:439
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2.h:57
void return_constant_C(const int it, real &Cval)
Definition: purification_sp2.h:220
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2.h:55
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2.h:56
void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2.h:342
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2.h:60
generalVector VectorType
Definition: purification_sp2.h:62
void get_poly(const int it, int &poly)
Definition: purification_sp2.h:121
void set_poly(const int it)
Definition: purification_sp2.h:91
Purification_sp2()
Definition: purification_sp2.h:64
void estimate_number_of_iterations(int &numit)
Definition: purification_sp2.h:230
void purify_X(const int it)
Definition: purification_sp2.h:130
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2.h:59
void purify_bounds(const int it)
Definition: purification_sp2.h:166
Definition: VectorGeneral.h:48
#define C_SP2
Constant used on the stopping criterion for the SP2 method.
Definition: constants.h:42
ergo_real real
Definition: test.cc:46
normType
Definition: matInclude.h:139
@ frobNorm
Definition: matInclude.h:139
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
#define LOG_AREA_DENSFROMF
Definition: output.h:61
#define LOG_CAT_INFO
Definition: output.h:49
Recursive density matrix expansion (or density matrix purification).
intervalType IntervalType
Definition: random_matrices.h:71
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)