ergo
purification_sp2acc.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
42 #include "purification_general.h"
43 
44 //#define DEBUG_OUTPUT
45 
50 template<typename MatrixType>
51 class Purification_sp2acc : public PurificationGeneral<MatrixType>
52 {
53 public:
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 
99 template<typename MatrixType>
102 
103 
104 
105 template<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 
165 template<typename MatrixType>
166 void 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 
179 template<typename MatrixType>
181 {
182 #ifdef DEBUG_OUTPUT
183  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Purify X...");
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 
249 template<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 
312 template<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 
391 template<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 
533 template<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 
548 template<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 
589 template<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 
618 template<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
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
C_SP2
#define C_SP2
Constant used on the stopping criterion for the SP2 method.
Definition: constants.h:42
PurificationGeneral
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:117
PurificationGeneral::info
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:128
Purification_sp2acc
Purification_sp2acc is a class which provides an interface for SP2ACC recursive expansion.
Definition: purification_sp2acc.h:52
PurificationGeneral::gammaStopEstim
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:503
purification_general.h
Recursive density matrix expansion (or density matrix purification).
PuriInfo::method
int method
Definition: puri_info.h:189
PurificationGeneral::real
ergo_real real
Definition: purification_general.h:119
IntervalType
intervalType IntervalType
Definition: random_matrices.h:71
Purification_sp2acc::apply_poly
virtual real apply_poly(const int it, real x)
Definition: purification_sp2acc.h:591
mat::frobNorm
@ frobNorm
Definition: matInclude.h:139
real
ergo_real real
Definition: test.cc:46
IterationInfo::alpha
real alpha
Definition: puri_info.h:94
template_blas_fabs
Treal template_blas_fabs(Treal x)
Purification_sp2acc::VecAlpha
VectorTypeReal VecAlpha
Definition: purification_sp2acc.h:93
Purification_sp2acc::purify_X
virtual void purify_X(const int it)
Definition: purification_sp2acc.h:180
IterationInfo::poly
int poly
Definition: puri_info.h:81
LOG_AREA_DENSFROMF
#define LOG_AREA_DENSFROMF
Definition: output.h:61
Purification_sp2acc::set_poly
virtual void set_poly(const int it)
Definition: purification_sp2acc.h:106
Purification_sp2acc::Purification_sp2acc
Purification_sp2acc()
Definition: purification_sp2acc.h:64
IterationInfo
Definition: puri_info.h:52
MatrixType
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
Purification_sp2acc::estimate_number_of_iterations
virtual void estimate_number_of_iterations(int &numit)
Definition: purification_sp2acc.h:392
mat::VectorGeneral
Definition: MatrixBase.h:55
Purification_sp2acc::VectorType
generalVector VectorType
Definition: purification_sp2acc.h:62
Purification_sp2acc::apply_inverse_poly_vector
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2acc.h:549
PurificationGeneral::VectorTypeInt
std::vector< int > VectorTypeInt
Definition: purification_general.h:122
LOG_CAT_INFO
#define LOG_CAT_INFO
Definition: output.h:49
Purification_sp2acc::save_other_iter_info
virtual void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2acc.h:534
mat::normType
normType
Definition: matInclude.h:139
mat::Interval< ergo_real >
max
#define max(a, b)
Definition: integrator.cc:87
Purification_sp2acc::purify_bounds
virtual void purify_bounds(const int it)
Definition: purification_sp2acc.h:250
Purification_sp2acc::deltaTurnOffAcc
static const real deltaTurnOffAcc
Definition: purification_sp2acc.h:96
PurificationGeneral::VectorTypeReal
std::vector< real > VectorTypeReal
Definition: purification_general.h:123
Purification_sp2acc::return_constant_C
virtual void return_constant_C(const int it, real &Cval)
Definition: purification_sp2acc.h:313
Purification_sp2acc::get_poly
virtual void get_poly(const int it, int &poly, real &alpha)
Definition: purification_sp2acc.h:166
IterationInfo::gap
real gap
Definition: puri_info.h:82
Purification_sp2acc::VectorTypeReal
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2acc.h:60
Purification_sp2acc::real
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2acc.h:55
do_output
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
PurificationGeneral::check_stopping_criterion_iter
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:482
Purification_sp2acc::IntervalType
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2acc.h:56
Purification_sp2acc::NormType
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2acc.h:57
Purification_sp2acc::compute_derivative
virtual real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2acc.h:620
Purification_sp2acc::set_init_params
virtual void set_init_params()
Definition: purification_sp2acc.h:66
Purification_sp2acc::VectorTypeInt
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2acc.h:59