diff --git a/test/triqs/gfs/bug1.cpp b/test/triqs/gfs/bug1.cpp new file mode 100644 index 00000000..ce94267c --- /dev/null +++ b/test/triqs/gfs/bug1.cpp @@ -0,0 +1,44 @@ +#include +#include "./common.hpp" + +using namespace triqs::gfs; + +int main() { + + using A = triqs::arrays::matrix_tensor_proxy, 3, void> &, true>; + + static_assert(std::is_constructible, matsubara_freq>::value, "oops"); + static_assert(triqs::arrays::is_scalar_or_convertible::value, "oops2"); + static_assert(triqs::arrays::is_scalar_for, A>::value, "oops2"); + static_assert(triqs::arrays::is_scalar_for::value, "oops2"); + + triqs::clef::placeholder<0> om_; + auto g = gf{{10, Fermion, 10}, {2, 2}}; + auto g2 = g; + auto g3 = g; + + g(om_) << om_ + 0.0; // Works + + g2(om_) << om_; // Did not compile : rhs = mesh_point + g3(om_) << om_ + om_; // Did not compile : rhs = matsubara_freq + + assert_equal_array(g.data(), g2.data(), "bug !"); + assert_equal_array(g.data(), g3.data()/2, "bug !"); + + std::cerr << g.data()(triqs::arrays::ellipsis(), 0,0) << std::endl; + std::cerr << g.data()(triqs::arrays::ellipsis(), 1,0) << std::endl; + std::cerr << g.data()(triqs::arrays::ellipsis(), 0,1) << std::endl; + std::cerr << g.data()(triqs::arrays::ellipsis(), 1,1) << std::endl; + + std::cerr << g2.data()(triqs::arrays::ellipsis(), 0,0) << std::endl; + std::cerr << g2.data()(triqs::arrays::ellipsis(), 1,0) << std::endl; + std::cerr << g2.data()(triqs::arrays::ellipsis(), 0,1) << std::endl; + std::cerr << g2.data()(triqs::arrays::ellipsis(), 1,1) << std::endl; + + std::cerr << g3.data()(triqs::arrays::ellipsis(), 0,0) << std::endl; + std::cerr << g3.data()(triqs::arrays::ellipsis(), 1,0) << std::endl; + std::cerr << g3.data()(triqs::arrays::ellipsis(), 0,1) << std::endl; + std::cerr << g3.data()(triqs::arrays::ellipsis(), 1,1) << std::endl; + + return 0; +} diff --git a/test/triqs/gfs/common.hpp b/test/triqs/gfs/common.hpp new file mode 100644 index 00000000..8b60231e --- /dev/null +++ b/test/triqs/gfs/common.hpp @@ -0,0 +1,11 @@ +template +void assert_equal(T const& x, T const& y, std::string mess) { + if (std::abs(x - y) > 1.e-13) TRIQS_RUNTIME_ERROR << mess; +} + +template +void assert_equal_array(T1 const& x, T2 const& y, std::string mess) { + if (max_element(abs(x - y)) > 1.e-13) TRIQS_RUNTIME_ERROR << mess << "\n" << x << "\n" << y << "\n" << max_element(abs(x - y)); +} + + diff --git a/test/triqs/gfs/multivar/curry_and_fourier.cpp b/test/triqs/gfs/multivar/curry_and_fourier.cpp index 0fec8274..674aedf3 100644 --- a/test/triqs/gfs/multivar/curry_and_fourier.cpp +++ b/test/triqs/gfs/multivar/curry_and_fourier.cpp @@ -2,6 +2,7 @@ #include #include #include +#include "../common.hpp" namespace h5 = triqs::h5; using namespace triqs::gfs; @@ -9,16 +10,6 @@ using namespace triqs::clef; using namespace triqs::arrays; using namespace triqs::lattice; -template -void assert_equal(T const& x, T const& y, std::string mess) { - if (std::abs(x - y) > 1.e-13) TRIQS_RUNTIME_ERROR << mess; -} - -template -void assert_equal_array(T1 const& x, T2 const& y, std::string mess) { - if (max_element(abs(x - y)) > 1.e-13) TRIQS_RUNTIME_ERROR << mess << "\n" << x << "\n" << y << "\n" << max_element(abs(x - y)); -} - #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; int main() { diff --git a/triqs/arrays/impl/assignment.hpp b/triqs/arrays/impl/assignment.hpp index 0d752bd7..2637daf3 100644 --- a/triqs/arrays/impl/assignment.hpp +++ b/triqs/arrays/impl/assignment.hpp @@ -135,31 +135,52 @@ namespace triqs { namespace arrays { TRIQS_REJECT_ASSIGN_TO_CONST; typedef typename LHS::value_type value_type; LHS & lhs; const RHS & rhs; - impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_){}//, p(*(lhs_.data_start())) {} + impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_){} template void operator()(Args const & ...args) const {_ops_::invoke(lhs(args...), rhs);} void invoke() { foreach(lhs,*this); } }; // ----------------- assignment for scalar RHS for Matrices -------------------------------------------------- - template bool kronecker(mini_vector const & key) { return ( (R==2) && (key[0]==key[1]));} + //template bool kronecker(mini_vector const & key) { return ( (R==2) && (key[0]==key[1]));} template bool kronecker(T const & x0, T const & x1) { return ( (x0==x1));} + // CONCEPT : reunifiy the 2 class, put require on operator() for the 2 cases // Specialisation for Matrix Classes : scalar is a unity matrix, and operation is E, A, S, but NOT M, D - template - struct impl::value && (MutableMatrix::value && (OP=='A'||OP=='S'||OP=='E')))> { - TRIQS_REJECT_ASSIGN_TO_CONST; - typedef typename LHS::value_type value_type; - LHS & lhs; const RHS & rhs; - impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_){} //, p(*(lhs_.data_start())) {} - // we MUST make off_diag like this, if value_type is a complicated type (i.e. gf, matrix) with a size - // off diagonal element is 0*rhs, i.e. a 0, but with the SAME SIZE as the diagonal part. - // otherwise further operation may fail later. - // TO DO : look at performance issue ?? (we can remote the multiplication by 0 using an auxiliary function) - template - void operator()(Args const & ... args) const {_ops_::invoke(lhs(args...), (kronecker(args...) ? rhs : RHS{0*rhs}));} - void invoke() { foreach(lhs,*this); } - }; + // First case : when it is a true scalar or convertible to + template + struct impl::value&& is_scalar_or_convertible::value&&( + MutableMatrix::value&&(OP == 'A' || OP == 'S' || OP == 'E')))> { + TRIQS_REJECT_ASSIGN_TO_CONST; + using value_type = typename LHS::value_type; + static_assert(is_scalar::value, "Internal error"); + LHS& lhs; + const RHS& rhs; + impl(LHS& lhs_, const RHS& rhs_) : lhs(lhs_), rhs(rhs_) {} + template void operator()(Args const&... args) const { + if (kronecker(args...)) + _ops_::invoke(lhs(args...), rhs); + else + _ops_::invoke(lhs(args...), 0); + } + void invoke() { foreach(lhs, *this); } + }; + + // Specialisation for Matrix Classes : scalar is a unity matrix, and operation is E, A, S, but NOT M, D + // Second generic case : we should introduce make_zero function ? + template + struct impl::value&&(!is_scalar_or_convertible::value) && + (MutableMatrix::value&&(OP == 'A' || OP == 'S' || OP == 'E')))> { + TRIQS_REJECT_ASSIGN_TO_CONST; + typedef typename LHS::value_type value_type; + LHS& lhs; + const RHS& rhs; + impl(LHS& lhs_, const RHS& rhs_) : lhs(lhs_), rhs(rhs_) {} + template void operator()(Args const&... args) const { + _ops_::invoke(lhs(args...), (kronecker(args...) ? rhs : RHS{0 * rhs})); + } + void invoke() { foreach(lhs, *this); } + }; #undef TRIQS_REJECT_MATRIX_COMPOUND_MUL_DIV_NON_SCALAR #undef TRIQS_REJECT_ASSIGN_TO_CONST diff --git a/triqs/arrays/impl/traits.hpp b/triqs/arrays/impl/traits.hpp index afea8f3d..a20c24ef 100644 --- a/triqs/arrays/impl/traits.hpp +++ b/triqs/arrays/impl/traits.hpp @@ -70,10 +70,13 @@ namespace arrays { template struct is_amv_value_or_view_class : _or, is_amv_view_class> {}; template struct is_scalar : _or, triqs::is_complex> {}; + template + struct is_scalar_or_convertible + : std::integral_constant::value || std::is_constructible, S>::value> {}; template struct is_scalar_for - : std::conditional::value, is_scalar, std::is_same>::type { + : std::conditional::value, is_scalar_or_convertible, std::is_same>::type { }; } } // namespace triqs::arrays