ergo
purification_sp2acc.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
39#ifndef HEADER_PURIFICATION_SP2ACC
40#define HEADER_PURIFICATION_SP2ACC
41
43
44//#define DEBUG_OUTPUT
45
50template<typename MatrixType>
51class Purification_sp2acc : public PurificationGeneral<MatrixType>
52{
53public:
54
58
61
63
65
66 virtual void set_init_params()
67 {
68 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen SP2 ACCELERATED purification method");
69 this->info.method = 2;
70
71 this->gammaStopEstim = 6 - 4 * template_blas_sqrt((real)2);
72
73 this->check_stopping_criterion_iter = -1; // will be changed during purification
74 }
75
76 virtual void get_poly(const int it, int& poly, real& alpha);
77 virtual void set_poly(const int it);
78
79 virtual void estimate_number_of_iterations(int& numit);
80 virtual void purify_X(const int it);
81 virtual void purify_bounds(const int it);
82 virtual void save_other_iter_info(IterationInfo& iter_info, int it);
83 virtual void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it);
84
85 virtual void return_constant_C(const int it, real& Cval);
86
87 virtual real apply_poly(const int it, real x);
88 virtual real compute_derivative(const int it, real x, real& DDf);
89
90
91 /* PARAMETERS */
92
94
95 // defined the iteration when we turn off acceleration
96 static const real deltaTurnOffAcc;
97};
98
99template<typename MatrixType>
102
103
104
105template<typename MatrixType>
107{
108 assert((int)this->VecPoly.size() > it);
109
110 // if cannot compute polynomial using homo and lumo eigevalues, compute using trace
111 if (this->VecPoly[it] == -1)
112 {
113 real Xtrace = this->X.trace();
114 real Xsqtrace = this->Xsq.trace();
115
116 real delta = deltaTurnOffAcc;
117
118 // Should we turn off acceleration or not
119 if ((this->check_stopping_criterion_iter == -1) && (this->lumo_bounds.low() < delta) && (this->homo_bounds.low() < delta))
120 {
121#ifdef DEBUG_OUTPUT
122 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Outer bounds of homo and lumo are less then %e: ", delta);
123 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo_out = %e, homo_out = %e ", this->lumo_bounds.low(), this->homo_bounds.low());
124#endif
125
126 this->lumo_bounds = IntervalType(0, this->lumo_bounds.upp());
127 this->homo_bounds = IntervalType(0, this->homo_bounds.upp());
128
129 // start to check stopping criterion
130 if (it == 1)
131 {
132 this->check_stopping_criterion_iter = it + 1; // in the it=0 we had the same eigenvalue bounds
133 }
134 else
135 {
136 this->check_stopping_criterion_iter = it + 2;
137 }
138 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
139 }
140
141 if ((template_blas_fabs(Xsqtrace - this->nocc) <
142 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
143 ||
144 (it % 2
145 &&
146 (template_blas_fabs(Xsqtrace - this->nocc) ==
147 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
148 ))
149 {
150 this->VecPoly[it] = 1;
151 VecAlpha[it] = 2 / (2 - this->lumo_bounds.low());
152 }
153 else
154 {
155 this->VecPoly[it] = 0;
156 VecAlpha[it] = 2 / (2 - this->homo_bounds.low());
157 }
158#ifdef DEBUG_OUTPUT
159 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Acceleration parameter: alpha = %lf", VecAlpha[it]);
160#endif
161 }
162}
163
164
165template<typename MatrixType>
166void Purification_sp2acc<MatrixType>::get_poly(const int it, int& poly, real& alpha)
167{
168 assert((int)this->VecPoly.size() > it);
169 assert(this->VecPoly[it] != -1);
170
171 //check also if alpha is computed
172 assert(this->VecAlpha[it] != -1);
173
174 poly = this->VecPoly[it];
175 alpha = VecAlpha[it];
176}
177
178
179template<typename MatrixType>
181{
182#ifdef DEBUG_OUTPUT
184#endif
185 real alpha_tmp;
186 int poly;
187
188 set_poly(it);
189
190 get_poly(it, poly, alpha_tmp);
191
192 /* It may happen that X2 has many more nonzeros than X, for
193 * example 5 times as many. Therefore it makes sense to try
194 * having only one "big" matrix in memory at a time. However,
195 * file operations have proved to be quite expensive and should
196 * be avoided if possible. Hence we want to achieve having only
197 * one big matrix in memory without unnecessary file
198 * operations. We are currently hoping that it will be ok to add
199 * a "small" matrix to a "big" one, that the memory usage after
200 * that operation will be like the memory usage for one big
201 * matrix + one small matrix. Therefore we are adding X to X2 (X
202 * is truncated, a "small" matrix) instead of the opposite.
203 */
204
205 if (poly == 1)
206 {
207 if (alpha_tmp != 1)
208 {
209 // (1-a+a*x)^2 = (1-a)^2 + 2*(1-a)*a*x + a^2*x^2
210 // this->X.mult_scalar((real)2.0 * (1 - alpha_tmp) * alpha_tmp);
211 // this->X.add_identity((real)(1 - alpha_tmp) * (1 - alpha_tmp));
212 // this->Xsq.mult_scalar((real)alpha_tmp * alpha_tmp);
213 // this->Xsq.add(this->X); // Xsq = (1-a+a*X)^2
214
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; // Xsq = (1-a+a*X)^2
219
220
221 }
222 else
223 {
224 // DO NOTHING
225 }
226 }
227 else
228 {
229 if (alpha_tmp != 1)
230 {
231 this->X *= ((real) - 2.0 * alpha_tmp);
232 this->Xsq *= ((real) - alpha_tmp * alpha_tmp);
233 this->Xsq -= this->X; // Xsq = 2*a*X - (a*X)^2
234 }
235 else
236 {
237 this->Xsq *= ((real) - 1.0);
238 this->X *= (real)2.0;
239 this->Xsq += this->X; // Xsq = -Xsq + 2X
240
241
242 }
243 } // if poly == 1
244
245 this->Xsq.transfer(this->X); // clear Xsq and old X
246}
247
248
249template<typename MatrixType>
251{
252#ifdef DEBUG_OUTPUT
253 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change homo and lumo bounds according to the chosen polynomial VecPoly = %d", this->VecPoly[it]);
254#endif
255
256 real homo_low, homo_upp, lumo_upp, lumo_low;
257 real alpha_tmp;
258 int poly;
259
260 get_poly(it, poly, alpha_tmp);
261
262 if (poly == 1)
263 {
264 // update bounds
265 homo_low = 2 * alpha_tmp * this->homo_bounds.low() - alpha_tmp * alpha_tmp * this->homo_bounds.low() * this->homo_bounds.low(); // 2*a*x - (a*x)^2
266 homo_upp = 2 * alpha_tmp * this->homo_bounds.upp() - alpha_tmp * alpha_tmp * this->homo_bounds.upp() * this->homo_bounds.upp(); // 2*a*x - (a*x)^2
267 lumo_low = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.low()); // (1-a+a*x)^2
268 lumo_low *= lumo_low;
269 lumo_upp = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.upp()); // (1-a+a*x)^2
270 lumo_upp *= lumo_upp;
271
272 this->homo_bounds = IntervalType(homo_low, homo_upp);
273 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
274 }
275 else
276 {
277 // update bounds
278 lumo_low = 2 * alpha_tmp * this->lumo_bounds.low() - alpha_tmp * alpha_tmp * this->lumo_bounds.low() * this->lumo_bounds.low(); // 2*a*x - (a*x)^2
279 lumo_upp = 2 * alpha_tmp * this->lumo_bounds.upp() - alpha_tmp * alpha_tmp * this->lumo_bounds.upp() * this->lumo_bounds.upp(); // 2*a*x - (a*x)^2
280 homo_low = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.low()); // (1-a+a*x)^2
281 homo_low *= homo_low;
282 homo_upp = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.upp()); // (1-a+a*x)^2
283 homo_upp *= homo_upp;
284
285 this->homo_bounds = IntervalType(homo_low, homo_upp);
286 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
287 }
288
289 IntervalType zero_one(0, 1);
290 this->homo_bounds.intersect(zero_one);
291 this->lumo_bounds.intersect(zero_one);
292
293#ifdef DEBUG_OUTPUT
294 if (this->homo_bounds.empty())
295 {
296 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
297 }
298 if (this->lumo_bounds.empty())
299 {
300 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
301 }
302
303
304 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "1-homo: [ %g , %g ],", this->homo_bounds.low(), this->homo_bounds.upp());
305 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo: [ %g , %g ].", this->lumo_bounds.low(), this->lumo_bounds.upp());
306#endif
307}
308
309
310/*****************************************************/
311
312template<typename MatrixType>
314{
315 assert(it >= 1);
316
317 real alpha1 = VecAlpha[it - 1];
318 real alpha2 = VecAlpha[it];
319
320 Cval = -1;
321
322#ifdef DEBUG_OUTPUT
323 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "alpha1 = %.4e , alpha2 = %.4e", alpha1, alpha2);
324#endif
325
326 if (it < 2)
327 {
328 return; // -1
329 }
330 // no acceleration
331 if (((alpha1 == 1) && (alpha2 == 1)))
332 {
333#ifdef DEBUG_OUTPUT
334 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Check SP2 stopping criterion.");
335#endif
336 Cval = C_SP2;
337 return;
338 }
339#ifdef DEBUG_OUTPUT
340 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Do not check stopping criterion.");
341#endif
342 // exit and return C = -1
343
344
345 // If we want to compute the constant C on every iterations, we can use the following:
346#if 0
347 // get bounds from the iteration it-2
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;
352 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo [%.16e, %e], homo [%.16e, %e]", lumo_low, lumo_upp, homo_low, homo_upp);
353
354
355 if ((homo_upp > 0.5) || (lumo_upp > 0.5))
356 {
357 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Inner bounds interval do not contain 0.5. Skip iteration. ");
358 Cval = -1;
359 return;
360 }
361
362 a = std::max(lumo_low, homo_low);
363
364 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen a = %g", a);
365
366 if (a <= 0)
367 {
368 // just in case
369 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot compute constant C since a = %g when alpha1 = %g"
370 " and alpha2 = %g", a, alpha1, alpha2);
371 Cval = -1;
372 return;
373 }
374
375 real C1;
376 C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
377 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Local maximum of g1: %g", C1);
378
379 Cval = C1;
380
381 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "*********** C = %g ************", Cval);
382#endif
383}
384
385
386
387/****************************************************************************************/
388
389
390
391template<typename MatrixType>
393{
394 int it = 1;
395 int maxit_total = this->maxit + this->additional_iterations;
396 int maxit_tmp = maxit_total;
397 real x, y, x_out, y_out;
398 real alpha_tmp;
399 real epsilon = this->get_epsilon();
400
401 int max_size = maxit_total + 1 + 2; // largest possible vector size, +1 is because we save initial iteration 0, +2 is because we might use Frobenius norm and then we add two extra iterations
402
403 this->VecPoly.clear();
404 this->VecPoly.resize(max_size, -1);
405
406 this->VecGap.clear();
407 this->VecGap.resize(max_size, -1);
408
409 VecAlpha.clear();
410 VecAlpha.resize(max_size, -1);
411
412 // we are interested in the inner bounds of gap
413 x = this->lumo_bounds.upp(); // = lumo
414 y = this->homo_bounds.upp(); // = 1 - homo
415
416// if VecGap is zero
417 if (1 - x - y <= 0)
418 {
419#ifdef DEBUG_OUTPUT
420 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap cannot be computed. Set estimated number of iteration to the maxit.");
421#endif
422 estim_num_iter = this->maxit;
423 return;
424 }
425
426 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT LUMO: [ %.12lf , %.12lf ]", (double)this->lumo_bounds.low(), (double)this->lumo_bounds.upp());
427 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT HOMO: [ %.12lf , %.12lf ]", (double)this->homo_bounds.low(), (double)this->homo_bounds.upp());
428
429
430 x_out = this->lumo_bounds.low();
431 y_out = this->homo_bounds.low();
432
433
434
435 real delta = deltaTurnOffAcc;
436
437 this->VecPoly[0] = -1;
438 this->VecGap[0] = 1 - x - y;
439
440 estim_num_iter = -1;
441
442 while (it <= maxit_tmp)
443 {
444 // note: avoid not-stopping in case of idempotent matrix
445 if ((x > y) || (it % 2 && (x == y))) // lumo > 1-homo
446 {
447 alpha_tmp = 2 / (2 - x_out);
448
449 x = (1 - alpha_tmp + alpha_tmp * x);
450 x *= x;
451 y = 2 * alpha_tmp * y - alpha_tmp * alpha_tmp * y * y;
452
453 x_out = (1 - alpha_tmp + alpha_tmp * x_out);
454 x_out *= x_out;
455 y_out = 2 * alpha_tmp * y_out - alpha_tmp * alpha_tmp * y_out * y_out;
456
457 this->VecPoly[it] = 1;
458 }
459 else
460 {
461 alpha_tmp = 2 / (2 - y_out);
462
463 x = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
464 y = (1 - alpha_tmp + alpha_tmp * y);
465 y *= y;
466
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);
469 y_out *= y_out;
470
471 this->VecPoly[it] = 0;
472 }
473
474 VecAlpha[it] = alpha_tmp;
475 this->VecGap[it] = 1 - x - y;
476
477 // find iteration where x_out < delta && y_out < delta
478 if ((x_out < delta) && (y_out < delta) && (this->check_stopping_criterion_iter == -1))
479 {
480 // start to check stopping criterion
481 if (it == 1)
482 {
483 this->check_stopping_criterion_iter = it + 1; // in the it=0 we had the same eigenvalue bounds
484 }
485 else
486 {
487 this->check_stopping_criterion_iter = it + 2;
488 }
489 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
490 x_out = 0;
491 y_out = 0;
492 }
493
494
495 // maybe we wish to perform some more iterations, then stopping criterion suggest
496 if ((estim_num_iter == -1) && (it >= this->check_stopping_criterion_iter) &&
497 (x - x * x < epsilon) && (y - y * y < epsilon) && // the eucledian norm is less then epsilon
498 (this->VecPoly[it] != this->VecPoly[it - 1])) // to apply stopping criterion, polynomials must alternate
499 {
500 estim_num_iter = it;
501 maxit_tmp = it + this->additional_iterations;
502
503 // if we use Frobenius norm, it seems that sometimes we need one or two iterations more than what is suggested by the stopping criterion (which assumes spectral norm)
504 if (this->normPuriStopCrit == mat::frobNorm)
505 {
506 estim_num_iter += 2;
507 maxit_tmp += 2;
508 }
509 }
510
511 ++it;
512 } //while
513
514 /*
515 Either we reached maxit number of iterations or due to the additional 2 iterations (because of the Frobenius norm) we overestimated number of iterations
516 */
517 if ( ((estim_num_iter == -1) && (it == maxit_tmp + 1) )
518 || (estim_num_iter > maxit_total))
519 {
520 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxit = %d number of iterations is reached in estimate_number_of_iteration()", this->maxit);
521 estim_num_iter = this->maxit;
522 maxit_tmp = maxit_total;
523 }
524
525 this->VecPoly.resize(maxit_tmp + 1); // make it less if needed
526 this->VecGap.resize(maxit_tmp + 1);
527 this->VecAlpha.resize(maxit_tmp + 1);
528}
529
530
531/************ SAVE INFORMATION ABOUT ITERATIONS SPECIFIC FOR ACC SP2 PURIFICATION **********x***/
532
533template<typename MatrixType>
535{
536 assert((int)this->VecPoly.size() > it);
537 assert((int)this->VecGap.size() > it);
538 assert((int)this->VecAlpha.size() > it);
539
540 iter_info.poly = this->VecPoly[it];
541 iter_info.gap = this->VecGap[it];
542 iter_info.alpha = VecAlpha[it];
543}
544
545
546/************ APPLY INVERSE POLYNOMIAL (FOR ESTIMATION OF EIGENVALUES) ************/
547
548template<typename MatrixType>
550{
551 real tau;
552 real alpha_tmp;
553 int poly;
554
555 for (int i = it; i >= 1; i--)
556 {
557 tau = 0;//this->info.Iterations[i].threshold_X;
558
559 get_poly(i, poly, alpha_tmp);
560
561 if (poly == 1)
562 {
563 bounds_from_it[0] = template_blas_sqrt(bounds_from_it[0]);
564 bounds_from_it[0] = (bounds_from_it[0] - 1 + alpha_tmp) / alpha_tmp - tau;
565 bounds_from_it[1] = template_blas_sqrt(bounds_from_it[1]);
566 bounds_from_it[1] = (bounds_from_it[1] - 1 + alpha_tmp) / alpha_tmp + tau;
567
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;
572 }
573 else
574 {
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;
579
580 bounds_from_it[2] = template_blas_sqrt(bounds_from_it[2]);
581 bounds_from_it[2] = (bounds_from_it[2] - 1 + alpha_tmp) / alpha_tmp - tau;
582 bounds_from_it[3] = template_blas_sqrt(bounds_from_it[3]);
583 bounds_from_it[3] = (bounds_from_it[3] - 1 + alpha_tmp) / alpha_tmp + tau;
584 }
585 }
586}
587
588
589template<typename MatrixType>
592{
593 assert(it >= 0);
594 if (it == 0)
595 {
596 return x;
597 }
598
599 real fx;
600 int poly;
601 real alpha_tmp;
602 get_poly(it, poly, alpha_tmp);
603
604 if (poly == 1)
605 {
606 fx = (1 - alpha_tmp + alpha_tmp * x) * (1 - alpha_tmp + alpha_tmp * x);
607 }
608 else
609 {
610 fx = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
611 }
612
613 return fx;
614}
615
616
617
618template<typename MatrixType>
621{
622 assert(it > 0);
623
624 real Df;
625 real temp, a, b;
626 int poly;
627 real alpha;
628
629 a = x;
630 Df = 1;
631 DDf = -1; // TODO
632
633 for (int i = 1; i <= it; i++)
634 {
635 temp = a;
636
637 get_poly(i, poly, alpha);
638
639 if (poly == 1)
640 {
641 a = ((1 - alpha) + alpha * temp) * ((1 - alpha) + alpha * temp);
642 b = 2 * alpha * ((1 - alpha) + alpha * temp);
643 }
644 else
645 {
646 a = 2 * alpha * temp - alpha * alpha * temp * temp;
647 b = 2 * alpha - 2 * alpha * alpha * temp;
648 }
649 Df *= b;
650 }
651
652 return Df;
653}
654
655
656#endif //HEADER_PURIFICATION_SP2ACC
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 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
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition purification_general.h:482
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 SP2ACC recursive expansion.
Definition purification_sp2acc.h:52
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
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
static const real deltaTurnOffAcc
Definition purification_sp2acc.h:96
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
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)