LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-25 05:43:40 +01:00
Code Issues Releases Wiki Activity

Fix #66: when mc_sign_type is complex

Browse Source 
- for complex type sign, we do not do the previous check
(to be improved), since isinf, signbit are not defined for
complex numbers.
This commit is contained in:
Olivier Parcollet 2014-03-28 15:37:37 +01:00
parent c996c3ff7d
commit 9c979ceb0b
1 changed files with 16 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
triqs/mc_tools/mc_move_set.hpp
View File

@ -162,6 +162,20 @@ namespace triqs { namespace mc_tools {
normaliseProba();// ready to run after each add !
}
private:
bool attempt_treat_infinite_ratio(std::complex<double>, double &) { return true; }
bool attempt_treat_infinite_ratio(double rate_ratio, double &abs_rate_ratio) {
bool is_inf = std::isinf(rate_ratio);
if (is_inf) { // in case the ratio is infinite
abs_rate_ratio = 100; // 1.e30; // >1 for metropolis
try_sign_ratio = (std::signbit(rate_ratio) ? -1 : 1); // signbit -> true iif the number is negative
}
return !is_inf;
}
public:
/**
* - Picks up one of the move at random (weighted by their proposition probability),
* - Call attempt method of that move
@ -184,15 +198,12 @@ namespace triqs { namespace mc_tools {
#endif
MCSignType rate_ratio = current->attempt();
double abs_rate_ratio;
if (std::isinf(rate_ratio)) { // in case the ratio is infinite
abs_rate_ratio = 100; //1.e30; // >1 for metropolis
try_sign_ratio = (std::signbit(rate_ratio) ? -1 : 1); // signbit -> true iif the number is negative
} else {
if (attempt_treat_infinite_ratio(rate_ratio, abs_rate_ratio)) { // in case the ratio is infinite
if (!std::isfinite(std::abs(rate_ratio)))
TRIQS_RUNTIME_ERROR << "Monte Carlo Error : the rate is not finite in move " << name_of_currently_selected();
abs_rate_ratio = std::abs(rate_ratio);
#ifdef TRIQS_TOOLS_MC_DEBUG
std::cerr << " Metropolis ratio " << rate_ratio << ". Abs(Metropolis ratio) " << abs_rate_ratio << std::endl;
std::cerr << " Metropolis ratio " << rate_ratio << ". Abs(Metropolis ratio) " << abs_rate_ratio << std::endl;
#endif
assert((abs_rate_ratio >= 0));
try_sign_ratio = (abs_rate_ratio > 1.e-14 ? rate_ratio / abs_rate_ratio : 1); // keep the sign