39#ifndef HEADER_PURIFICATION_SP2ACC
40#define HEADER_PURIFICATION_SP2ACC
50template<
typename MatrixType>
69 this->
info.method = 2;
76 virtual void get_poly(
const int it,
int& poly,
real& alpha);
99template<
typename MatrixType>
105template<
typename MatrixType>
108 assert((
int)this->
VecPoly.size() > it);
113 real Xtrace = this->
X.trace();
114 real Xsqtrace = this->
Xsq.trace();
165template<
typename MatrixType>
168 assert((
int)this->
VecPoly.size() > it);
169 assert(this->
VecPoly[it] != -1);
179template<
typename MatrixType>
215 this->
X *= ((
real)2.0 * (1 - alpha_tmp) * alpha_tmp);
216 this->
X.add_identity((
real)(1 - alpha_tmp) * (1 - alpha_tmp));
217 this->
Xsq *= ((
real)alpha_tmp * alpha_tmp);
218 this->
Xsq += this->
X;
231 this->
X *= ((
real) - 2.0 * alpha_tmp);
232 this->
Xsq *= ((
real) - alpha_tmp * alpha_tmp);
233 this->
Xsq -= this->
X;
238 this->
X *= (
real)2.0;
239 this->
Xsq += this->
X;
245 this->
Xsq.transfer(this->
X);
249template<
typename MatrixType>
256 real homo_low, homo_upp, lumo_upp, lumo_low;
267 lumo_low = (1 - alpha_tmp + alpha_tmp * this->
lumo_bounds.low());
268 lumo_low *= lumo_low;
269 lumo_upp = (1 - alpha_tmp + alpha_tmp * this->
lumo_bounds.upp());
270 lumo_upp *= lumo_upp;
280 homo_low = (1 - alpha_tmp + alpha_tmp * this->
homo_bounds.low());
281 homo_low *= homo_low;
282 homo_upp = (1 - alpha_tmp + alpha_tmp * this->
homo_bounds.upp());
283 homo_upp *= homo_upp;
312template<
typename MatrixType>
331 if (((alpha1 == 1) && (alpha2 == 1)))
348 real homo_low = this->
info.Iterations[it - 2].homo_bound_low;
349 real homo_upp = this->
info.Iterations[it - 2].homo_bound_upp;
350 real lumo_low = this->
info.Iterations[it - 2].lumo_bound_low;
351 real lumo_upp = this->
info.Iterations[it - 2].lumo_bound_upp;
355 if ((homo_upp > 0.5) || (lumo_upp > 0.5))
362 a = std::max(lumo_low, homo_low);
370 " and alpha2 = %g", a, alpha1, alpha2);
376 C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
391template<
typename MatrixType>
396 int maxit_tmp = maxit_total;
397 real x, y, x_out, y_out;
401 int max_size = maxit_total + 1 + 2;
404 this->
VecPoly.resize(max_size, -1);
407 this->
VecGap.resize(max_size, -1);
422 estim_num_iter = this->
maxit;
438 this->
VecGap[0] = 1 - x - y;
442 while (it <= maxit_tmp)
445 if ((x > y) || (it % 2 && (x == y)))
447 alpha_tmp = 2 / (2 - x_out);
449 x = (1 - alpha_tmp + alpha_tmp * x);
451 y = 2 * alpha_tmp * y - alpha_tmp * alpha_tmp * y * y;
453 x_out = (1 - alpha_tmp + alpha_tmp * x_out);
455 y_out = 2 * alpha_tmp * y_out - alpha_tmp * alpha_tmp * y_out * y_out;
461 alpha_tmp = 2 / (2 - y_out);
463 x = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
464 y = (1 - alpha_tmp + alpha_tmp * y);
467 x_out = 2 * alpha_tmp * x_out - alpha_tmp * alpha_tmp * x_out * x_out;
468 y_out = (1 - alpha_tmp + alpha_tmp * y_out);
475 this->
VecGap[it] = 1 - x - y;
497 (x - x * x < epsilon) && (y - y * y < epsilon) &&
501 maxit_tmp = it + this->additional_iterations;
517 if ( ((estim_num_iter == -1) && (it == maxit_tmp + 1) )
518 || (estim_num_iter > maxit_total))
521 estim_num_iter = this->
maxit;
522 maxit_tmp = maxit_total;
525 this->
VecPoly.resize(maxit_tmp + 1);
526 this->
VecGap.resize(maxit_tmp + 1);
527 this->
VecAlpha.resize(maxit_tmp + 1);
533template<
typename MatrixType>
536 assert((
int)this->
VecPoly.size() > it);
537 assert((
int)this->
VecGap.size() > it);
538 assert((
int)this->
VecAlpha.size() > it);
548template<
typename MatrixType>
555 for (
int i = it; i >= 1; i--)
564 bounds_from_it[0] = (bounds_from_it[0] - 1 + alpha_tmp) / alpha_tmp - tau;
566 bounds_from_it[1] = (bounds_from_it[1] - 1 + alpha_tmp) / alpha_tmp + tau;
568 bounds_from_it[2] = bounds_from_it[2] / (1 +
template_blas_sqrt(1 - bounds_from_it[2]));
569 bounds_from_it[2] = bounds_from_it[2] / alpha_tmp - tau;
570 bounds_from_it[3] = bounds_from_it[3] / (1 +
template_blas_sqrt(1 - bounds_from_it[3]));
571 bounds_from_it[3] = bounds_from_it[3] / alpha_tmp + tau;
575 bounds_from_it[0] = bounds_from_it[0] / (1 +
template_blas_sqrt(1 - bounds_from_it[0]));
576 bounds_from_it[0] = bounds_from_it[0] / alpha_tmp - tau;
577 bounds_from_it[1] = bounds_from_it[1] / (1 +
template_blas_sqrt(1 - bounds_from_it[1]));
578 bounds_from_it[1] = bounds_from_it[1] / alpha_tmp + tau;
581 bounds_from_it[2] = (bounds_from_it[2] - 1 + alpha_tmp) / alpha_tmp - tau;
583 bounds_from_it[3] = (bounds_from_it[3] - 1 + alpha_tmp) / alpha_tmp + tau;
589template<
typename MatrixType>
606 fx = (1 - alpha_tmp + alpha_tmp * x) * (1 - alpha_tmp + alpha_tmp * x);
610 fx = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
618template<
typename MatrixType>
633 for (
int i = 1; i <= it; i++)
641 a = ((1 - alpha) + alpha * temp) * ((1 - alpha) + alpha * temp);
642 b = 2 * alpha * ((1 - alpha) + alpha * temp);
646 a = 2 * alpha * temp - alpha * alpha * temp * temp;
647 b = 2 * alpha - 2 * alpha * alpha * temp;
Definition puri_info.h:52
real gap
Definition puri_info.h:82
int poly
Definition puri_info.h:81
real alpha
Definition puri_info.h:94
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition purification_general.h:482
VectorTypeReal VecGap
Gap computed using inner homo and lumo bounds on each iteration.
Definition purification_general.h:513
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
int maxit
Maximum number of iterations.
Definition purification_general.h:481
NormType normPuriStopCrit
Norm used in the stopping criterion Can be mat::frobNorm, mat::mixedNorm, or mat::euclNorm.
Definition purification_general.h:492
IntervalType homo_bounds
(1-homo) bounds for Xi in iteration i
Definition purification_general.h:531
static real get_epsilon()
Get machine epsilon.
Definition purification_general.h:196
MatrixType Xsq
Matrix X^2.
Definition purification_general.h:132
ergo_real real
Definition purification_general.h:119
PurificationGeneral()
Definition purification_general.h:135
PuriInfo info
Fill in during purification with useful information.
Definition purification_general.h:128
VectorTypeInt VecPoly
Polynomials computed in the function estimated_number_of_iterations() VecPoly[i] = 1 if we use X=X^2 ...
Definition purification_general.h:508
IntervalType lumo_bounds
Lumo bounds for Xi in iteration i.
Definition purification_general.h:532
intervalType IntervalType
Definition purification_general.h:120
int additional_iterations
Do a few more iterations after convergence.
Definition purification_general.h:479
MatrixType X
Matrix X.
Definition purification_general.h:130
int nocc
Number of occupied orbitals.
Definition purification_general.h:484
mat::normType NormType
Definition purification_general.h:121
VectorTypeReal VecAlpha
Definition purification_sp2acc.h:93
PurificationGeneral< MatrixType >::NormType NormType
Definition purification_sp2acc.h:57
virtual real compute_derivative(const int it, real x, real &DDf)
Definition purification_sp2acc.h:620
virtual void purify_X(const int it)
Definition purification_sp2acc.h:180
virtual void get_poly(const int it, int &poly, real &alpha)
Definition purification_sp2acc.h:166
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition purification_sp2acc.h:59
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition purification_sp2acc.h:56
static const real deltaTurnOffAcc
Definition purification_sp2acc.h:96
virtual void set_init_params()
Definition purification_sp2acc.h:66
virtual real apply_poly(const int it, real x)
Definition purification_sp2acc.h:591
virtual void estimate_number_of_iterations(int &numit)
Definition purification_sp2acc.h:392
virtual void save_other_iter_info(IterationInfo &iter_info, int it)
Definition purification_sp2acc.h:534
PurificationGeneral< MatrixType >::real real
Definition purification_sp2acc.h:55
Purification_sp2acc()
Definition purification_sp2acc.h:64
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition purification_sp2acc.h:60
virtual void purify_bounds(const int it)
Definition purification_sp2acc.h:250
generalVector VectorType
Definition purification_sp2acc.h:62
virtual void set_poly(const int it)
Definition purification_sp2acc.h:106
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition purification_sp2acc.h:549
virtual void return_constant_C(const int it, real &Cval)
Definition purification_sp2acc.h:313
#define C_SP2
Constant used on the stopping criterion for the SP2 method.
Definition constants.h:42
ergo_real real
Definition test.cc:46
mat::VectorGeneral< ergo_real, Vectorrr > generalVector
Definition matrix_typedefs.h:76
@ 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)