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
30 /** @file purification_sp2.h
31
32 @brief SP2 recursive density matrix expansion (or density matrix
33 purification).
34
35 @author Anastasia Kruchinina <em>responsible</em>
36 */
37
38
39 #ifndef HEADER_PURIFICATION_SP2
40 #define HEADER_PURIFICATION_SP2
41
42 #include "purification_general.h"
43
44 //#define DEBUG_OUTPUT
45
46 /** Purification_sp2acc is a class which provides an interface for SP2
47 * recursive expansion.
48 *
49 * \tparam MatrixType Type of a matrix (ex. symmMatrix). */
50 template<typename MatrixType>
51 class Purification_sp2 : public PurificationGeneral<MatrixType>
52 {
53 public:
54
55 typedef typename PurificationGeneral<MatrixType>::real real;
56 typedef typename PurificationGeneral<MatrixType>::IntervalType IntervalType;
57 typedef typename PurificationGeneral<MatrixType>::NormType NormType;
58
59 typedef typename PurificationGeneral<MatrixType>::VectorTypeInt VectorTypeInt;
60 typedef typename PurificationGeneral<MatrixType>::VectorTypeReal VectorTypeReal;
61
62 typedef generalVector VectorType;
63
Purification_sp2()64 Purification_sp2() : PurificationGeneral<MatrixType>() {}
65
set_init_params()66 void set_init_params()
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>
set_poly(const int it)91 void Purification_sp2<MatrixType>::set_poly(const int it)
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>
get_poly(const int it,int & poly)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>
purify_X(const int it)130 void Purification_sp2<MatrixType>::purify_X(const int it)
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>
purify_bounds(const int it)166 void Purification_sp2<MatrixType>::purify_bounds(const int it)
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>
return_constant_C(const int it,real & Cval)220 void Purification_sp2<MatrixType>::return_constant_C(const int it, real& Cval)
221 {
222 Cval = C_SP2;
223 }
224
225
226 /**************************************************************************************/
227
228
229 template<typename MatrixType>
estimate_number_of_iterations(int & estim_num_iter)230 void Purification_sp2<MatrixType>::estimate_number_of_iterations(int& estim_num_iter)
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 }
334 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "");
335 #endif
336 }
337
338
339 /************ SAVE INFORMATION ABOUT ITERATIONS SPECIFIC FOR SP2 PURIFICATION *************/
340
341 template<typename MatrixType>
save_other_iter_info(IterationInfo & iter_info,int it)342 void Purification_sp2<MatrixType>::save_other_iter_info(IterationInfo& iter_info, int it)
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>
apply_inverse_poly_vector(const int it,VectorTypeReal & bounds_from_it)356 void Purification_sp2<MatrixType>::apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it)
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>
409 typename Purification_sp2<MatrixType>::real
apply_poly(const int it,real x)410 Purification_sp2<MatrixType>::apply_poly(const int it, real x)
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>
438 typename Purification_sp2<MatrixType>::real
compute_derivative(const int it,real x,real & DDf)439 Purification_sp2<MatrixType>::compute_derivative(const int it, real x, real& DDf)
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
476