1 //////////////////////////////////////////////////////////////////////////////
2 /// \file formic/utils/lmyengine/block_detai.cpp
3 ///
4 /// \brief implementation file for block-recursive algorithm functions
5 ///
6 //////////////////////////////////////////////////////////////////////////////
7
8 #include <vector>
9 #include <list>
10 #include <string>
11 #include <numeric>
12 #include <cassert>
13 #include <algorithm>
14 #include <cmath>
15 #include <sstream>
16 //#include <mpi.h>
17
18 #include <boost/format.hpp>
19 #include <boost/shared_ptr.hpp>
20
21 #include "formic/utils/matrix.h"
22 #include "formic/utils/numeric.h"
23 #include "formic/utils/mpi_interface.h"
24 #include "block_detail.h"
25 #include "formic/utils/lmyengine/var_dependencies.h"
26 #include "formic/utils/lmyengine/eigen_solver.h"
27 #include "formic/utils/lmyengine/davidson_solver.h"
28 #include "formic/utils/lmyengine/spam_solver.h"
29
30 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
31 /// \brief Get the begining and end indices and the length for each block of variables
32 ///
33 /// \param[in] nvar the number of variables that will be divided into blocks
34 /// \param[in] nblock the number of blocks
35 /// \param[out] block_beg on exit, the length nblock vector of indices marking the begining of each block
36 /// \param[out] block_end on exit, the length nblock vector of block lengths
37 ///
38 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
brlm_get_block_info(const int nvar,const int nblock,std::vector<int> & block_beg,std::vector<int> & block_end,std::vector<int> & block_len)39 void cqmc::engine::brlm_get_block_info(const int nvar, const int nblock, std::vector<int> & block_beg, std::vector<int> & block_end, std::vector<int> & block_len) {
40
41 // initialize output vectors
42 block_beg.assign(nblock, 0);
43 block_end.assign(nblock, 0);
44 block_len.assign(nblock, 0);
45
46 // get the length of most blocks
47 const int base_block_len = nvar / nblock;
48
49 // get the residual
50 const int base_block_rem = nvar % nblock;
51
52 int beg = 0;
53
54 // loop over blocks
55 for (int i = 0; i < nblock; i++) {
56 block_beg.at(i) = beg;
57 block_len.at(i) = base_block_len + ( i < base_block_rem ? 1 : 0 );
58 beg += block_len.at(i);
59 block_end.at(i) = beg;
60 }
61 if ( beg != nvar )
62 throw formic::Exception("after block index steup, beg = %i but should have been %i in cqmc::engin::brlm_get_block_info") % beg % nvar;
63 }
64
65 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
66 /// \brief Computes the first so many important update directions in the basis of those that define
67 /// the basis for the provided Hamiltinonian and overlap matrices
68 ///
69 /// \param[in] nkeep the number of important update directions to keep
70 /// \param[in] curr_e the current energy
71 /// \param[in] shift_i magnitude of the identity shift
72 /// \param[in] shift_s magnitude of the overlap shift
73 /// \param[in] ovl_thresh threshold below which overlap eigenvalues will be truncated
74 /// \param[in] hh the Hamiltonian matrix in the basis of possible directions.
75 /// the first of which is the current wavefunction
76 /// \param[in] ss the overlap matrix in the basis of possible directions,
77 /// the first of which is the current wavefunction
78 /// \param[in] dd the identity shift matrix in the basis of possible directions,
79 /// the first of which is the current wavefunction
80 /// \param[in] ostream the output stream
81 ///
82 /// \return a matrix holding the best nkeep directions in its column, the first of which
83 /// is always the current wavefunction
84 ///
85 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
get_important_brlm_dirs(const int nkeep,const double curr_e,const double shift_i,const double shift_s,const double ovl_thresh,const formic::Matrix<double> & hh,const formic::Matrix<double> & ss,const formic::Matrix<double> & dd,std::ostream & output)86 formic::Matrix<double> cqmc::engine::get_important_brlm_dirs(const int nkeep,
87 const double curr_e,
88 const double shift_i,
89 const double shift_s,
90 const double ovl_thresh,
91 const formic::Matrix<double> & hh,
92 const formic::Matrix<double> & ss,
93 const formic::Matrix<double> & dd,
94 std::ostream & output) {
95
96 ////////////////////////////////////////////////////////////////////////
97 // Note in this function we deal with four sets of basis vectors
98 //
99 // Basis 1: The initial basis, in which hh and ss are defined.
100 // The current wave function is the first basis vector.
101 //
102 // Basis 2: Basis after projecting out the current wavefunction from
103 // the second through through last basis vectors of basis 1.
104 //
105 // Basis 3: The basis consisting of the second through last vectors
106 // of basis 2, i.e. everything except the current wfn.
107 //
108 // Basis 4: The basis after projecting out all previous updates from
109 // the vectors in basis 2 and then orthonormalizing, dropping
110 // vectors with pre-normalized norms below ovl_thresh.
111 // The first vector in this basis is still the current wfn.
112 //
113 /////////////////////////////////////////////////////////////////////////
114
115 // get the number of vectors in basis 1
116 const int nd = hh.rows();
117
118
119 // sanity checks
120 if ( hh.rows() != hh.cols() )
121 throw formic::Exception("hh (%i by %i) should be square in cqmc::engine::get_important_brlm_dirs") % hh.rows() % hh.cols();
122 if ( hh.rows() != ss.rows() || hh.cols() != ss.cols() )
123 throw formic::Exception("ss dimension (%i by %i) should match hh dimension (%i by %i) in cqmc::engine::get_important_brlm_dirs")
124 % ss.rows() % ss.cols() % hh.rows() % hh.cols();
125 //if ( std::abs( 1.0 - std::abs(ss.at(0,0)) ) > 1.0e-14 )
126 // throw formic::Exception("expected ss(0,0) to be 1.0 in cqmc::engine::get_important_brlm_dirs");
127
128 // put the vectors defining basis 2 in the columns of a matrix
129 formic::Matrix<double> basis2 = formic::identity_matrix<double>(nd);
130 basis2.at(0,0) = 1.0/std::abs(std::sqrt(ss(0,0)));
131 for (int i = 1; i < nd; i++)
132 basis2.at(0,i) -= (ss(i,0)/ss(0,0));
133
134 // get H and S in basis 2, including the identity shift effect on H
135 formic::Matrix<double> hh2 = basis2.t() * ( hh + shift_i * dd ) * basis2;
136 formic::Matrix<double> ss2 = basis2.t() * ss * basis2;
137 //output << ss2.print("%12.6f", "ss2");
138
139 //output << hh2.print("%12.4f", "hh2") << std::endl;
140
141 // extract overlap matrix for basis 3
142 formic::Matrix<double> ss3(nd-1,nd-1);
143 for (int i = 0; i < nd-1; i++) {
144 for (int j = 0; j < nd-1; j++) {
145 ss3.at(i,j) = ss2.at(i+1,j+1);
146 }
147 }
148
149 // initialize matrix of update direction coefficients in basis 2
150 formic::Matrix<double> update_dirs(nd, nkeep, 1.0);
151
152 // ietratively find the best nkeep directions
153 for (int nkept = 0; nkept < nkeep; nkept++) {
154
155 // make sure previous update directions are excluded by projecting them out
156 formic::Matrix<double> proj_dirs = formic::identity_matrix<double>(nd-1);
157 if ( nkept > 0 ) {
158
159 // get the matrix of update coefficients that excludes the coeff on the currenr wfn
160 formic::Matrix<double> bigV(nd-1, nkept);
161 for (int j = 0; j < nkept; j++) {
162 for (int i = 0; i < nd-1; i++) {
163 bigV.at(i,j) = update_dirs.at(i+1, j);
164 }
165 formic::ColVec<double> current_col(nd-1, bigV.col_begin(j));
166 current_col /= std::abs(std::sqrt(current_col.norm2()));
167 }
168
169 // get "bare" overlap between these vectors
170 formic::Matrix<double> bares = bigV.t() * bigV;
171
172 // check the smallest eigenvalue of the bare overlap
173 formic::ColVec<double> bare_evals;
174 formic::Matrix<double> bare_evecs;
175 bares.sym_eig(bare_evals, bare_evecs);
176
177 // project these directions out of the full set of directions
178 proj_dirs -= bigV * bares.inv() * bigV.t();
179 }
180
181 // get overlap in these projected directions
182 formic::Matrix<double> ps = proj_dirs.t() * ss3 * proj_dirs;
183 //output << ps.print("%12.6f", "ps") << std::endl;
184
185 // diagonalize the projected overlap
186 formic::ColVec<double> ps_evals;
187 formic::Matrix<double> ps_evecs;
188 ps.sym_eig(ps_evals, ps_evecs);
189 //output << ps_evecs.print("%12.6f", "ps_evecs") << std::endl;
190 //output << ps_evals.print("%12.6f", "ps_evals") << std::endl;
191
192 // count number of projected overlap eigenvalues above the threshold
193 //const int nabove = std::accumulate(ps_evals.begin(), ps_evals.end(), int(0), [ovl_thresh] (int i, double x) { return i + ( x > ovl_thresh ? 1 : 0 ); });
194 const int nabove = std::accumulate(ps_evals.begin(), ps_evals.end(), int(0), [ovl_thresh] (int i, double x) { return i + ( x > ovl_thresh ? 1 : 0 ); });
195
196 const int nbelow = ps_evals.size() - nabove;
197
198 // get coefficients for basis 4 vectors in terms of basis 2 vectors
199 formic::Matrix<double> basis4(nd, nabove+1, 0.0);
200
201 // initial wfn
202 basis4.at(0,0) = 1.0;
203
204 for (int i = 0; i < nabove; i++) {
205 const double eval_sqrt = std::abs(std::sqrt(ps_evals.at(nbelow+i)));
206 for (int j = 0; j < nd-1; j++) {
207 basis4.at(j+1, i+1) = ps_evecs.at(j, i+nbelow) / eval_sqrt;
208 }
209 }
210 //output << basis4.print("%12.6f", "basis4") << std::endl;
211
212 // compute overlap in basis 4 and check that it is the identity matrix
213 formic::Matrix<double> ss4 = basis4.t() * ss2 * basis4;
214 //output << ss4.print("%12.6f", "ss4");
215 for (int i = 0; i < ss4.rows(); i++) {
216 for (int j = 0; j < ss4.cols(); j++) {
217 if ( std::abs( ss4.at(i,j) - ( i == j ? 1.0 : 0.0) ) > 1.0e-6 )
218 throw formic::Exception("failure to move to orthonormal basis in cqmc::engine::get_important_brlm_dirs, ss4(%i,%i) = %.8e")
219 % i % j % ss4.at(i,j);
220 }
221 }
222
223 // compute Hamiltonian in basis 4
224 formic::Matrix<double> hh4 = basis4.t() * hh2 * basis4;
225 //output << hh4.print("%12.6f", "hh4") << std::endl;
226
227 // apply overlap shift to the Hamiltonian
228 for (int i = 1; i < hh4.rows(); i++)
229 hh4.at(i,i) += shift_s;
230
231 // solve for the lowest energy in eigenvector in basis 4
232 formic::ColVec<std::complex<double> > e_evals;
233 formic::Matrix<std::complex<double> > e_evecs;
234 hh4.nonsym_eig(e_evals, e_evecs);
235
236 // find the lowest energy eigenvalue
237 double lowest_eval = e_evals.at(0).real();
238 int lowest_index = 0;
239 for (int i = 1; i < e_evals.size(); i++) {
240 if ( e_evals.at(i).real() < lowest_eval ) {
241 lowest_eval = e_evals.at(i).real();
242 lowest_index = i;
243 }
244 }
245
246 output << boost::format("shift_s %.2e vec %i target = %20.12f + %16.2f i delta = %17.12f")
247 % shift_s % nkept % e_evals.at(lowest_index).real() % e_evals.at(lowest_index).imag() % (e_evals.at(lowest_index).real() - curr_e)
248 << std::endl;
249
250 // check that the eigenvalue is real
251 if ( std::abs( e_evals.at(lowest_index).imag() ) > 1.0e-6 )
252 throw formic::Exception("lowest_eval has an imaginary component of %.2e in cqmc::engine::get_lm_brnr_update") % e_evals.at(lowest_index).imag();
253
254 // check that the eigenvector has weight on the initial wavefunction
255 if ( std::abs( e_evecs.at(0, lowest_index) ) < 1.0e-6 )
256 throw formic::Exception("lowest_evec has a small weight of %.2e on the initial wave function in cqmc::engine::get_lm_brnr_update") % std::abs(e_evecs.at(0,lowest_index));
257
258 // attempt to get the eigenvector in terms of real numbers
259 formic::ColVec<double> evec_real(e_evecs.rows(), 0.0);
260 for (int i = 0; i < e_evecs.rows(); i++) {
261 std::complex<double> val = e_evecs.at(i, lowest_index) / e_evecs.at(0, lowest_index);
262 if ( std::abs(val.imag()) > 1.0e-6 )
263 throw formic::Exception("lowest_evec element %i has an un-removed imaginary component of %.2e in cqmc::engine::get_lm_brnr_update") % i % val.imag();
264 evec_real.at(i) = val.real();
265 }
266
267 // convert the eigenvector to basis 2
268 ( basis4 * evec_real ) >>= update_dirs.col_begin(nkept);
269
270 // check that this eigenvector still has weight on the initial wfn
271 if ( std::abs( update_dirs.at(0, nkept) ) < 1.0e-6 )
272 throw formic::Exception("basis 2 eigenvector has a small weight of %.2e on the initial wave function in cqmc::engine::get_lm_brnr_update")
273 % std::abs( update_dirs.at(0, nkept) );
274
275 // scale the eigenvector in basis 2 so that the initial wavefunction has unit weight
276 formic::ColVec<double>(update_dirs.rows(), update_dirs.col_begin(nkept)) /= *update_dirs.col_begin(nkept);
277
278 }
279
280 // convert the eigenvectors to basis 1
281 update_dirs = basis2 * update_dirs;
282
283 // scale the eigenvectors so that in basis 1 they have unit weight on the initial wave function
284 for (int nkept = 0; nkept < nkeep; nkept++) {
285 if ( std::abs( *update_dirs.col_begin(nkept) ) < 1.0e-6 )
286 throw formic::Exception("basis 1 eigenvector has an small weight of %.2e on the initial wave function in cqmc::engine::get_lm_brnr_update")
287 % std::abs( *update_dirs.col_begin(nkept) );
288 formic::ColVec<double>(update_dirs.rows(), update_dirs.col_begin(nkept)) /= (*update_dirs.col_begin(nkept));
289 }
290
291 // clean up memory used by matrices and vectors
292 formic::reusable_array_garbage_collect();
293
294 // return the eigenvectors in basis 1
295 return update_dirs;
296
297 }
298
299 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
300 /// \brief Computes the first so many important update directions in the basis of those that define
301 /// the basis for the provided Hamiltinonian and overlap matrices
302 ///
303 /// \param[in] nkeep the number of important update directions to keep
304 /// \param[in] curr_e the current energy
305 /// \param[in] omega energy shift used for harmonic davidson
306 /// \param[in] shift_i magnitude of the identity shift
307 /// \param[in] shift_s magnitude of the overlap shift
308 /// \param[in] ovl_thresh threshold below which overlap eigenvalues will be truncated
309 /// \param[in] hh the Hamiltonian matrix in the basis of possible directions.
310 /// the first of which is the current wavefunction
311 /// \param[in] ss the overlap matrix in the basis of possible directions,
312 /// the first of which is the current wavefunction
313 /// \param[in] dd the identity shift matrix in the basis of possible directions,
314 /// the first of which is the current wavefunction
315 /// \param[in] ostream the output stream
316 ///
317 /// \return a matrix holding the best nkeep directions in its column, the first of which
318 /// is always the current wavefunction
319 ///
320 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
get_important_brlm_dirs_davidson(const formic::VarDeps * dep_ptr,const int nkeep,const double omega,const double curr_cost,const double shift_i,const double shift_s,const double ovl_thresh,formic::Matrix<double> & hh,formic::Matrix<double> & ss,formic::Matrix<double> & dd,std::ostream & output)321 formic::Matrix<double> cqmc::engine::get_important_brlm_dirs_davidson(const formic::VarDeps * dep_ptr,
322 const int nkeep,
323 const double omega,
324 const double curr_cost,
325 const double shift_i,
326 const double shift_s,
327 const double ovl_thresh,
328 formic::Matrix<double> & hh,
329 formic::Matrix<double> & ss,
330 formic::Matrix<double> & dd,
331 std::ostream & output) {
332
333
334 // get rank number and number of ranks
335 int my_rank = formic::mpi::rank();
336
337 // check to see whether nkeep is one and throw out an error if not
338 if ( nkeep != 1 )
339 throw formic::Exception("nkeep must be 1 in get_important_brlm_dirs_davidson!");
340
341 // get the number of vectors in basis 1
342 int nd = hh.rows();
343 formic::mpi::bcast(&nd, 1);
344 //MPI_Bcast(&nd, 1, MPI_INT, 0, MPI_COMM_WORLD);
345
346 // get the flag that whether the calculation is ground or targeted excited
347 bool ground = true;
348 if ( my_rank == 0 )
349 ground = ( std::abs(ss.at(0,0) - 1.0) < 1.0e-6 );
350 formic::mpi::bcast(&ground, 1);
351 //MPI_Bcast(&ground, 1, MPI::BOOL, 0, MPI_COMM_WORLD);
352
353 // initial energy/target fn value
354 double init_cost = 0.0;
355 if ( my_rank == 0 )
356 init_cost = hh.at(0,0)/ss.at(0,0);
357 formic::mpi::bcast(&init_cost, 1);
358 //MPI_Bcast(&init_cost, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
359
360 // sanity checks
361 if ( my_rank == 0 ) {
362 if ( hh.rows() != hh.cols() )
363 throw formic::Exception("hh (%i by %i) should be square in cqmc::engine::get_important_brlm_dirs") % hh.rows() % hh.cols();
364 if ( hh.rows() != ss.rows() || hh.cols() != ss.cols() )
365 throw formic::Exception("ss dimension (%i by %i) should match hh dimension (%i by %i) in cqmc::engine::get_important_brlm_dirs")
366 % ss.rows() % ss.cols() % hh.rows() % hh.cols();
367 }
368
369 // initialize matrix of update direction coefficients
370 formic::Matrix<double> update_dirs(nd, nkeep, 1.0);
371
372 // add identity shift to the Hamiltonian
373 formic::Matrix<double> m_hh = hh + shift_i * dd;
374
375 // create eigensolver
376 std::vector<double> vf_var;
377 boost::shared_ptr< cqmc::engine::EigenSolver<double> > eigensolver(new cqmc::engine::DavidsonLMHD<double>(dep_ptr,
378 nd,
379 60,
380 1.0e-7,
381 ovl_thresh,
382 false, // var deps use?
383 true, // chase lowest
384 false, // chase closest
385 ground, // ground state?
386 false, // variance correct
387 true, // build matrix
388 vf_var,
389 init_cost, // initial energy
390 0.0, // initial variance
391 omega, // omega
392 0.0, // var weight
393 20.0, // maximun energy change
394 0.0, // total weight
395 0.0, // total weight
396 hh,
397 hh,
398 m_hh,
399 ss,
400 ss));
401
402
403 // iteratively find the best nkeep directions CURRENTLY ONLY WORK FOR NKEEP=1!!!!!
404 for (int nkept = 0; nkept < nkeep; nkept++) {
405
406 // reset the eigensolver
407 eigensolver -> reset();
408
409 // update shift
410 eigensolver -> update_lm_shift(0.0, shift_s);
411
412 // add first krylov vector
413 {
414 formic::ColVec<double> temp(hh.cols());
415 for (int j = 0; j < temp.size(); j++)
416 temp.at(j) = ( j == 0 ? 1.0 : 0.0);
417 eigensolver -> add_krylov_vector(temp);
418 }
419
420
421 // solve the eigenvalue problem
422 double davidson_eval;
423 bool solve_shift = eigensolver -> solve(davidson_eval, output);
424
425 // scale the eigenvector so that the initial wavefunction has the unit weight
426 if ( my_rank == 0 ) {
427 eigensolver -> convert_to_wf_coeff();
428 formic::ColVec<double> evec_eigen = eigensolver -> wf_coeff();
429
430 // compute the largest weight on the derivative vector
431 double max_update_abs_value = std::abs(evec_eigen.at(1));
432
433 // store this eigenvector
434 evec_eigen >>= update_dirs.col_begin(nkept);
435
436 for (int k = 0; k < nd; k++) {
437 if ( k != 0 )
438 max_update_abs_value = std::abs(evec_eigen.at(k)) > max_update_abs_value ? std::abs(evec_eigen.at(k)) : max_update_abs_value;
439 }
440 output << boost::format("The largest weight on the derivative vector for shift %.4e is %.6e") % shift_i % max_update_abs_value << std::endl << std::endl;
441 }
442 }
443
444 // clean up memory used by matrices and vectors
445 formic::reusable_array_garbage_collect();
446
447 // return the eigenvectors in basis 1
448 return update_dirs;
449 }
450
451