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