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);
99 template<
typename MatrixType>
105 template<
typename MatrixType>
108 assert((
int)this->VecPoly.size() > it);
111 if (this->VecPoly[it] == -1)
113 real Xtrace = this->X.trace();
114 real Xsqtrace = this->Xsq.trace();
116 real delta = deltaTurnOffAcc;
119 if ((this->check_stopping_criterion_iter == -1) && (this->lumo_bounds.low() < delta) && (this->homo_bounds.low() < delta))
126 this->lumo_bounds =
IntervalType(0, this->lumo_bounds.upp());
127 this->homo_bounds =
IntervalType(0, this->homo_bounds.upp());
132 this->check_stopping_criterion_iter = it + 1;
136 this->check_stopping_criterion_iter = it + 2;
150 this->VecPoly[it] = 1;
151 VecAlpha[it] = 2 / (2 - this->lumo_bounds.low());
155 this->VecPoly[it] = 0;
156 VecAlpha[it] = 2 / (2 - this->homo_bounds.low());
165 template<
typename MatrixType>
168 assert((
int)this->VecPoly.size() > it);
169 assert(this->VecPoly[it] != -1);
172 assert(this->VecAlpha[it] != -1);
174 poly = this->VecPoly[it];
175 alpha = VecAlpha[it];
179 template<
typename MatrixType>
190 get_poly(it, poly, alpha_tmp);
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;
237 this->Xsq *= ((
real) - 1.0);
238 this->X *= (
real)2.0;
239 this->Xsq += this->X;
245 this->Xsq.transfer(this->X);
249 template<
typename MatrixType>
256 real homo_low, homo_upp, lumo_upp, lumo_low;
260 get_poly(it, poly, alpha_tmp);
265 homo_low = 2 * alpha_tmp * this->homo_bounds.low() - alpha_tmp * alpha_tmp * this->homo_bounds.low() * this->homo_bounds.low();
266 homo_upp = 2 * alpha_tmp * this->homo_bounds.upp() - alpha_tmp * alpha_tmp * this->homo_bounds.upp() * this->homo_bounds.upp();
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;
278 lumo_low = 2 * alpha_tmp * this->lumo_bounds.low() - alpha_tmp * alpha_tmp * this->lumo_bounds.low() * this->lumo_bounds.low();
279 lumo_upp = 2 * alpha_tmp * this->lumo_bounds.upp() - alpha_tmp * alpha_tmp * this->lumo_bounds.upp() * this->lumo_bounds.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;
290 this->homo_bounds.intersect(zero_one);
291 this->lumo_bounds.intersect(zero_one);
294 if (this->homo_bounds.empty())
298 if (this->lumo_bounds.empty())
312 template<
typename MatrixType>
317 real alpha1 = VecAlpha[it - 1];
318 real alpha2 = VecAlpha[it];
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))
370 " and alpha2 = %g", a, alpha1, alpha2);
376 C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
391 template<
typename MatrixType>
395 int maxit_total = this->maxit + this->additional_iterations;
396 int maxit_tmp = maxit_total;
397 real x, y, x_out, y_out;
399 real epsilon = this->get_epsilon();
401 int max_size = maxit_total + 1 + 2;
403 this->VecPoly.clear();
404 this->VecPoly.resize(max_size, -1);
406 this->VecGap.clear();
407 this->VecGap.resize(max_size, -1);
410 VecAlpha.resize(max_size, -1);
413 x = this->lumo_bounds.upp();
414 y = this->homo_bounds.upp();
422 estim_num_iter = this->maxit;
430 x_out = this->lumo_bounds.low();
431 y_out = this->homo_bounds.low();
435 real delta = deltaTurnOffAcc;
437 this->VecPoly[0] = -1;
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;
457 this->VecPoly[it] = 1;
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);
471 this->VecPoly[it] = 0;
474 VecAlpha[it] = alpha_tmp;
475 this->VecGap[it] = 1 - x - y;
478 if ((x_out < delta) && (y_out < delta) && (this->check_stopping_criterion_iter == -1))
483 this->check_stopping_criterion_iter = it + 1;
487 this->check_stopping_criterion_iter = it + 2;
496 if ((estim_num_iter == -1) && (it >= this->check_stopping_criterion_iter) &&
497 (x - x * x < epsilon) && (y - y * y < epsilon) &&
498 (this->VecPoly[it] != this->VecPoly[it - 1]))
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);
533 template<
typename MatrixType>
536 assert((
int)this->VecPoly.size() > it);
537 assert((
int)this->VecGap.size() > it);
538 assert((
int)this->VecAlpha.size() > it);
540 iter_info.
poly = this->VecPoly[it];
541 iter_info.
gap = this->VecGap[it];
542 iter_info.
alpha = VecAlpha[it];
548 template<
typename MatrixType>
555 for (
int i = it; i >= 1; i--)
559 get_poly(i, poly, alpha_tmp);
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;
589 template<
typename MatrixType>
602 get_poly(it, poly, alpha_tmp);
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;
618 template<
typename MatrixType>
633 for (
int i = 1; i <= it; i++)
637 get_poly(i, poly, alpha);
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;
656 #endif //HEADER_PURIFICATION_SP2ACC