diff --git a/CMakeLists.txt b/CMakeLists.txt index 922a0cf..6618e89 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,8 @@ if(USE_QUAD_PRECISION) endif() #ALPSCore disable debug for gf library -set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DBOOST_DISABLE_ASSERTS -DNDEBUG") +#(please do not set NDEBUG for DEBUG build. That would would disable runtime checks in the solver) +set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -DBOOST_DISABLE_ASSERTS -DNDEBUG") #include directories set(CHYB_LIBRARY_INCLUDE_DIRS ${EIGEN3_INCLUDE_DIR} ${MPI_CXX_INCLUDE_PATH}) #rest taken care of by libraries dependencies @@ -26,9 +27,21 @@ include_directories(${CHYB_LIBRARY_INCLUDE_DIRS}) #source files set(LIB_FILES ./src/sliding_window.cpp ./src/legendre.cpp ./src/model.cpp ./src/model_detail/clustering.cpp src/operator_util.cpp ./src/util.cpp) - ADD_LIBRARY(hyb ${LIB_FILES}) -target_link_libraries(hyb ${ALPSCore_LIBRARIES} ${MPI_CXX_LIBRARIES} ${Boost_LIBRARIES}) + +#Compiler dependent libraries +set(EXTRA_LIBS "") +if (USE_QUAD_PRECISION) + if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + list(APPEND EXTRA_LIBS "quadmath") + endif() + if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") + set(CMAKE_CXX_FLAGS " -Qoption,cpp,--extended_float_type ${CMAKE_CXX_FLAGS}") + endif() +endif() + +#Set link libraries +target_link_libraries(hyb ${ALPSCore_LIBRARIES} ${MPI_CXX_LIBRARIES} ${Boost_LIBRARIES} ${EXTRA_LIBS}) #executable ADD_EXECUTABLE(hybmat ./src/main.cpp) diff --git a/src/impurity.hpp b/src/impurity.hpp index 2c67276..9f886a0 100755 --- a/src/impurity.hpp +++ b/src/impurity.hpp @@ -124,7 +124,7 @@ class HybridizationSimulation : public alps::mcbase long sweeps; // sweeps done alps::ResizableMatrix M; SCALAR sign; // the sign of w=Z_k_up*Z_k'_down*trace - EXTENDED_SCALAR det; + //EXTENDED_SCALAR det; EXTENDED_SCALAR trace; // matrix trace typedef typename std::iterator_traits::iterator>::value_type mytpe; diff --git a/src/impurity.ipp b/src/impurity.ipp index a77f7ac..e1f545c 100755 --- a/src/impurity.ipp +++ b/src/impurity.ipp @@ -58,7 +58,7 @@ HybridizationSimulation::HybridizationSimulation(parameters_type cons sweeps(0), //sweeps done up to now M(0,0), sign(1), - det(1), + //det(1), trace(std::numeric_limits::max()), N_shift_flavor(FLAVORS, static_cast(p["N_SHIFT"])), sliding_window(p_model.get(), BETA), @@ -163,7 +163,7 @@ void HybridizationSimulation::update() { for (int flavor = 0; flavor < FLAVORS; flavor++) { const int flavor_target = (int) (random() * FLAVORS); boost::tuple r = - insert_remove_pair_flavor(random, flavor_target,flavor_target,det, BETA, + insert_remove_pair_flavor(random, flavor_target,flavor_target, BETA, order_creation_flavor,order_annihilation_flavor, creation_operators, annihilation_operators, M, sign, trace, operators, max_distance_pair, sliding_window, max_order); @@ -188,7 +188,7 @@ void HybridizationSimulation::update() { int a_flavor = (int) (random() * FLAVORS); insert_remove_pair_flavor(random, c_flavor, a_flavor, - det, BETA, + BETA, order_creation_flavor, order_annihilation_flavor, creation_operators, @@ -203,7 +203,7 @@ void HybridizationSimulation::update() { /**** shift an operator ****/ for (int ns = 0; ns < N_shift*FLAVORS; ns++) { sanity_check(); - boost::tuple r = shift_lazy(random, det, BETA, + boost::tuple r = shift_lazy(random, BETA, creation_operators, annihilation_operators, M, sign, trace, operators, max_distance_shift, @@ -220,9 +220,9 @@ void HybridizationSimulation::update() { } } - if (my_isnan(det)) { - throw std::runtime_error("det is NaN after insertion/removal/shift update"); - } + //if (my_isnan(det)) { + //throw std::runtime_error("det is NaN after insertion/removal/shift update"); + //} if (my_isnan(trace)) { throw std::runtime_error("trace is NaN after insertion/removal/shift update"); } @@ -237,9 +237,9 @@ void HybridizationSimulation::update() { //Perform global updates which might cost O(beta) expensive_updates(); - if (my_isnan(det)) { - throw std::runtime_error("det is NaN after global update"); - } + //if (my_isnan(det)) { + //throw std::runtime_error("det is NaN after global update"); + //} if (my_isnan(trace)) { throw std::runtime_error("trace is NaN after global update"); } @@ -295,6 +295,9 @@ void HybridizationSimulation::measure() { measurements["order"] << order_creation_meas; } + //measurements["AbsDeterminant"] << convert_to_double(myabs(det)); + measurements["AbsTrace"] << convert_to_double(myabs(trace)); + // measure acceptance rate measurements["Insertion_attempted"] << to_std_vector(weight_vs_distance.get_counter()); measurements["Shift_attempted"] << to_std_vector(weight_vs_distance_shift.get_counter()); @@ -438,7 +441,7 @@ void HybridizationSimulation::expensive_updates() { for (int itry = 0; itry < swap_vector.size(); ++itry) { const int iupdate = execute_ordering[itry]; - const bool accepted = exchange_flavors(random, det, BETA, creation_operators, + const bool accepted = exchange_flavors(random, BETA, creation_operators, annihilation_operators, order_creation_flavor, order_annihilation_flavor, M, sign, trace, operators, @@ -449,9 +452,9 @@ void HybridizationSimulation::expensive_updates() { const int source_template = swap_vector[iupdate].second; if (accepted) { swap_acc_rate[iupdate].accepted(); - if (my_isnan(det)) { - throw std::runtime_error("det is NaN after swap update"+boost::lexical_cast(itry)); - } + //if (my_isnan(det)) { + //throw std::runtime_error("det is NaN after swap update"+boost::lexical_cast(itry)); + //} if (my_isnan(trace)) { throw std::runtime_error("trace is NaN after swap update"+boost::lexical_cast(itry)); } @@ -466,15 +469,15 @@ void HybridizationSimulation::expensive_updates() { //Shift operators to restore translational symmetry if (do_global_shift) { - const bool accepted = global_shift(random, det, BETA, creation_operators, + const bool accepted = global_shift(random, BETA, creation_operators, annihilation_operators, order_creation_flavor, order_annihilation_flavor, M, sign, trace, operators, sliding_window); if (accepted) { global_shift_acc_rate.accepted(); - if (my_isnan(det)) { - throw std::runtime_error("det is NaN after global shift"); - } + //if (my_isnan(det)) { + //throw std::runtime_error("det is NaN after global shift"); + //} if (my_isnan(trace)) { throw std::runtime_error("trace is NaN after global shift"); } @@ -631,7 +634,13 @@ void HybridizationSimulation::sanity_check() const { } // compute determinants - EXTENDED_SCALAR det_new = cal_det(creation_operators, annihilation_operators, M_new, BETA, p_model->get_F()); + SCALAR sign_det = 1.0; + if (creation_operators.size() > 0) { + const std::vector& det_vec = cal_det_as_vector(creation_operators, annihilation_operators, M_new, BETA, p_model->get_F()); + for (int i=0; i::sanity_check() const { throw std::runtime_error("trace != trace_new"); } - if (std::abs((EXTENDED_SCALAR(det/det_new).template convert_to()-1.0))>1E-5) { - throw std::runtime_error("det_new != det"); - } + //if (std::abs((EXTENDED_SCALAR(det/det_new).template convert_to()-1.0))>1E-5) { + //throw std::runtime_error("det_new != det"); + //} SCALAR sign_overall_new = dsign(sign_new) *dsign(trace).template convert_to() - *dsign(det_new).template convert_to(); + *sign_det; + //std::cout << "debug det "<< det << " " << det_new << std::endl; //std::cout << "debug sign "<< sign_overall_new << " " << sign << std::endl; if (std::abs(sign_overall_new/sign-1.0)>1E-5) { throw std::runtime_error("sign_overall_new != sign"); diff --git a/src/impurity_init.ipp b/src/impurity_init.ipp index 1c6e579..7cca6b8 100755 --- a/src/impurity_init.ipp +++ b/src/impurity_init.ipp @@ -11,6 +11,8 @@ void HybridizationSimulation::create_observables() { measurements << alps::accumulators::LogBinningAccumulator >("n"); measurements << alps::accumulators::LogBinningAccumulator("Sign"); + //measurements << alps::accumulators::NoBinningAccumulator("AbsDeterminant"); + measurements << alps::accumulators::NoBinningAccumulator("AbsTrace"); measurements << alps::accumulators::NoBinningAccumulator >("order"); measurements << alps::accumulators::NoBinningAccumulator >("Insertion_attempted"); diff --git a/src/model.cpp b/src/model.cpp index 6f200d7..0772a55 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1,6 +1,7 @@ +#include + #include "operator.hpp" #include "model.hpp" -#include #include "./model_detail/model.ipp" #include "./model_detail/eigenbasis.ipp" diff --git a/src/model_detail/eigenbasis.ipp b/src/model_detail/eigenbasis.ipp index 353aac0..2e1f5af 100644 --- a/src/model_detail/eigenbasis.ipp +++ b/src/model_detail/eigenbasis.ipp @@ -187,7 +187,8 @@ void ImpurityModelEigenBasis::build_basis(const alps::params& par) { } else { overflow_prevention = 0; } - Base::reference_energy_ = (0.5*(eigenvalue_max+eigenvalue_min)-overflow_prevention); + //Base::reference_energy_ = (0.5*(eigenvalue_max+eigenvalue_min)-overflow_prevention); + Base::reference_energy_ = eigenvalue_min; if (Base::verbose_) { std::cout << "Reference energy " << Base::reference_energy_ << std::endl; } diff --git a/src/moves.hpp b/src/moves.hpp index ab9b99d..80e263d 100755 --- a/src/moves.hpp +++ b/src/moves.hpp @@ -72,7 +72,7 @@ inline std::pair operators_remove_nocopy(operator_container_t &operato // update_type, accepted, distance, acceptance probability, valid_move_generated // Note: this function is too long. Better to split into several functions. template -boost::tuple insert_remove_pair_flavor(R& rng, int creation_flavor, int annihilation_flavor, EXTENDED_SCALAR & det, double BETA, +boost::tuple insert_remove_pair_flavor(R& rng, int creation_flavor, int annihilation_flavor, double BETA, std::vector & order_creation_flavor, std::vector & order_annihilation_flavor, operator_container_t& creation_operators, operator_container_t& annihilation_operators, M_TYPE& M, SCALAR & sign, EXTENDED_SCALAR & trace, operator_container_t& operators, double cutoff, @@ -185,7 +185,6 @@ boost::tuple insert_remove_pair_flavor(R& rng, int compute_M_row_column_up(row, column, M, sign_Fs, Fe_M, det_rat); assert(!my_isnan(det_rat)); - det *= det_rat; sign *= prob/std::abs(prob); trace=trace_new; return boost::make_tuple(0,true,std::abs(t_ins-t_rem),prob,valid_move_generated); @@ -256,7 +255,6 @@ boost::tuple insert_remove_pair_flavor(R& rng, int } if (r_th < std::abs(prob)) { // move accepted - det *= det_rat; sign *= prob/std::abs(prob); @@ -282,7 +280,7 @@ boost::tuple insert_remove_pair_flavor(R& rng, int template boost::tuple -shift_lazy(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t & creation_operators, +shift_lazy(R & rng, double BETA, operator_container_t & creation_operators, operator_container_t & annihilation_operators, M_TYPE & M, SCALAR &sign, EXTENDED_SCALAR &trace, operator_container_t & operators, double distance, SLIDING_WINDOW& sliding_window) @@ -433,20 +431,11 @@ shift_lazy(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t & c sign *= prob/std::abs(prob); trace = trace_new; //Note that det_rat is computed without taking into account exchanges of rows and columns in the matrix. This yields a sign flip. - if (num_row_or_column_swaps % 2 == 1) { - det *= -det_rat; - } else { - det *= det_rat; - } -#ifndef NDEBUG - if (!equal_det(EXTENDED_SCALAR(det), static_cast(EXTENDED_SCALAR(1.0)/M.template safe_determinant())) ) { - std::cout << "ERROR IN SHIFT UPDATE " - << "det (fast update) = " << det - << "det (M^-1) = " << 1./M.determinant() - << std::endl; - exit(1); - } -#endif + //if (num_row_or_column_swaps % 2 == 1) { + //det *= -det_rat; + //} else { + //det *= det_rat; + //} } else { safe_erase(operators,new_operator); safe_insert(operators,removed_op); @@ -463,7 +452,7 @@ shift_lazy(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t & c */ template bool -global_shift(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t & creation_operators, operator_container_t & annihilation_operators, +global_shift(R & rng, double BETA, operator_container_t & creation_operators, operator_container_t & annihilation_operators, std::vector & order_creation_flavor, std::vector & order_annihilation_flavor, M_TYPE & M, SCALAR &sign, EXTENDED_SCALAR &trace, operator_container_t & operators, SLIDING_WINDOW& sliding_window @@ -502,7 +491,6 @@ global_shift(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t if (rng() < std::abs(prob)) { sign *= prob/std::abs(prob); trace = trace_new; - det *= det_rat; std::swap(operators, operators_new); std::swap(creation_operators, creation_operators_new); std::swap(annihilation_operators, annihilation_operators_new); @@ -518,7 +506,7 @@ global_shift(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t template bool -exchange_flavors(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container_t & creation_operators, operator_container_t & annihilation_operators, +exchange_flavors(R & rng, double BETA, operator_container_t & creation_operators, operator_container_t & annihilation_operators, std::vector & order_creation_flavor, std::vector & order_annihilation_flavor, M_TYPE & M, SCALAR &sign, EXTENDED_SCALAR &trace, operator_container_t & operators, @@ -530,6 +518,10 @@ exchange_flavors(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container if (creation_operators.size() == 0) { return false; } + const int pert_order = creation_operators.size(); + + const std::vector& det_vec = cal_det_as_vector(creation_operators, annihilation_operators, M, BETA, + sliding_window.get_p_model()->get_F()); //compute new trace operator_container_t operators_new; @@ -544,26 +536,28 @@ exchange_flavors(R & rng, EXTENDED_SCALAR & det, double BETA, operator_container copy_exchange_flavors_ops(creation_operators, creation_operators_new, new_flavors_first); copy_exchange_flavors_ops(annihilation_operators, annihilation_operators_new, new_flavors_first); M_TYPE M_new; - const EXTENDED_SCALAR det_new = cal_det(creation_operators_new, annihilation_operators_new, M_new, BETA, - sliding_window.get_p_model()->get_F()); - - const bool isnan_tmp = my_isnan(det_new); - if (isnan_tmp) { - std::cerr << "Warning: determinant of a new configuration is NaN. This may be because BETA is too large (overflow in computing determinant)."; + //const EXTENDED_SCALAR det_new = cal_det(creation_operators_new, annihilation_operators_new, M_new, BETA, + //sliding_window.get_p_model()->get_F()); + const std::vector& det_vec_new = cal_det_as_vector(creation_operators_new, annihilation_operators_new, M_new, BETA, + sliding_window.get_p_model()->get_F()); + SCALAR det_rat = 1.0; + for (int i=0; i +std::vector +cal_det_as_vector(const O& creation_operators, const O& annihilation_operators, MAT& M, double BETA, const G& F) { + typedef typename MAT::type SCALAR; + + const int size1 = creation_operators.size(); + + if (size1 == 0) { + throw std::runtime_error("size is zero in cal_det_as_vector"); + } + + assert(creation_operators.size() == annihilation_operators.size()); + + Eigen::Matrix matrix(size1, size1); + construct_blas_matrix(matrix, creation_operators, annihilation_operators, BETA, F); + + Eigen::FullPivLU > lu(matrix); + std::vector results(size1); + for (int i = 0; i < size1; ++i) { + results[i] = lu.matrixLU()(i,i); + } + results[0] *= lu.permutationP().determinant()*lu.permutationQ().determinant(); + + M.destructive_resize(size1, size1); + if (size1 > 0) { + M.block() = matrix; + M.safe_invert(); + } + + return results; +} + template int copy_exchange_flavors_ops(const O& operators, O& operators_new, Iterator new_flavors_begin) { int count = 0; diff --git a/src/util.hpp b/src/util.hpp index 112e733..a1389f8 100644 --- a/src/util.hpp +++ b/src/util.hpp @@ -173,11 +173,11 @@ double spectral_norm_diag(const M& mat) { Eigen::SelfAdjointEigenSolver esolv(mat_tmp,false); //Eigen::Matrix abs_evals = esolv.eigenvalues().cwiseAbs(); const double norm = std::sqrt(esolv.eigenvalues().cwiseAbs().maxCoeff())/coeff; - if (isnan(norm)) { - std::cout << "debug " << max_abs << std::endl; - std::cout << mat_tmp; - } - assert(!isnan(norm)); + //if (isnan(norm)) { + //std::cout << "debug " << max_abs << std::endl; + //std::cout << mat_tmp; + //} + //assert(!isnan(norm)); return norm; } } diff --git a/src/wide_scalar.hpp b/src/wide_scalar.hpp index 119b28e..45ff024 100644 --- a/src/wide_scalar.hpp +++ b/src/wide_scalar.hpp @@ -43,13 +43,18 @@ inline std::complex convert_to_complex(const std::complex& x) { #else -#include +#if (!defined(__clang__) && __GNUG__>0) + #include + typedef boost::multiprecision::float128 EXTENDED_REAL; +#elif (__INTEL_COMPILER>0) + #include + typedef boost::multiprecision::float128 EXTENDED_REAL; +#else + #include + typedef boost::multiprecision::cpp_bin_float_quad EXTENDED_REAL; +#endif template class wcomplex; - -//typedef boost::multiprecision::number > cpp_dec_float_15; -//typedef cpp_dec_float_15 EXTENDED_REAL; -typedef boost::multiprecision::cpp_bin_float_quad EXTENDED_REAL; typedef wcomplex EXTENDED_COMPLEX; template