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 recexp_sprandsym.cc
31
32 \brief Test serial recursive expansion and computation of homo and lumo
33 eigenvectors of a sparse matrix with a given eigenspectrum. Matrix is generated using Givens rotations starting from a diagonal matrix with a required eigenvalues, so eigenvectors of a result matrix are explicitly known and saved into the matrix Q.
34
35 @author Anastasia Kruchinina <em>responsible</em>
36 */
37
38
39 #ifndef USE_CHUNKS_AND_TASKS
40
41 #include "purification_sp2.h"
42 #include "purification_sp2acc.h"
43 #include "matrix_typedefs.h" // definitions of matrix types and interval type (source)
44 #include "realtype.h" // definitions of types (utilities_basic)
45 #include "matrix_utilities.h"
46 #include "integral_matrix_wrappers.h"
47 #include "SizesAndBlocks.h"
48 #include "Matrix.h"
49 #include "Vector.h"
50 #include "MatrixSymmetric.h"
51 #include "MatrixTriangular.h"
52 #include "MatrixGeneral.h"
53 #include "VectorGeneral.h"
54 #include "output.h"
55
56 #include <iostream>
57 #include <fstream>
58 #include <sstream>
59 #include <string.h>
60
61 #include "random_matrices.h"
62 #include "get_eigenvectors.h"
63
64
65
66 typedef ergo_real real;
67
68 #define SQRT_EPSILON_REAL template_blas_sqrt(mat::getMachineEpsilon<real>())
69
70 real TOL_ERR_SUBS_DEFAULT = SQRT_EPSILON_REAL;
71 real TOL_TRACE_ERROR_DEFAULT = SQRT_EPSILON_REAL;
72 real SCALAR_TOL = SQRT_EPSILON_REAL;
73
74 #ifdef PRECISION_SINGLE
75 real TOL_EIGENSOLVER_ACC_DEFAULT = 1e-6;
76 #elif PRECISION_DOUBLE
77 real TOL_EIGENSOLVER_ACC_DEFAULT = 1e-12;
78 #elif PRECISION_LONG_DOUBLE
79 real TOL_EIGENSOLVER_ACC_DEFAULT = 1e-16;
80 #elif PRECISION_QUAD_FLT128
81 real TOL_EIGENSOLVER_ACC_DEFAULT = 1e-24;
82 #endif
83
84
abs_compare(real a,real b)85 static bool abs_compare(real a, real b)
86 {return (template_blas_fabs(a) < template_blas_fabs(b));}
87
88 typedef ergo_real real;
89 typedef symmMatrix MatrixType;
90 typedef normalMatrix MatrixTypeGeneral;
91 typedef MatrixType::VectorType VectorType;
92
93
main(int argc,char * argv[])94 int main(int argc, char *argv[])
95 {
96 printf("Program performing the recursive expansion on a given matrix.\n");
97 printf("Written by Anastasia Kruchinina, Nov 2017\n");
98 printf("Usage: %s N N_occ rand_seed gap gap_around sparsity err_sub [option filename]\n", argv[0]);
99 printf(" where option is: \n");
100 printf(" 1 - equidistant eigenvalues outside gap (default) in [0,1] \n");
101 printf(" 2 - random spectrum in [0,1]\n");
102 printf(" 3 - read vector with eigs in [a,b] from file 'filename'\n");
103 printf("Note: matrix sparsity is given in %%.\n");
104
105 printf("\n");
106
107 #ifdef _OPENMP
108 int defThreads;
109 const char *env = getenv("OMP_NUM_THREADS");
110 if ( !(env && (defThreads=atoi(env)) > 0) ) {
111 defThreads = 1;
112 }
113
114 mat::Params::setNProcs(defThreads);
115 mat::Params::setMatrixParallelLevel(2);
116 std::cout<<"OpenMP is used, number of threads set to "
117 <<mat::Params::getNProcs()<<". Matrix parallel level: "
118 <<mat::Params::getMatrixParallelLevel()<<"."<<std::endl;
119 #endif
120
121
122 int N = 200;
123 int N_occ = 50;
124 int rand_seed = 10000;
125 double gap = 0.01;
126 double gap_around = 0.5;
127 double MATRIX_SPARSITY = 20;
128 double err_sub = TOL_ERR_SUBS_DEFAULT;
129 int eig_option = 1;
130
131 real homo_in = MAX_DOUBLE;
132 real homo_out = -MAX_DOUBLE;
133 real lumo_in = -MAX_DOUBLE;
134 real lumo_out = MAX_DOUBLE;
135
136 char filename[100];
137
138
139 if(argc == 9 || argc == 10)
140 {
141 N = atoi(argv[1]);
142 N_occ = atoi(argv[2]);
143 rand_seed = atoi(argv[3]);
144 gap = atof(argv[4]);
145 gap_around = atof(argv[5]);
146 MATRIX_SPARSITY = atof(argv[6]);
147 err_sub = atof(argv[7]);
148 if(argc >= 9)
149 eig_option = atoi(argv[8]);
150 else
151 eig_option = 1;
152
153 if(argc == 10)
154 strcpy(filename, argv[9]);
155
156 if(eig_option == 3 && argc < 10)
157 {
158 printf("Filename is not specified!\n");
159 return EXIT_FAILURE;
160 }
161 }
162
163 printf("Data for the matrix F:\n");
164 printf("N = %4d\n", N );
165 printf("N_occ = %4d\n", N_occ );
166 printf("rand_seed = %4d\n", rand_seed);
167 printf("gap = %lf\n", gap);
168 printf("gap_around %lf\n", gap_around);
169 printf("err_sub = %e\n", (double)err_sub);
170 printf("sparsity = %lf\n", MATRIX_SPARSITY);
171 if(eig_option == 1)
172 printf("Equidistant eigenvalues in [0,1] outside gap\n");
173 else if(eig_option == 2)
174 printf("Random spectrum in [0,1]\n");
175 else if(eig_option == 3)
176 printf("Read vector with eigenvalues from file %s\n", filename);
177 else
178 printf("Random spectrum\n");
179
180
181 srand(rand_seed);
182
183 // Create random symmetric matrix F with eigenvalues eigvalList
184 std::vector<real> eigvalList(N);
185
186 if(eig_option == 1)
187 {
188 // [0, gap_around-gap/2]
189 for(int i = 1; i < N_occ; ++i)
190 eigvalList[i] = (real)i/(N_occ-1)*(gap_around-gap/2);
191
192 for(int i = 0; i < N-N_occ; ++i)
193 eigvalList[N_occ+i] = (real)i/(N-N_occ-1)*(1-(gap_around+gap/2)) + gap_around+gap/2;
194
195 }
196
197 if(eig_option == 2)
198 {
199 // [0, gap_around-gap/2]
200 for(int i = 1; i < N_occ; ++i)
201 eigvalList[i] = (real)rand()/RAND_MAX*(gap_around-gap/2);
202 eigvalList[N_occ-1] = gap_around-gap/2;
203
204 eigvalList[N_occ] = gap_around+gap/2;
205 for(int i = 1; i < N-N_occ-1; ++i)
206 eigvalList[N_occ+i] = (real)rand()/RAND_MAX*(1-(gap_around+gap/2)) + gap_around+gap/2;
207 eigvalList[N-1] = 1;
208
209 sort(eigvalList.begin(), eigvalList.end());
210
211 write_vector("vector_eigs.txt", eigvalList);
212 }
213
214
215 if(eig_option == 3)
216 {
217 if(read_vector(filename, eigvalList, N, false) == -1)
218 return EXIT_FAILURE;
219
220 gap = eigvalList[N_occ] - eigvalList[N_occ-1];
221 gap_around = (eigvalList[N_occ] + eigvalList[N_occ-1])/2;
222 }
223
224
225
226 // for(int i = 0; i < N; ++i)
227 // printf("%lf ", (double)eigvalList[i]);
228
229
230 printf("Generating matrix...\n");
231
232 MatrixType F;
233 MatrixTypeGeneral Q;
234 sprandsym(N, F, Q, eigvalList, MATRIX_SPARSITY);
235
236 printf("Matrices generated:\n");
237 printf(" F - Fock matrix, Q - matrix of eigenvectors of F\n");
238 printf("NNZ in F: %lu, sparsity is %lf %%\n", F.nnz(), (double)F.nnz()/(N*N) * 100);
239
240 // // write F to the file
241 // char filenameF[] = "F.bin";
242 // std::vector<real> vals;
243 // F.fullMatrix(vals);
244 // ofstream f (filenameF, ios::out | ios::binary);
245 // f.write ((char*)&vals[0], vals.size() * sizeof(real));
246 // f.close();
247 //
248 // // write Q to the file
249 // char filenameQ[] = "Q.bin";
250 // vals.clear();
251 // Q.fullMatrix(vals);
252 // ofstream fQ (filenameQ, ios::out | ios::binary);
253 // fQ.write ((char*)&vals[0], vals.size() * sizeof(real));
254 // fQ.close();
255 //
256 // printf("Matrices F and Q are saved into files.\n");
257
258 printf("Use spectral norm for the truncation and the stopping criterion.\n");
259 mat::normType normPuri = mat::euclNorm;
260 mat::normType normPuriStopCrit = mat::euclNorm;
261
262 printf("\n\n");
263
264
265 ergo_real homo = eigvalList[N_occ-1];
266 ergo_real lumo = eigvalList[N_occ ];
267
268 ergo_real epsilon_for_homo_lumo_intervals = 1e-4;
269 homo_out = homo-epsilon_for_homo_lumo_intervals;
270 homo_in = homo+epsilon_for_homo_lumo_intervals;
271 lumo_in = lumo-epsilon_for_homo_lumo_intervals;
272 lumo_out = lumo+epsilon_for_homo_lumo_intervals;
273
274 // SET HOMO AND LUMO BOUNDS
275 IntervalType homo_bounds(homo_out, homo_in);
276 IntervalType lumo_bounds(lumo_in, lumo_out);
277
278 printf("homo bounds: [%lf, %lf]\n", (double)homo_out, (double)homo_in);
279 printf("lumo bounds: [%lf, %lf]\n", (double)lumo_in, (double)lumo_out);
280
281
282 // JUST A CHECK...
283 if( homo_bounds.empty() )
284 {
285 printf("Interval homo_bounds is empty.\n");
286 return EXIT_FAILURE;
287 }
288 if ( lumo_bounds.empty() )
289 {
290 printf("Interval lumo_bounds is empty.\n");
291 return EXIT_FAILURE;
292 }
293
294
295 int maxit = 100;
296 real err_eig = 1e-3; // it is not needed for the new stopping criterion
297
298
299 Purification_sp2<MatrixType> Puri;
300 //Purification_sp2acc<MatrixType> Puri;
301
302 string eigenvectors_method = "square";
303 string eigenvectors_iterative_method = "lanczos";
304 real eigensolver_accuracy = TOL_EIGENSOLVER_ACC_DEFAULT;
305 int eigensolver_maxiter = 200;
306 VectorType eigVecHOMO;
307 VectorType eigVecLUMO;
308
309 IntervalType dummy;
310
311 Puri.set_eigenvectors_params(eigenvectors_method,
312 eigenvectors_iterative_method,
313 eigensolver_accuracy,
314 eigensolver_maxiter,
315 0,
316 0
317 );
318
319
320 // enable_printf_output(); // write more debug info about each iteration
321
322 Puri.initialize(F,
323 lumo_bounds,
324 homo_bounds,
325 maxit,
326 err_sub,
327 err_eig,
328 1, // 1 = new, 0 = old stopping criterion
329 normPuri,
330 normPuriStopCrit,
331 N_occ);
332
333
334
335 // RUN RECURSIVE EXPANSION
336 printf("Start recursive expansion...\n");
337 Puri.PurificationStart();
338 Puri.info.print_collected_info_printf();
339
340 // CHECK RESULT OF THE RECURSIVE EXPANSION
341 if (Puri.info.converged != 1)
342 {
343 throw std::runtime_error("SP2 did not converge!");
344 }
345 MatrixType X(Puri.X);
346 ergo_real traceX = X.trace();
347 if (template_blas_fabs(traceX - N_occ) > TOL_TRACE_ERROR_DEFAULT)
348 {
349 throw std::runtime_error("SP2: Wrong value of trace! (abs(traceX - N_occ) > TOL_TRACE_ERROR_DEFAULT)");
350 }
351
352 printf("\n");
353
354 std::vector<VectorType> eigVecOCC, eigVecUNOCC;
355 std::vector<real> eigValOCC, eigValUNOCC;
356 Puri.extract_computed_eigenpairs(eigVecUNOCC, eigVecOCC, eigValUNOCC, eigValOCC);
357 assert(eigVecOCC.size() > 0);
358 assert(eigVecUNOCC.size() > 0);
359 assert(eigVecOCC.size() == eigValOCC.size());
360 assert(eigVecUNOCC.size() == eigValUNOCC.size());
361 eigVecHOMO = eigVecOCC[0];
362 eigVecLUMO = eigVecUNOCC[0];
363
364 real homo_computed = eigValOCC[0];
365 real lumo_computed = eigValUNOCC[0];
366 printf("Computed HOMO eigenvalue is %.8lf, exact HOMO is %.8lf, homo bounds: [%lf, %lf]\n", (double)homo_computed, (double)homo, (double)homo_out, (double)homo_in);
367 printf("Computed LUMO eigenvalue is %.8lf, exact LUMO is %.8lf, lumo bounds: [%lf, %lf]\n", (double)lumo_computed, (double)lumo, (double)lumo_in, (double)lumo_out);
368
369 assert(template_blas_fabs(lumo_computed - lumo) < SCALAR_TOL);
370 assert(template_blas_fabs(homo_computed - homo) < SCALAR_TOL);
371
372 printf("\n");
373
374
375 // CHECK EIGENVECTORS
376
377 vector<int> I(N);
378 vector<int> JH(N), JL(N);
379 vector<real> eigVecHOMOExact(N), eigVecLUMOExact(N);
380 // get exact homo vector
381 for(int i = 0; i < N; ++i)
382 {
383 I[i] = i;
384 JH[i] = N_occ-1; JL[i] = N_occ;
385 }
386 Q.get_values(I, JH, eigVecHOMOExact);
387 Q.get_values(I, JL, eigVecLUMOExact);
388
389 vector<real> eigVecHOMO_vec, eigVecLUMO_vec;
390 eigVecHOMO.fullvector(eigVecHOMO_vec);
391 eigVecLUMO.fullvector(eigVecLUMO_vec);
392
393 // let maximum component will be positive
394 {
395 vector<real>::iterator ind;
396 ind = max_element(eigVecLUMOExact.begin(), eigVecLUMOExact.end(), abs_compare);
397 real a = eigVecLUMOExact[distance(eigVecLUMOExact.begin(), ind)];
398 if( a < 0 )
399 std::transform(eigVecLUMOExact.begin(), eigVecLUMOExact.end(), eigVecLUMOExact.begin(), std::negate<real>());
400
401 ind = max_element(eigVecHOMOExact.begin(), eigVecHOMOExact.end(), abs_compare);
402 a = eigVecHOMOExact[distance(eigVecHOMOExact.begin(), ind)];
403 if( a < 0 )
404 std::transform(eigVecHOMOExact.begin(), eigVecHOMOExact.end(), eigVecHOMOExact.begin(), std::negate<real>());
405
406 ind = max_element(eigVecLUMO_vec.begin(), eigVecLUMO_vec.end(), abs_compare);
407 a = eigVecLUMO_vec[distance(eigVecLUMO_vec.begin(), ind)];
408 if( a < 0 )
409 std::transform(eigVecLUMO_vec.begin(), eigVecLUMO_vec.end(), eigVecLUMO_vec.begin(), std::negate<real>());
410
411 ind = max_element(eigVecHOMO_vec.begin(), eigVecHOMO_vec.end(), abs_compare);
412 a = eigVecHOMO_vec[distance(eigVecHOMO_vec.begin(), ind)];
413 if( a < 0 )
414 std::transform(eigVecHOMO_vec.begin(), eigVecHOMO_vec.end(), eigVecHOMO_vec.begin(), std::negate<real>());
415
416
417 }
418
419 printf("Checking eigenvectors...\n");
420
421 for(int i = 0; i < N; ++i)
422 {
423 assert(template_blas_fabs(eigVecHOMO_vec[i] - eigVecHOMOExact[i]) < SCALAR_TOL);
424 assert(template_blas_fabs(eigVecLUMO_vec[i] - eigVecLUMOExact[i]) < SCALAR_TOL);
425 }
426
427 printf("\n");
428
429 // OUTPUT FIGURES
430 // ostringstream name;
431 // name << "out_norm" << ".m";
432 // Puri.gen_matlab_file_norm_diff(name.str().c_str());
433 // name.str("");
434 //
435 // name << "out_tau" << ".m";
436 // Puri.gen_matlab_file_threshold(name.str().c_str());
437 // name.str("");
438 //
439 // name << "out_nnz" << ".m";
440 // Puri.gen_matlab_file_nnz(name.str().c_str());
441 // name.str("");
442 //
443 // name << "out_eigs" << ".m";
444 // Puri.gen_matlab_file_eigs(name.str().c_str());
445 // name.str("");
446 //
447 // name << "out_cond" << ".m";
448 // Puri.gen_matlab_file_cond_num(name.str().c_str());
449 // name.str("");
450 //
451 // name << "out_time" << ".m";
452 // Puri.gen_matlab_file_time(name.str().c_str());
453
454
455 F.clear();
456 Puri.clear();
457
458 printf("DONE!\n");
459
460
461 return EXIT_SUCCESS;
462 }
463
464 #endif
465