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 puri_info.cc
31
32 @brief File containing definitions of functions declared in classes IterationInfo and PuriInfo.
33 IterationInfo is a class with the information stored for each iteration of the recursive expansion.
34 PuriInfo is a class containing general information about the recursive expansion.
35
36 @author Anastasia Kruchinina <em>responsible</em>
37 @sa puri_info.h
38 */
39
40
41
42 #include "puri_info.h"
43
44
45 typedef PuriInfo::real real;
46
47
get_total_Xsquare_time()48 real PuriInfo::get_total_Xsquare_time()
49 {
50 real total_sum = 0;
51 for(int it = 0; it <= total_it; ++it)
52 total_sum += Iterations[it].Xsquare_time;
53 return total_sum;
54 }
55
get_total_Xtrunc_time()56 real PuriInfo::get_total_Xtrunc_time()
57 {
58 real total_sum = 0;
59 for(int it = 0; it <= total_it; ++it)
60 total_sum += Iterations[it].trunc_time;
61 return total_sum;
62 }
63
get_total_nnz_time()64 real PuriInfo::get_total_nnz_time()
65 {
66 real total_sum = 0;
67 for(int it = 0; it <= total_it; ++it)
68 total_sum += Iterations[it].nnz_time;
69 return total_sum;
70 }
71
72
get_total_inf_diff_time()73 real PuriInfo::get_total_inf_diff_time()
74 {
75 real total_sum = 0;
76 for(int it = 0; it <= total_it; ++it)
77 total_sum += Iterations[it].inf_diff_time;
78 return total_sum;
79 }
get_total_eucl_diff_time()80 real PuriInfo::get_total_eucl_diff_time()
81 {
82 real total_sum = 0;
83 for(int it = 0; it <= total_it; ++it)
84 total_sum += Iterations[it].eucl_diff_time;
85 return total_sum;
86 }
get_total_mixed_diff_time()87 real PuriInfo::get_total_mixed_diff_time()
88 {
89 real total_sum = 0;
90 for(int it = 0; it <= total_it; ++it)
91 total_sum += Iterations[it].mixed_diff_time;
92 return total_sum;
93 }
get_total_frob_diff_time()94 real PuriInfo::get_total_frob_diff_time()
95 {
96 real total_sum = 0;
97 for(int it = 0; it <= total_it; ++it)
98 total_sum += Iterations[it].frob_diff_time;
99 return total_sum;
100 }
101
get_total_stopping_criterion_time()102 real PuriInfo::get_total_stopping_criterion_time()
103 {
104 real total_sum = 0;
105 for(int it = 0; it <= total_it; ++it)
106 total_sum += Iterations[it].stopping_criterion_time;
107 return total_sum;
108 }
109
get_total_trace_diff_time()110 real PuriInfo::get_total_trace_diff_time()
111 {
112 real total_sum = 0;
113 for(int it = 0; it <= total_it; ++it)
114 total_sum += Iterations[it].trace_diff_time;
115 return total_sum;
116 }
117
get_total_purify_time()118 real PuriInfo::get_total_purify_time()
119 {
120 real total_sum = 0;
121 for(int it = 0; it <= total_it; ++it)
122 total_sum += Iterations[it].purify_time;
123 return total_sum;
124 }
125
126
127
get_poly_seq(std::vector<int> & poly)128 void PuriInfo::get_poly_seq(std::vector<int> &poly)
129 {
130 poly.clear();
131 poly.resize(total_it+1);
132
133 for(int it = 0; it <= total_it; ++it)
134 poly[it] = Iterations[it].poly;
135 }
136
get_vec_frob_norms(std::vector<real> & norms)137 void PuriInfo::get_vec_frob_norms(std::vector<real> &norms)
138 {
139 norms.clear();
140 norms.resize(total_it+1);
141
142 for(int it = 0; it <= total_it; ++it)
143 norms[it] = Iterations[it].XmX2_fro_norm;
144 }
145
get_vec_infty_norms(std::vector<real> & norms)146 void PuriInfo::get_vec_infty_norms(std::vector<real> &norms)
147 {
148 norms.clear();
149 norms.resize(total_it+1);
150
151 for(int it = 0; it <= total_it; ++it)
152 norms[it] = Iterations[it].XmX2_infty_norm;
153 }
154
get_vec_mixed_norms(std::vector<real> & norms)155 void PuriInfo::get_vec_mixed_norms(std::vector<real> &norms)
156 {
157 norms.clear();
158 norms.resize(total_it+1);
159
160 for(int it = 0; it <= total_it; ++it)
161 norms[it] = Iterations[it].XmX2_mixed_norm;
162 }
163
164
get_vec_traces(std::vector<real> & traces)165 void PuriInfo::get_vec_traces(std::vector<real> & traces)
166 {
167 traces.clear();
168 traces.resize(total_it+1);
169
170 for(int it = 0; it <= total_it; ++it)
171 traces[it] = Iterations[it].XmX2_trace;
172 }
173
174
get_spectrum_bounds(real & lower_spectrum_bound_,real & upper_spectrum_bound_) const175 void PuriInfo::get_spectrum_bounds(real &lower_spectrum_bound_, real &upper_spectrum_bound_) const
176 {
177 lower_spectrum_bound_ = lower_spectrum_bound;
178 upper_spectrum_bound_ = upper_spectrum_bound;
179 }
180
set_spectrum_bounds(const real lower_spectrum_bound_,const real upper_spectrum_bound_)181 void PuriInfo::set_spectrum_bounds(const real lower_spectrum_bound_, const real upper_spectrum_bound_)
182 {
183 lower_spectrum_bound = lower_spectrum_bound_;
184 upper_spectrum_bound = upper_spectrum_bound_;
185 }
186
print_collected_info()187 void PuriInfo::print_collected_info()
188 {
189
190 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"=========================================================================================================================================");
191 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF," SHORT OUTPUT ");
192 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"=========================================================================================================================================");
193
194 // if used new stopping criterion
195 if( stopping_criterion == 1)
196 {
197 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Used the new stopping criterion");
198 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Estimated number of iterations:\t %d\n", estim_total_it);
199
200 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Total number of iterations: iter + addit = %d + %d\n", total_it - additional_iterations, additional_iterations);
201 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Iter \t Poly \t Order Frob \t Trace \t\t NNZ(X) \t NNZ(X2) \t Treshold \t Total time");
202 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"------------------------------------------------------------------------------------------------------------------------");
203 for(int i = 1; i <= total_it - additional_iterations; i++)
204 {
205 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"%3d \t %1d \t %6.2lf\t %12.5e \t %12.5e \t %6.2lf \t %6.2lf \t %9.5e \t %9.5lf", Iterations[i].it, Iterations[i].poly, Iterations[i].order,
206 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, Iterations[i].NNZ_X, Iterations[i].NNZ_X2,
207 (double)Iterations[i].threshold_X, (double)Iterations[i].total_time);
208 }
209 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"-----------------------------------------------------------------");
210 for(int i = total_it - additional_iterations + 1; i <= total_it; i++)
211 {
212 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"%3d \t %1d \t %6.2lf\t %12.5e \t %12.5e \t %6.2lf \t %6.2lf \t %9.5e \t %9.5lf", Iterations[i].it, Iterations[i].poly, Iterations[i].order,
213 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, Iterations[i].NNZ_X, Iterations[i].NNZ_X2,
214 (double)Iterations[i].threshold_X, (double)Iterations[i].total_time);
215 }
216 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Total time of the purification: \t %lf sec\n", (double)total_time);
217
218 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"");
219 if(method == 2) // accelerated
220 {
221 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Iter \t Poly \t Alpha \t\t Lumo \t\t\t 1-Homo \t\t C ");
222 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"-----------------------------------------------------------------");
223 for(int i = 1; i <= total_it - additional_iterations; i++)
224 {
225 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"%d \t %d \t %lf \t [%lf, %lf] \t [%lf, %lf] \t %lf \n", Iterations[i].it, Iterations[i].poly, (double)Iterations[i].alpha, (double)Iterations[i].lumo_bound_low, (double)Iterations[i].lumo_bound_upp, (double)Iterations[i].homo_bound_low, (double)Iterations[i].homo_bound_upp, (double)Iterations[i].constantC);
226 }
227 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"-----------------------------------------------------------------");
228 for(int i = total_it - additional_iterations + 1; i <= total_it; i++)
229 {
230 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"%d \t %d \t %lf \t [%lf, %lf] \t [%lf, %lf] \t %lf \n", Iterations[i].it, Iterations[i].poly, (double)Iterations[i].alpha, (double)Iterations[i].lumo_bound_low, (double)Iterations[i].lumo_bound_upp, (double)Iterations[i].homo_bound_low, (double)Iterations[i].homo_bound_upp, (double)Iterations[i].constantC);
231 }
232 }
233
234
235 #ifdef CHECK_IF_STOPPED_TOO_LATE_OR_TOO_EARLY
236 double XmX2_fro_norm_at_stop = Iterations[total_it - additional_iterations ].XmX2_fro_norm;
237 double XmX2_fro_norm_before_stop = Iterations[total_it - additional_iterations - 2].XmX2_fro_norm;
238 double XmX2_fro_norm_after_stop = Iterations[total_it - additional_iterations + 2].XmX2_fro_norm;
239 if(XmX2_fro_norm_at_stop > XmX2_fro_norm_before_stop*2)
240 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"FAILED (STOPPED TOO LATE)! XmX2_fro_norm_before_stop = %e, XmX2_fro_norm_at_stop = %e.\n", XmX2_fro_norm_before_stop, XmX2_fro_norm_at_stop);
241 if(XmX2_fro_norm_at_stop > XmX2_fro_norm_after_stop*2 && XmX2_fro_norm_after_stop != 0.0)
242 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"FAILED (STOPPED TOO EARLY)! XmX2_fro_norm_after_stop = %e, XmX2_fro_norm_at_stop = %e.\n", XmX2_fro_norm_after_stop, XmX2_fro_norm_at_stop);
243 #endif
244 }
245 else // old stopping criterion
246 {
247 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Used the old stopping criterion");
248 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Estimated number of iterations:\t %d\n", estim_total_it);
249
250 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Total number of iterations: iter + addit = %d + %d\n", total_it - additional_iterations, additional_iterations);
251 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"");
252 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Iter \t Frob \t\t Trace \t\t NNZ(X) NNZ(X2) Total time ");
253 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"-----------------------------------------------------------------------------");
254 for(int i = 1; i <= total_it - additional_iterations; i++)
255 {
256 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"%d \t %e \t %e \t %d \t %d \t %lf \n", Iterations[i].it,
257 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, Iterations[i].NNZ_X, Iterations[i].NNZ_X2,
258 (double)Iterations[i].total_time);
259 }
260 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"");
261 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"Total time of the purification: \t %lf sec\n", (double)total_time);
262 }
263 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF,"=========================================================================================================================================");
264
265 }
266
267
268
print_collected_info_printf()269 void PuriInfo::print_collected_info_printf()
270 {
271
272 printf("=========================================================================================================================================\n");
273 printf("= SHORT OUTPUT =\n");
274 printf("=========================================================================================================================================\n");
275
276 // if used new stopping criterion
277 if( stopping_criterion == 1)
278 {
279 printf("Used the new stopping criterion\n");
280 printf("Estimated number of iterations:\t %d\n\n", estim_total_it);
281
282 printf("Total number of iterations: iter + addit = %d + %d\n", total_it - additional_iterations, additional_iterations);
283 printf("Iter \t Poly \t Order Frob \t Trace \t\t NNZ(X) \t NNZ(X2) \t Treshold \t Total time\n");
284 printf("-------------------------------------------------------------------------------------------------------------------------------------\n");
285 for(int i = 1; i <= total_it - additional_iterations; i++)
286 {
287 printf("%3d \t %1d \t %6.2lf\t %12.5e \t %12.5e \t %6.2lf \t %6.2lf \t %9.5e \t %9.5lf\n", Iterations[i].it,
288 Iterations[i].poly, (double)Iterations[i].order,
289 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, (double)Iterations[i].NNZ_X, (double)Iterations[i].NNZ_X2,
290 (double)Iterations[i].threshold_X, (double)Iterations[i].total_time);
291 }
292 printf("-------------------------------------------------------------------------------------------------------------------------------------\n");
293 for(int i = total_it - additional_iterations + 1; i <= total_it; i++)
294 {
295 printf("%3d \t %1d \t %6.2lf\t %12.5e \t %12.5e \t %6.2lf \t %6.2lf \t %9.5e \t %9.5lf\n", Iterations[i].it,
296 Iterations[i].poly, (double)Iterations[i].order,
297 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, (double)Iterations[i].NNZ_X, (double)Iterations[i].NNZ_X2,
298 (double)Iterations[i].threshold_X, (double)Iterations[i].total_time);
299 }
300 printf("Total time of the purification: \t %lf sec\n", (double)total_time);
301
302 printf("\n");
303 if(method == 2) // accelerated
304 {
305 printf("Iter \t Poly \t Alpha \t\t Lumo \t\t\t 1-Homo \t\t C \n");
306 printf("-----------------------------------------------------------------\n");
307 for(int i = 1; i <= total_it - additional_iterations; i++)
308 {
309 printf("%d \t %d \t %lf \t [%lf, %lf] \t [%lf, %lf] \t %lf \n", Iterations[i].it, Iterations[i].poly,
310 (double)Iterations[i].alpha, (double)Iterations[i].lumo_bound_low, (double)Iterations[i].lumo_bound_upp,
311 (double)Iterations[i].homo_bound_low, (double)Iterations[i].homo_bound_upp, (double)Iterations[i].constantC);
312 }
313 printf("-------------------------------------------------------------------------------------------------------------------------------------\n");
314 for(int i = total_it - additional_iterations + 1; i <= total_it; i++)
315 {
316 printf("%d \t %d \t %lf \t [%lf, %lf] \t [%lf, %lf] \t %lf \n", Iterations[i].it, Iterations[i].poly,
317 (double)Iterations[i].alpha, (double)Iterations[i].lumo_bound_low, (double)Iterations[i].lumo_bound_upp,
318 (double)Iterations[i].homo_bound_low, (double)Iterations[i].homo_bound_upp, (double)Iterations[i].constantC);
319 }
320 }
321
322
323 #ifdef CHECK_IF_STOPPED_TOO_LATE_OR_TOO_EARLY
324 double XmX2_fro_norm_at_stop = Iterations[total_it - additional_iterations ].XmX2_fro_norm;
325 double XmX2_fro_norm_before_stop = Iterations[total_it - additional_iterations - 2].XmX2_fro_norm;
326 double XmX2_fro_norm_after_stop = Iterations[total_it - additional_iterations + 2].XmX2_fro_norm;
327 if(XmX2_fro_norm_at_stop > XmX2_fro_norm_before_stop*2)
328 printf("FAILED (STOPPED TOO LATE)! XmX2_fro_norm_before_stop = %e, XmX2_fro_norm_at_stop = %e.\n", XmX2_fro_norm_before_stop, XmX2_fro_norm_at_stop);
329 if(XmX2_fro_norm_at_stop > XmX2_fro_norm_after_stop*2 && XmX2_fro_norm_after_stop != 0.0)
330 printf("FAILED (STOPPED TOO EARLY)! XmX2_fro_norm_after_stop = %e, XmX2_fro_norm_at_stop = %e.\n", XmX2_fro_norm_after_stop, XmX2_fro_norm_at_stop);
331 #endif
332 }
333 else // old stopping criterion
334 {
335 printf("Used the old stopping criterion\n");
336 printf("Estimated number of iterations:\t %d\n", estim_total_it);
337
338 printf("Total number of iterations: iter + addit = %d + %d\n", total_it - additional_iterations, additional_iterations);
339 printf("\n");
340 printf("Iter \t Frob \t\t Trace \t\t NNZ(X) NNZ(X2) Total time \n");
341 printf("-----------------------------------------------------------------------------\n");
342 for(int i = 1; i <= total_it - additional_iterations; i++)
343 {
344 printf("%d \t %e \t %e \t %.2lf \t %.2lf \t %lf \n", Iterations[i].it,
345 (double)Iterations[i].XmX2_fro_norm, (double)Iterations[i].XmX2_trace, (double)Iterations[i].NNZ_X, (double)Iterations[i].NNZ_X2,
346 (double)Iterations[i].total_time);
347 }
348 printf("\n");
349 printf("Total time of the purification: \t %lf sec\n", (double)total_time);
350 }
351 printf("=========================================================================================================================================\n");
352
353 }
354
355