ergo
purification_sp2.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_SP2
40 #define HEADER_PURIFICATION_SP2
41 
42 #include "purification_general.h"
43 
44 //#define DEBUG_OUTPUT
45 
50 template<typename MatrixType>
51 class Purification_sp2 : public PurificationGeneral<MatrixType>
52 {
53 public:
54 
58 
61 
63 
65 
67  {
68  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen SP2 purification method");
69  this->info.method = 1;
70 
71  this->gammaStopEstim = (3 - template_blas_sqrt((real)5)) / 2.0;
72  }
73 
74  void get_poly(const int it, int& poly);
75  void set_poly(const int it);
76 
77  void estimate_number_of_iterations(int& numit);
78  void purify_X(const int it);
79  void purify_bounds(const int it);
80  void save_other_iter_info(IterationInfo& iter_info, int it);
81  void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it);
82 
83  void return_constant_C(const int it, real& Cval);
84 
85  //real apply_inverse_poly(const int it, real x);
86  real apply_poly(const int it, real x);
87  real compute_derivative(const int it, real x, real& DDf);
88 };
89 
90 template<typename MatrixType>
92 {
93  assert((int)this->VecPoly.size() > it);
94 
95  // if cannot compute polynomial using homo and lumo eigevalues, compute using trace
96  if (this->VecPoly[it] == -1)
97  {
98  real Xtrace = this->X.trace();
99  real Xsqtrace = this->Xsq.trace();
100 
101  if ((template_blas_fabs(Xsqtrace - this->nocc) <
102  template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
103  ||
104  (it % 2
105  &&
106  (template_blas_fabs(Xsqtrace - this->nocc) ==
107  template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
108  ))
109  {
110  this->VecPoly[it] = 1;
111  }
112  else
113  {
114  this->VecPoly[it] = 0;
115  }
116  }
117 }
118 
119 
120 template<typename MatrixType>
121 void Purification_sp2<MatrixType>::get_poly(const int it, int& poly)
122 {
123  assert((int)this->VecPoly.size() > it);
124  assert(this->VecPoly[it] != -1);
125  poly = this->VecPoly[it];
126 }
127 
128 
129 template<typename MatrixType>
131 {
132 #ifdef DEBUG_OUTPUT
133  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Purify X...");
134 #endif
135  int poly;
136 
137  set_poly(it);
138  get_poly(it, poly);
139 
140  /* It may happen that X2 has many more nonzeros than X, for
141  * example 5 times as many. Therefore it makes sense to try
142  * having only one "big" matrix in memory at a time. However,
143  * file operations have proved to be quite expensive and should
144  * be avoided if possible. Hence we want to achieve having only
145  * one big matrix in memory without unnecessary file
146  * operations. We are currently hoping that it will be ok to add
147  * a "small" matrix to a "big" one, that the memory usage after
148  * that operation will be like the memory usage for one big
149  * matrix + one small matrix. Therefore we are adding X to X2 (X
150  * is truncated, a "small" matrix) instead of the opposite when
151  * the 2*X-X*X polynomial is evaluated.
152  */
153 
154  if (poly == 0)
155  {
156  this->Xsq *= ((real) - 1.0);
157  this->X *= (real)2.0;
158  this->Xsq += this->X; // Xsq = -Xsq + 2X
159  }
160 
161  this->Xsq.transfer(this->X); // clear Xsq and old X
162 }
163 
164 
165 template<typename MatrixType>
167 {
168 #ifdef DEBUG_OUTPUT
169  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change homo and lumo bounds according to the chosen polynomial VecPoly = %d", this->VecPoly[it]);
170 #endif
171  real homo_low, homo_upp, lumo_upp, lumo_low;
172  int poly;
173 
174  get_poly(it, poly);
175 
176  if (poly == 1)
177  {
178  // update bounds
179  homo_low = 2 * this->homo_bounds.low() - this->homo_bounds.low() * this->homo_bounds.low(); // 2x - x^2
180  homo_upp = 2 * this->homo_bounds.upp() - this->homo_bounds.upp() * this->homo_bounds.upp(); // 2x - x^2
181  lumo_low = this->lumo_bounds.low() * this->lumo_bounds.low(); // x^2
182  lumo_upp = this->lumo_bounds.upp() * this->lumo_bounds.upp(); // x^2
183  this->homo_bounds = IntervalType(homo_low, homo_upp);
184  this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
185  }
186  else
187  {
188  // update bounds
189  lumo_low = 2 * this->lumo_bounds.low() - this->lumo_bounds.low() * this->lumo_bounds.low(); // 2x - x^2
190  lumo_upp = 2 * this->lumo_bounds.upp() - this->lumo_bounds.upp() * this->lumo_bounds.upp(); // 2x - x^2
191  homo_low = this->homo_bounds.low() * this->homo_bounds.low(); // x^2
192  homo_upp = this->homo_bounds.upp() * this->homo_bounds.upp(); // x^2
193  this->homo_bounds = IntervalType(homo_low, homo_upp);
194  this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
195  }
196 
197  IntervalType zero_one(0, 1);
198  this->homo_bounds.intersect(zero_one);
199  this->lumo_bounds.intersect(zero_one);
200 
201 #ifdef DEBUG_OUTPUT
202  if (this->homo_bounds.empty())
203  {
204  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
205  }
206  if (this->lumo_bounds.empty())
207  {
208  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
209  }
210 
211  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "1-homo: [ %lf , %lf ],", this->homo_bounds.low(), this->homo_bounds.upp());
212  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo: [ %lf , %lf ].", this->lumo_bounds.low(), this->lumo_bounds.upp());
213 #endif
214 }
215 
216 
217 /**************************************************************************************/
218 
219 template<typename MatrixType>
221 {
222  Cval = C_SP2;
223 }
224 
225 
226 /**************************************************************************************/
227 
228 
229 template<typename MatrixType>
231 {
232  int it = 1;
233  int maxit_total = this->maxit + this->additional_iterations;
234  int maxit_tmp = maxit_total;
235  real x, y;
236  real epsilon = this->get_epsilon();
237 
238  // we need values on iteration it-2 for the stopping criterion
239  this->check_stopping_criterion_iter = 2;
240 
241  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
242  this->VecPoly.clear();
243  this->VecPoly.resize(max_size, -1);
244 
245  this->VecGap.clear();
246  this->VecGap.resize(max_size, -1);
247 
248  // we are interested in the inner bounds of gap
249  x = this->lumo_bounds.upp(); // = lumo
250  y = this->homo_bounds.upp(); // = 1 - homo
251 
252 
253  // if VecGap is zero
254  if (1 - x - y <= 0)
255  {
256 #ifdef DEBUG_OUTPUT
257  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap cannot be computed. Set estimated number of iteration to the maxit.");
258 #endif
259  estim_num_iter = this->maxit;
260  return;
261  }
262 
263  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT LUMO: [ %.12lf , %.12lf ]", (double)this->lumo_bounds.low(), (double)this->lumo_bounds.upp());
264  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT HOMO: [ %.12lf , %.12lf ]", (double)this->homo_bounds.low(), (double)this->homo_bounds.upp());
265 
266 
267 
268  this->VecPoly[0] = -1;
269  this->VecGap[0] = 1 - x - y;
270 
271  estim_num_iter = -1;
272 
273  while (it <= maxit_tmp)
274  {
275  // note: avoid not-stopping in case of idempotent matrix
276  if ((x > y) || (it % 2 && (x == y))) // lumo > 1-homo
277  {
278  x *= x;
279  y = 2 * y - y * y;
280  this->VecPoly[it] = 1;
281  }
282  else
283  {
284  x = 2 * x - x * x;
285  y *= y;
286  this->VecPoly[it] = 0;
287  }
288 
289  this->VecGap[it] = 1 - x - y;
290 
291  // maybe we wish to perform some more iterations, then stopping criterion suggest
292  if ((estim_num_iter == -1) &&
293  (x - x * x < epsilon) && (y - y * y < epsilon)
294  // the eucledian norm is less then epsilon
295  &&
296  (this->VecPoly[it] != this->VecPoly[it - 1])) // to apply stopping criterion, polynomials must alternate
297  {
298  estim_num_iter = it;
299  maxit_tmp = it + this->additional_iterations;
300 
301  // 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)
302  if (this->normPuriStopCrit == mat::frobNorm)
303  {
304  estim_num_iter += 2;
305  maxit_tmp += 2;
306  }
307  }
308 
309  ++it;
310  } //while
311 
312 
313  /*
314  Either we reached maxit number of iterations or due to the additional 2 iterations (because of the Frobenius norm) we overestimated number of iterations
315  */
316  if ( ((estim_num_iter == -1) && (it == maxit_tmp + 1) )
317  || (estim_num_iter > maxit_total))
318  {
319  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxit = %d number of iterations is reached in estimate_number_of_iteration()", this->maxit);
320  estim_num_iter = this->maxit;
321  maxit_tmp = maxit_total;
322  }
323 
324 
325  this->VecPoly.resize(maxit_tmp + 1); // make it less if needed
326  this->VecGap.resize(maxit_tmp + 1);
327 
328 #ifdef DEBUG_OUTPUT
329  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Sequence of polynomials VecPoly: ");
330  for (int i = 0; i < (int)this->VecPoly.size(); ++i)
331  {
332  do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "%d ", this->VecPoly[i]);
333  }
335 #endif
336 }
337 
338 
339 /************ SAVE INFORMATION ABOUT ITERATIONS SPECIFIC FOR SP2 PURIFICATION *************/
340 
341 template<typename MatrixType>
343 {
344  assert((int)this->VecPoly.size() > it);
345  assert((int)this->VecGap.size() > it);
346 
347  iter_info.poly = this->VecPoly[it];
348  iter_info.gap = this->VecGap[it];
349 }
350 
351 
352 /************ APPLY INVERSE POLYNOMIAL (FOR ESTIMATION OF EIGENVALUES) ************/
353 
354 
355 template<typename MatrixType>
357 {
358  int poly;
359  for (int i = it; i >= 1; i--)
360  {
361  get_poly(i, poly);
362 
363  if (poly == 1)
364  {
365  bounds_from_it[0] = template_blas_sqrt(bounds_from_it[0]);
366  bounds_from_it[1] = template_blas_sqrt(bounds_from_it[1]);
367 
368  bounds_from_it[2] = bounds_from_it[2] / (1 + template_blas_sqrt(1 - bounds_from_it[2]));
369  bounds_from_it[3] = bounds_from_it[3] / (1 + template_blas_sqrt(1 - bounds_from_it[3]));
370  }
371  else
372  {
373  bounds_from_it[0] = bounds_from_it[0] / (1 + template_blas_sqrt(1 - bounds_from_it[0]));
374  bounds_from_it[1] = bounds_from_it[1] / (1 + template_blas_sqrt(1 - bounds_from_it[1]));
375 
376  bounds_from_it[2] = template_blas_sqrt(bounds_from_it[2]);
377  bounds_from_it[3] = template_blas_sqrt(bounds_from_it[3]);
378  }
379  }
380 }
381 
382 
383 /*
384  *
385  *
386  * template<typename MatrixType>
387  * typename Purification_sp2<MatrixType>::real
388  * Purification_sp2<MatrixType>::apply_inverse_poly(const int it, real x)
389  * {
390  * if( it == 0 ) return -1; // no polynomials was applied in 0 iteration
391  *
392  * int poly;
393  *
394  * real finvx = x;
395  * for(int i = it; i >= 1; i--)
396  * {
397  * get_poly(i, poly);
398  *
399  * if(poly == 1)
400  * finvx = template_blas_sqrt(finvx);
401  * else
402  * finvx = finvx/(1+template_blas_sqrt(1-finvx));
403  * }
404  * return finvx;
405  * }
406  */
407 
408 template<typename MatrixType>
411 {
412  assert(it >= 0);
413  if (it == 0)
414  {
415  return x;
416  }
417 
418  real fx;
419  int poly;
420  get_poly(it, poly);
421 
422  if (poly == 1)
423  {
424  fx = x * x;
425  }
426  else
427  {
428  fx = 2 * x - x * x;
429  }
430 
431  return fx;
432 }
433 
434 
435 
436 
437 template<typename MatrixType>
440 {
441  assert(it > 0);
442 
443  real Df;
444  real a;
445  int poly;
446 
447  a = x;
448  Df = 1;
449  DDf = 1;
450 
451  for (int i = 1; i <= it; i++)
452  {
453  get_poly(i, poly);
454 
455  if (poly == 1)
456  {
457  DDf = 2 * Df * Df + (2 * a) * DDf;
458  Df *= 2 * a;
459  a = a * a;
460  }
461  else
462  {
463  DDf = -2 * Df * Df + (2 - 2 * a) * DDf;
464  Df *= 2 - 2 * a;
465  a = 2 * a - a * a;
466  }
467  }
468 
469  //do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Derivative (it = %d) = %lf\n", it, Df);
470 
471  return Df;
472 }
473 
474 
475 #endif //HEADER_PURIFICATION_SP2
Purification_sp2::apply_poly
real apply_poly(const int it, real x)
Definition: purification_sp2.h:410
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
Purification_sp2::save_other_iter_info
void save_other_iter_info(IterationInfo &iter_info, int it)
Definition: purification_sp2.h:342
Purification_sp2::purify_X
void purify_X(const int it)
Definition: purification_sp2.h:130
PurificationGeneral::info
PuriInfo info
Fill in during purification with useful information.
Definition: purification_general.h:128
Purification_sp2::Purification_sp2
Purification_sp2()
Definition: purification_sp2.h:64
PurificationGeneral::gammaStopEstim
real gammaStopEstim
Used on the stopping criterion for estimation of eigenvalues from purification.
Definition: purification_general.h:503
Purification_sp2::VectorTypeInt
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition: purification_sp2.h:59
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_sp2::get_poly
void get_poly(const int it, int &poly)
Definition: purification_sp2.h:121
Purification_sp2::apply_inverse_poly_vector
void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition: purification_sp2.h:356
mat::frobNorm
@ frobNorm
Definition: matInclude.h:139
real
ergo_real real
Definition: test.cc:46
template_blas_fabs
Treal template_blas_fabs(Treal x)
Purification_sp2::return_constant_C
void return_constant_C(const int it, real &Cval)
Definition: purification_sp2.h:220
Purification_sp2::VectorTypeReal
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition: purification_sp2.h:60
IterationInfo::poly
int poly
Definition: puri_info.h:81
Purification_sp2::compute_derivative
real compute_derivative(const int it, real x, real &DDf)
Definition: purification_sp2.h:439
LOG_AREA_DENSFROMF
#define LOG_AREA_DENSFROMF
Definition: output.h:61
Purification_sp2
Purification_sp2acc is a class which provides an interface for SP2 recursive expansion.
Definition: purification_sp2.h:52
Purification_sp2::NormType
PurificationGeneral< MatrixType >::NormType NormType
Definition: purification_sp2.h:57
IterationInfo
Definition: puri_info.h:52
MatrixType
symmMatrix MatrixType
Definition: recexp_diag_test.cc:66
mat::VectorGeneral
Definition: MatrixBase.h:55
PurificationGeneral::VectorTypeInt
std::vector< int > VectorTypeInt
Definition: purification_general.h:122
LOG_CAT_INFO
#define LOG_CAT_INFO
Definition: output.h:49
mat::normType
normType
Definition: matInclude.h:139
mat::Interval< ergo_real >
Purification_sp2::set_poly
void set_poly(const int it)
Definition: purification_sp2.h:91
PurificationGeneral::VectorTypeReal
std::vector< real > VectorTypeReal
Definition: purification_general.h:123
IterationInfo::gap
real gap
Definition: puri_info.h:82
Purification_sp2::real
PurificationGeneral< MatrixType >::real real
Definition: purification_sp2.h:55
Purification_sp2::VectorType
generalVector VectorType
Definition: purification_sp2.h:62
Purification_sp2::set_init_params
void set_init_params()
Definition: purification_sp2.h:66
Purification_sp2::estimate_number_of_iterations
void estimate_number_of_iterations(int &numit)
Definition: purification_sp2.h:230
do_output
void do_output(int logCategory, int logArea, const char *format,...)
Definition: output.cc:53
Purification_sp2::IntervalType
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition: purification_sp2.h:56
Purification_sp2::purify_bounds
void purify_bounds(const int it)
Definition: purification_sp2.h:166