ergo
purification_sp2acc.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 truncate_matrix(real& threshold, int it);
80 
81  virtual void estimate_number_of_iterations(int& numit);
82  virtual void purify_X(const int it);
83  virtual void purify_bounds(const int it);
84  virtual void save_other_iter_info(IterationInfo& iter_info, int it);
85  virtual void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it);
86 
87  virtual void return_constant_C(const int it, real& Cval);
88 
89  virtual real apply_poly(const int it, real x);
90  virtual void apply_poly_to_eigs(const int it, real& homo, real& lumo);
91  virtual real compute_derivative(const int it, real x, real& DDf);
92 
93 
94  /* PARAMETERS */
95 
97 
98  // defined the iteration when we turn off acceleration
99  static const real deltaTurnOffAcc;
100 };
101 
102 template<typename MatrixType>
105 
106 
107 
108 template<typename MatrixType>
110 {
111  assert((int)this->VecPoly.size() > it);
112 
113  // if cannot compute polynomial using homo and lumo eigevalues, compute using trace
114  if (this->VecPoly[it] == -1)
115  {
116  real Xtrace = this->X.trace();
117  real Xsqtrace = this->Xsq.trace();
118 
119  real delta = deltaTurnOffAcc;
120 
121  // Should we turn off acceleration or not
122  if ((this->check_stopping_criterion_iter == -1) && (this->lumo_bounds.low() < delta) && (this->homo_bounds.low() < delta))
123  {
124 #ifdef DEBUG_OUTPUT
125  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Outer bounds of homo and lumo are less then %e: ", delta);
126  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo_out = %e, homo_out = %e ", this->lumo_bounds.low(), this->homo_bounds.low());
127 #endif
128 
129  this->lumo_bounds = IntervalType(0, this->lumo_bounds.upp());
130  this->homo_bounds = IntervalType(0, this->homo_bounds.upp());
131 
132  // start to check stopping criterion
133  if (it == 1)
134  {
135  this->check_stopping_criterion_iter = it + 1; // in the it=0 we had the same eigenvalue bounds
136  }
137  else
138  {
139  this->check_stopping_criterion_iter = it + 2;
140  }
141  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
142  }
143 
144  if ((template_blas_fabs(Xsqtrace - this->nocc) <
145  template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
146  ||
147  (it % 2
148  &&
149  (template_blas_fabs(Xsqtrace - this->nocc) ==
150  template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
151  ))
152  {
153  this->VecPoly[it] = 1;
154  VecAlpha[it] = 2 / (2 - this->lumo_bounds.low());
155  }
156  else
157  {
158  this->VecPoly[it] = 0;
159  VecAlpha[it] = 2 / (2 - this->homo_bounds.low());
160  }
161 #ifdef DEBUG_OUTPUT
162  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Acceleration parameter: alpha = %lf", VecAlpha[it]);
163 #endif
164  }
165 }
166 
167 
168 template<typename MatrixType>
169 void Purification_sp2acc<MatrixType>::get_poly(const int it, int& poly, real& alpha)
170 {
171  assert((int)this->VecPoly.size() > it);
172  assert(this->VecPoly[it] != -1);
173 
174  //check also if alpha is computed
175  assert(this->VecAlpha[it] != -1);
176 
177  poly = this->VecPoly[it];
178  alpha = VecAlpha[it];
179 }
180 
181 
182 template<typename MatrixType>
184 {
185 #ifdef DEBUG_OUTPUT
186  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Purify X...");
187 #endif
188  real alpha_tmp;
189  int poly;
190 
191  set_poly(it);
192 
193  get_poly(it, poly, alpha_tmp);
194 
195  /* It may happen that X2 has many more nonzeros than X, for
196  * example 5 times as many. Therefore it makes sense to try
197  * having only one "big" matrix in memory at a time. However,
198  * file operations have proved to be quite expensive and should
199  * be avoided if possible. Hence we want to achieve having only
200  * one big matrix in memory without unnecessary file
201  * operations. We are currently hoping that it will be ok to add
202  * a "small" matrix to a "big" one, that the memory usage after
203  * that operation will be like the memory usage for one big
204  * matrix + one small matrix. Therefore we are adding X to X2 (X
205  * is truncated, a "small" matrix) instead of the opposite.
206  */
207 
208  if (poly == 1)
209  {
210  if (alpha_tmp != 1)
211  {
212  // (1-a+a*x)^2 = (1-a)^2 + 2*(1-a)*a*x + a^2*x^2
213  // this->X.mult_scalar((real)2.0 * (1 - alpha_tmp) * alpha_tmp);
214  // this->X.add_identity((real)(1 - alpha_tmp) * (1 - alpha_tmp));
215  // this->Xsq.mult_scalar((real)alpha_tmp * alpha_tmp);
216  // this->Xsq.add(this->X); // Xsq = (1-a+a*X)^2
217 
218  this->X *= ((real)2.0 * (1 - alpha_tmp) * alpha_tmp);
219  this->X.add_identity((real)(1 - alpha_tmp) * (1 - alpha_tmp));
220  this->Xsq *= ((real)alpha_tmp * alpha_tmp);
221  this->Xsq += this->X; // Xsq = (1-a+a*X)^2
222 
223 
224  }
225  else
226  {
227  // DO NOTHING
228  }
229  }
230  else
231  {
232  if (alpha_tmp != 1)
233  {
234  this->X *= ((real) - 2.0 * alpha_tmp);
235  this->Xsq *= ((real) - alpha_tmp * alpha_tmp);
236  this->Xsq -= this->X; // Xsq = 2*a*X - (a*X)^2
237  }
238  else
239  {
240  this->Xsq *= ((real) - 1.0);
241  this->X *= (real)2.0;
242  this->Xsq += this->X; // Xsq = -Xsq + 2X
243 
244 
245  }
246  } // if poly == 1
247 
248  this->Xsq.transfer(this->X); // clear Xsq and old X
249 }
250 
251 
252 template<typename MatrixType>
254 {
255 #ifdef DEBUG_OUTPUT
256  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change homo and lumo bounds according to the chosen polynomial VecPoly = %d", this->VecPoly[it]);
257 #endif
258 
259  real homo_low, homo_upp, lumo_upp, lumo_low;
260  real alpha_tmp;
261  int poly;
262 
263  get_poly(it, poly, alpha_tmp);
264 
265  if (poly == 1)
266  {
267  // update bounds
268  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
269  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
270  lumo_low = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.low()); // (1-a+a*x)^2
271  lumo_low *= lumo_low;
272  lumo_upp = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.upp()); // (1-a+a*x)^2
273  lumo_upp *= lumo_upp;
274 
275  this->homo_bounds = IntervalType(homo_low, homo_upp);
276  this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
277  }
278  else
279  {
280  // update bounds
281  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
282  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
283  homo_low = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.low()); // (1-a+a*x)^2
284  homo_low *= homo_low;
285  homo_upp = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.upp()); // (1-a+a*x)^2
286  homo_upp *= homo_upp;
287 
288  this->homo_bounds = IntervalType(homo_low, homo_upp);
289  this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
290  }
291 
292  IntervalType zero_one(0, 1);
293  this->homo_bounds.intersect(zero_one);
294  this->lumo_bounds.intersect(zero_one);
295 
296 #ifdef DEBUG_OUTPUT
297  if (this->homo_bounds.empty())
298  {
299  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
300  }
301  if (this->lumo_bounds.empty())
302  {
303  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
304  }
305 
306 
307  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "1-homo: [ %g , %g ],", this->homo_bounds.low(), this->homo_bounds.upp());
308  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo: [ %g , %g ].", this->lumo_bounds.low(), this->lumo_bounds.upp());
309 #endif
310 }
311 
312 
313 /*****************************************************/
314 
315 template<typename MatrixType>
317 {
318  assert(it >= 1);
319 
320  real alpha1 = VecAlpha[it - 1];
321  real alpha2 = VecAlpha[it];
322 
323  Cval = -1;
324 
325 #ifdef DEBUG_OUTPUT
326  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "alpha1 = %.4e , alpha2 = %.4e", alpha1, alpha2);
327 #endif
328 
329  if (it < 2)
330  {
331  return; // -1
332  }
333  // no acceleration
334  if (((alpha1 == 1) && (alpha2 == 1)))
335  {
336 #ifdef DEBUG_OUTPUT
337  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Check SP2 stopping criterion.");
338 #endif
339  Cval = C_SP2;
340  return;
341  }
342 #ifdef DEBUG_OUTPUT
343  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Do not check stopping criterion.");
344 #endif
345  // exit and return C = -1
346 
347 
348  // If we want to compute the constant C on every iterations, we can use the following:
349 #if 0
350  // get bounds from the iteration it-2
351  real homo_low = this->info.Iterations[it - 2].homo_bound_low;
352  real homo_upp = this->info.Iterations[it - 2].homo_bound_upp;
353  real lumo_low = this->info.Iterations[it - 2].lumo_bound_low;
354  real lumo_upp = this->info.Iterations[it - 2].lumo_bound_upp;
355  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo [%.16e, %e], homo [%.16e, %e]", lumo_low, lumo_upp, homo_low, homo_upp);
356 
357 
358  if ((homo_upp > 0.5) || (lumo_upp > 0.5))
359  {
360  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Inner bounds interval do not contain 0.5. Skip iteration. ");
361  Cval = -1;
362  return;
363  }
364 
365  a = std::max(lumo_low, homo_low);
366 
367  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen a = %g", a);
368 
369  if (a <= 0)
370  {
371  // just in case
372  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot compute constant C since a = %g when alpha1 = %g"
373  " and alpha2 = %g", a, alpha1, alpha2);
374  Cval = -1;
375  return;
376  }
377 
378  real C1;
379  C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
380  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Local maximum of g1: %g", C1);
381 
382  Cval = C1;
383 
384  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "*********** C = %g ************", Cval);
385 #endif
386 }
387 
388 
389 /**************************************************************************************/
390 
391 template<typename MatrixType>
393 {
394  real allowed_error = this->error_per_it;
395 
396  threshold = 0;
397  if (allowed_error == 0)
398  {
399  return;
400  }
401 
402  assert((int)this->VecGap.size() > it);
403 #ifdef DEBUG_OUTPUT
404  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Truncate matrix X: ");
405 #endif
406  real tau;
407 
408  if (this->VecGap[it] > 0) // TODO: zero or some small value?
409  {
410 #ifdef DEBUG_OUTPUT
411  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap[ %d ] = %e , ", it, this->VecGap[it]);
412 #endif
413  tau = (allowed_error * this->VecGap[it]) / (1 + allowed_error); // control error in the occupied subspace
414  }
415  else // if gap is not known
416  {
417  tau = allowed_error * 0.01;
418  }
419 
420 
421 // TODO
422  #ifdef USE_CHUNKS_AND_TASKS
423  threshold = (this->X).thresh_frob(tau);
424  #else
425  threshold = (this->X).thresh(tau, this->normPuriTrunc);
426  #endif
427 
428 #ifdef DEBUG_OUTPUT
429  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "tau = %e", tau);
430 #endif
431 }
432 
433 
434 /****************************************************************************************/
435 
436 
437 
438 template<typename MatrixType>
440 {
441  int it = 1;
442  int maxit_tmp = this->maxit + this->additional_iterations;
443  real x, y, x_out, y_out;
444  real alpha_tmp;
445  real epsilon = this->get_epsilon();
446 
447  int max_size = maxit_tmp + 1 + this->additional_iterations + 2; // largest possible size
448 
449  this->VecPoly.clear();
450  this->VecPoly.resize(max_size, -1);
451 
452  this->VecGap.clear();
453  this->VecGap.resize(max_size, -1);
454 
455  VecAlpha.clear();
456  VecAlpha.resize(max_size, -1);
457 
458  // we are interested in the inner bounds of gap
459  x = this->lumo_bounds.upp(); // = lumo
460  y = this->homo_bounds.upp(); // = 1 - homo
461 
462 // if VecGap is zero
463  if (1 - x - y <= 0)
464  {
465 #ifdef DEBUG_OUTPUT
466  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap cannot be computed. Set estimated number of iteration to the maxit.");
467 #endif
468  estim_num_iter = this->maxit;
469  return;
470  }
471 
472  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT LUMO: [ %.12lf , %.12lf ]", (double)this->lumo_bounds.low(), (double)this->lumo_bounds.upp());
473  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT HOMO: [ %.12lf , %.12lf ]", (double)this->homo_bounds.low(), (double)this->homo_bounds.upp());
474 
475 
476  x_out = this->lumo_bounds.low();
477  y_out = this->homo_bounds.low();
478 
479 
480 
481  real delta = deltaTurnOffAcc;
482 
483  this->VecPoly[0] = -1;
484  this->VecGap[0] = 1 - x - y;
485 
486  estim_num_iter = -1;
487 
488  while (it <= maxit_tmp)
489  {
490  // note: avoid not-stopping in case of idempotent matrix
491  if ((x > y) || (it % 2 && (x == y))) // lumo > 1-homo
492  {
493  alpha_tmp = 2 / (2 - x_out);
494 
495  x = (1 - alpha_tmp + alpha_tmp * x);
496  x *= x;
497  y = 2 * alpha_tmp * y - alpha_tmp * alpha_tmp * y * y;
498 
499  x_out = (1 - alpha_tmp + alpha_tmp * x_out);
500  x_out *= x_out;
501  y_out = 2 * alpha_tmp * y_out - alpha_tmp * alpha_tmp * y_out * y_out;
502 
503  this->VecPoly[it] = 1;
504  }
505  else
506  {
507  alpha_tmp = 2 / (2 - y_out);
508 
509  x = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
510  y = (1 - alpha_tmp + alpha_tmp * y);
511  y *= y;
512 
513  x_out = 2 * alpha_tmp * x_out - alpha_tmp * alpha_tmp * x_out * x_out;
514  y_out = (1 - alpha_tmp + alpha_tmp * y_out);
515  y_out *= y_out;
516 
517  this->VecPoly[it] = 0;
518  }
519 
520  VecAlpha[it] = alpha_tmp;
521  this->VecGap[it] = 1 - x - y;
522 
523  // find iteration where x_out < delta && y_out < delta
524  if ((x_out < delta) && (y_out < delta) && (this->check_stopping_criterion_iter == -1))
525  {
526  // start to check stopping criterion
527  if (it == 1)
528  {
529  this->check_stopping_criterion_iter = it + 1; // in the it=0 we had the same eigenvalue bounds
530  }
531  else
532  {
533  this->check_stopping_criterion_iter = it + 2;
534  }
535  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
536  x_out = 0;
537  y_out = 0;
538  }
539 
540 
541  // maybe we wish to perform some more iterations, then stopping criterion suggest
542  if ((estim_num_iter == -1) && (it >= this->check_stopping_criterion_iter) &&
543  (x - x * x < epsilon) && (y - y * y < epsilon) && // the eucledian norm is less then epsilon
544  (this->VecPoly[it] != this->VecPoly[it - 1])) // to apply stopping criterion, polynomials must alternate
545  {
546  estim_num_iter = it;
547  maxit_tmp = it + this->additional_iterations;
548 
549  if (this->normPuriStopCrit == mat::frobNorm)
550  {
551  estim_num_iter += 2;
552  maxit_tmp += 2;
553  }
554  }
555 
556  ++it;
557  } //while
558 
559  if ((estim_num_iter == -1) && (it == maxit_tmp + 1))
560  {
561 #ifdef DEBUG_OUTPUT
562  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxit = %d number of iterations is reached in estimate_number_of_iteration()", this->maxit);
563 #endif
564  estim_num_iter = this->maxit;
565  }
566 
567  this->VecPoly.resize(maxit_tmp + 1); // make it less if needed
568  this->VecGap.resize(maxit_tmp + 1);
569  this->VecAlpha.resize(maxit_tmp + 1);
570 }
571 
572 
573 /************ SAVE INFORMATION ABOUT ITERATIONS SPECIFIC FOR ACC SP2 PURIFICATION **********x***/
574 
575 template<typename MatrixType>
577 {
578  assert((int)this->VecPoly.size() > it);
579  assert((int)this->VecGap.size() > it);
580  assert((int)this->VecAlpha.size() > it);
581 
582  iter_info.poly = this->VecPoly[it];
583  iter_info.gap = this->VecGap[it];
584  iter_info.alpha = VecAlpha[it];
585 }
586 
587 
588 /************ APPLY INVERSE POLYNOMIAL (FOR ESTIMATION OF EIGENVALUES) ************/
589 
590 template<typename MatrixType>
592 {
593  real tau;
594  real alpha_tmp;
595  int poly;
596 
597  for (int i = it; i >= 1; i--)
598  {
599  tau = 0;//this->info.Iterations[i].threshold_X;
600 
601  get_poly(i, poly, alpha_tmp);
602 
603  if (poly == 1)
604  {
605  bounds_from_it[0] = template_blas_sqrt(bounds_from_it[0]);
606  bounds_from_it[0] = (bounds_from_it[0] - 1 + alpha_tmp) / alpha_tmp - tau;
607  bounds_from_it[1] = template_blas_sqrt(bounds_from_it[1]);
608  bounds_from_it[1] = (bounds_from_it[1] - 1 + alpha_tmp) / alpha_tmp + tau;
609 
610  bounds_from_it[2] = bounds_from_it[2] / (1 + template_blas_sqrt(1 - bounds_from_it[2]));
611  bounds_from_it[2] = bounds_from_it[2] / alpha_tmp - tau;
612  bounds_from_it[3] = bounds_from_it[3] / (1 + template_blas_sqrt(1 - bounds_from_it[3]));
613  bounds_from_it[3] = bounds_from_it[3] / alpha_tmp + tau;
614  }
615  else
616  {
617  bounds_from_it[0] = bounds_from_it[0] / (1 + template_blas_sqrt(1 - bounds_from_it[0]));
618  bounds_from_it[0] = bounds_from_it[0] / alpha_tmp - tau;
619  bounds_from_it[1] = bounds_from_it[1] / (1 + template_blas_sqrt(1 - bounds_from_it[1]));
620  bounds_from_it[1] = bounds_from_it[1] / alpha_tmp + tau;
621 
622  bounds_from_it[2] = template_blas_sqrt(bounds_from_it[2]);
623  bounds_from_it[2] = (bounds_from_it[2] - 1 + alpha_tmp) / alpha_tmp - tau;
624  bounds_from_it[3] = template_blas_sqrt(bounds_from_it[3]);
625  bounds_from_it[3] = (bounds_from_it[3] - 1 + alpha_tmp) / alpha_tmp + tau;
626  }
627  }
628 }
629 
630 
631 template<typename MatrixType>
634 {
635  assert(it >= 0);
636  if (it == 0)
637  {
638  return x;
639  }
640 
641  real fx;
642  int poly;
643  real alpha_tmp;
644  get_poly(it, poly, alpha_tmp);
645 
646  if (poly == 1)
647  {
648  fx = (1 - alpha_tmp + alpha_tmp * x) * (1 - alpha_tmp + alpha_tmp * x);
649  }
650  else
651  {
652  fx = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
653  }
654 
655  return fx;
656 }
657 
658 
659 template<typename MatrixType>
661 {
662  assert(it >= 0);
663  if (it == 0)
664  {
665  return;
666  }
667 
668  int poly;
669  real alpha_tmp;
670  get_poly(it, poly, alpha_tmp);
671 
672  if (poly == 1)
673  {
674  homo = 2 * alpha_tmp * homo - alpha_tmp * alpha_tmp * homo * homo;
675  lumo = (1 - alpha_tmp + alpha_tmp * lumo) * (1 - alpha_tmp + alpha_tmp * lumo);
676  }
677  else
678  {
679  homo = (1 - alpha_tmp + alpha_tmp * homo) * (1 - alpha_tmp + alpha_tmp * homo);
680  lumo = 2 * alpha_tmp * lumo - alpha_tmp * alpha_tmp * lumo * lumo;
681  }
682 }
683 
684 
685 template<typename MatrixType>
688 {
689  assert(it > 0);
690 
691  real Df;
692  real temp, a, b;
693  int poly;
694  real alpha;
695 
696  a = x;
697  Df = 1;
698  DDf = -1; // TODO
699 
700  for (int i = 1; i <= it; i++)
701  {
702  temp = a;
703 
704  get_poly(i, poly, alpha);
705 
706  if (poly == 1)
707  {
708  a = ((1 - alpha) + alpha * temp) * ((1 - alpha) + alpha * temp);
709  b = 2 * alpha * ((1 - alpha) + alpha * temp);
710  }
711  else
712  {
713  a = 2 * alpha * temp - alpha * alpha * temp * temp;
714  b = 2 * alpha - 2 * alpha * alpha * temp;
715  }
716  Df *= b;
717  }
718 
719  return Df;
720 }
721 
722 
723 #endif //HEADER_PURIFICATION_SP2ACC
Treal upp() const
Definition: Interval.h:145
ergo_real real
Definition: purification_general.h:114
#define C_SP2
Constant used on the stopping criterion for the SP2 method.
Definition: constants.h:42
VectorTypeReal VecAlpha
Definition: purification_sp2acc.h:96
virtual void purify_X(const int it)
Definition: purification_sp2acc.h:183
virtual void apply_poly_to_eigs(const int it, real &homo, real &lumo)
Definition: purification_sp2acc.h:660
intervalType IntervalType
Definition: random_matrices.h:67
virtual void set_poly(const int it)
Definition: purification_sp2acc.h:109
PurificationGeneral is an abstract class which provides an interface for SP2, SP2ACC and possibly oth...
Definition: purification_general.h:111
generalVector VectorType
Definition: purification_sp2acc.h:62
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:426
real alpha
Definition: puri_info.h:91
virtual void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2acc.h:591
int check_stopping_criterion_iter
Iteration when to start to check stopping criterion.
Definition: purification_general.h:405
int poly
Definition: puri_info.h:78
#define LOG_AREA_DENSFROMF
Definition: output.h:61
virtual void set_init_params()
Definition: purification_sp2acc.h:66
std::vector< int > VectorTypeInt
Definition: purification_general.h:117
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2acc.h:56
virtual void return_constant_C(const int it, real &Cval)
Definition: purification_sp2acc.h:316
Treal template_blas_fabs(Treal x)
ergo_real real
Definition: test.cc:46
virtual void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2acc.h:576
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2acc.h:59
#define max(a, b)
Definition: integrator.cc:87
Purification_sp2acc is a class which provides an interface for SP2ACC recursive expansion.
Definition: purification_sp2acc.h:51
#define LOG_CAT_INFO
Definition: output.h:49
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
static const real deltaTurnOffAcc
Definition: purification_sp2acc.h:99
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2acc.h:55
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2acc.h:57
int method
Definition: puri_info.h:171
std::vector< real > VectorTypeReal
Definition: purification_general.h:118
Recursive density matrix expansion (or density matrix purification).
Definition: puri_info.h:52
virtual void truncate_matrix(real &threshold, int it)
Definition: purification_sp2acc.h:392
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:246
Definition: matInclude.h:139
virtual real apply_poly(const int it, real x)
Definition: purification_sp2acc.h:633
virtual real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2acc.h:687
Purification_sp2acc()
Definition: purification_sp2acc.h:64
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
real gap
Definition: puri_info.h:79
virtual void estimate_number_of_iterations(int &numit)
Definition: purification_sp2acc.h:439
virtual void get_poly(const int it, int &poly, real &alpha)
Definition: purification_sp2acc.h:169
virtual void purify_bounds(const int it)
Definition: purification_sp2acc.h:253
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2acc.h:60
Treal template_blas_sqrt(Treal x)
normType
Definition: matInclude.h:139