ergo
purification_sp2.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_SP2
40#define HEADER_PURIFICATION_SP2
41
43
44//#define DEBUG_OUTPUT
45
50template<typename MatrixType>
51class Purification_sp2 : public PurificationGeneral<MatrixType>
52{
53public:
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
90template<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
120template<typename MatrixType>
121void 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
129template<typename MatrixType>
131{
132#ifdef DEBUG_OUTPUT
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
165template<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
219template<typename MatrixType>
221{
222 Cval = C_SP2;
223}
224
225
226/**************************************************************************************/
227
228
229template<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
341template<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
355template<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
408template<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
437template<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
Definition puri_info.h:52
real gap
Definition puri_info.h:82
int poly
Definition puri_info.h:81
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
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 SP2 recursive expansion.
Definition purification_sp2.h:52
void set_init_params()
Definition purification_sp2.h:66
void apply_inverse_poly_vector(const int it, VectorTypeReal &bounds_from_it)
Definition purification_sp2.h:356
real apply_poly(const int it, real x)
Definition purification_sp2.h:410
real compute_derivative(const int it, real x, real &DDf)
Definition purification_sp2.h:439
PurificationGeneral< MatrixType >::NormType NormType
Definition purification_sp2.h:57
void return_constant_C(const int it, real &Cval)
Definition purification_sp2.h:220
PurificationGeneral< MatrixType >::real real
Definition purification_sp2.h:55
PurificationGeneral< MatrixType >::IntervalType IntervalType
Definition purification_sp2.h:56
void save_other_iter_info(IterationInfo &iter_info, int it)
Definition purification_sp2.h:342
PurificationGeneral< MatrixType >::VectorTypeReal VectorTypeReal
Definition purification_sp2.h:60
generalVector VectorType
Definition purification_sp2.h:62
void get_poly(const int it, int &poly)
Definition purification_sp2.h:121
void set_poly(const int it)
Definition purification_sp2.h:91
Purification_sp2()
Definition purification_sp2.h:64
void estimate_number_of_iterations(int &numit)
Definition purification_sp2.h:230
void purify_X(const int it)
Definition purification_sp2.h:130
PurificationGeneral< MatrixType >::VectorTypeInt VectorTypeInt
Definition purification_sp2.h:59
void purify_bounds(const int it)
Definition purification_sp2.h:166
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)