From 4d0bb56790427fae6774d0fc245e18d769877944 Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Sun, 29 Sep 2013 21:05:17 +0200 Subject: [PATCH] arrays. Clean expression templates - clean array, matrix, vector expression template they take const & of objects, or move && objects no more views. -> C++11 modernisation - Fix a bug in array resize : it was resetting the indexmap to C memory layout e.g. for a fortran array - Fix a bug in read h5 array when not in C order (forgot an else, the array was read twice). --- test/triqs/arrays/expr_matrix.cpp | 6 +- test/triqs/arrays/expr_matrix.output | 6 +- test/triqs/arrays/inverse2.cpp | 12 +-- .../expression_template/array_algebra.hpp | 95 +++++++++++-------- .../expression_template/matrix_algebra.hpp | 59 +++++++----- triqs/arrays/expression_template/tools.hpp | 54 ++++++----- .../expression_template/vector_algebra.hpp | 54 ++++++----- triqs/arrays/h5/simple_read_write.hpp | 7 +- triqs/arrays/impl/indexmap_storage_pair.hpp | 3 +- triqs/arrays/linalg/det_and_inverse.hpp | 42 ++++---- 10 files changed, 194 insertions(+), 144 deletions(-) diff --git a/test/triqs/arrays/expr_matrix.cpp b/test/triqs/arrays/expr_matrix.cpp index 31151d58..591249ba 100644 --- a/test/triqs/arrays/expr_matrix.cpp +++ b/test/triqs/arrays/expr_matrix.cpp @@ -72,15 +72,15 @@ int main(int argc, char **argv) { std::cout<<" matrix( Af * (Bf + Cf) )"<(Af*(Bf+ Cf))< V(3); V(0) = 1; V(1)= 2; V(2) = 3; tqa::vector V2(3); V2(0) = 10; V2(1)= 20; V2(2) = 30; - TEST (tqa::make_vector( V2 + 2.0 *V)); + TEST (make_vector( V2 + 2.0 *V)); // test the division by a matrix TEST(Af); diff --git a/test/triqs/arrays/expr_matrix.output b/test/triqs/arrays/expr_matrix.output index ac1081c1..47016006 100644 --- a/test/triqs/arrays/expr_matrix.output +++ b/test/triqs/arrays/expr_matrix.output @@ -38,19 +38,19 @@ Cuboid of rank 2 and dimensions (2 2) (A) ---> [[0,1] [10,11]] -(tqa::make_matrix( 2*A )) ---> +(make_matrix( 2*A )) ---> [[0,2] [20,22]] (A+ 2) ---> ( [[0,1] [10,11]] + 2) -(tqa::make_matrix(A+2 )) ---> +(make_matrix(A+2 )) ---> [[2,1] [10,13]] (make_matrix(1 + A )) ---> [[1,1] [10,12]] -(tqa::make_vector( V2 + 2.0 *V)) ---> [12,24,36] +(make_vector( V2 + 2.0 *V)) ---> [12,24,36] (Af) ---> [[0,1] [10,11]] diff --git a/test/triqs/arrays/inverse2.cpp b/test/triqs/arrays/inverse2.cpp index c5174446..2e6a27cf 100644 --- a/test/triqs/arrays/inverse2.cpp +++ b/test/triqs/arrays/inverse2.cpp @@ -30,14 +30,15 @@ using std::cout; using std::endl; using namespace triqs::arrays; template -matrix_view -eval_as_matrix( Expr const & e) { return matrix(e);} +matrix +eval_as_matrix( Expr const & e) { return e;} int main(int argc, char **argv) { try { triqs::arrays::matrix W(3,3,FORTRAN_LAYOUT),Wi(3,3,FORTRAN_LAYOUT),Wkeep(3,3,FORTRAN_LAYOUT),A(FORTRAN_LAYOUT); + //triqs::arrays::matrix W(3,3),Wi(3,3),Wkeep(3,3),A{}; for (int i =0; i<3; ++i) for (int j=0; j<3; ++j) W(i,j) = (i>j ? i+2.5*j : i*0.8-j); @@ -64,9 +65,10 @@ int main(int argc, char **argv) { assert ( (abs(should_be_one(i,j) - (i==j ? 1 : 0))) <1.e-10 ); std::cerr<< "W* inverse(W)" << " = "<< triqs::arrays::matrix (W*Wi)<::value> struct get_call_const_return_type; + template struct get_call_const_return_type { typedef decltype(std::declval()()) type; }; + template struct get_call_const_return_type { typedef decltype(std::declval()(0)) type; }; + template struct get_call_const_return_type { typedef decltype(std::declval()(0,0)) type; }; + template struct get_call_const_return_type { typedef decltype(std::declval()(0,0,0)) type; }; + template struct get_call_const_return_type { typedef decltype(std::declval()(0,0,0,0)) type; }; + template struct get_call_const_return_type { typedef decltype(std::declval()(0,0,0,0,0)) type; }; + template struct array_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) { - typedef typename keeper_type::type L_t; - typedef typename keeper_type::type R_t; + typedef typename std::remove_reference::type L_t; + typedef typename std::remove_reference::type R_t; static_assert( get_rank::value==0 || get_rank::value==0 || get_rank::value == get_rank::value, "rank mismatch in array operations"); - typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; + typedef typename std::result_of(typename get_call_const_return_type::type,typename get_call_const_return_type::type)>::type value_type; + //typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; typedef typename std::remove_reference::type>::type domain_type; - - L_t l; R_t r; + + L l; R r; template array_expr(LL && l_, RR && r_) : l(std::forward(l_)), r(std::forward(r_)) {} - + domain_type domain() const { return combine_domain()(l,r); } - - //template value_type operator[](KeyType && key) const { return operation()(l[std::forward(key)] , r[std::forward(key)]);} - template value_type operator()(Args && ... args) const { return operation()(l(std::forward(args)...) , r(std::forward(args)...));} - friend std::ostream &operator <<(std::ostream &sout, array_expr const &expr){return sout << "("<::name << " "< auto operator()(Args && ... args) const DECL_AND_RETURN( utility::operation()(l(std::forward(args)...) , r(std::forward(args)...))); + template value_type operator()(Args && ... args) const { return utility::operation()(l(std::forward(args)...) , r(std::forward(args)...));} + + friend std::ostream &operator <<(std::ostream &sout, array_expr const &expr){return sout << "("<::name << " "< make_array(array_expr const & e) { return e;} }; // a special case : the unary operator ! + template + struct array_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) { + typedef typename std::remove_reference::type L_t; + typedef typename L_t::value_type value_type; + typedef typename L_t::domain_type domain_type; + + L l; + template array_unary_m_expr(LL && l_) : l(std::forward(l_)) {} + + domain_type domain() const { return l.domain(); } + template value_type operator()(Args && ... args) const { return -l(std::forward(args)...);} + + friend std::ostream &operator <<(std::ostream &sout, array_unary_m_expr const &expr){return sout << '-'< make_array(array_unary_m_expr const & e) { return e;} + }; - template struct array_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) { - typedef typename keeper_type::type L_t; - typedef typename L_t::value_type value_type; - typedef typename L_t::domain_type domain_type; - - L_t l; - template array_unary_m_expr(LL && l_) : l(std::forward(l_)) {} - - domain_type domain() const { return l.domain(); } - - //template value_type operator[](KeyType&& key) const {return -l[key];} - template value_type operator()(Args && ... args) const { return -l(std::forward(args)...);} - friend std::ostream &operator <<(std::ostream &sout, array_unary_m_expr const &expr){return sout << '-'<\ - typename std::enable_if::value && TRAIT2 ::value, array_expr>::type\ - operator OP (A1 const & a1, A2 const & a2) { return array_expr(a1,a2);} + typename std::enable_if::value && TRAIT2 ::value, \ + array_expr::type, typename node_t::type>>::type\ + operator OP (A1 && a1, A2 && a2) { return {std::forward(a1),std::forward(a2)};} -DEFINE_OPERATOR(plus, +, ImmutableArray,ImmutableArray); -DEFINE_OPERATOR(minus, -, ImmutableArray,ImmutableArray); -DEFINE_OPERATOR(multiplies, *, ImmutableArray,ImmutableArray); -DEFINE_OPERATOR(multiplies, *, is_in_ZRC,ImmutableArray); -DEFINE_OPERATOR(multiplies, *, ImmutableArray,is_in_ZRC); -DEFINE_OPERATOR(divides, /, ImmutableArray,ImmutableArray); -DEFINE_OPERATOR(divides, /, is_in_ZRC,ImmutableArray); -DEFINE_OPERATOR(divides, /, ImmutableArray,is_in_ZRC); + DEFINE_OPERATOR(plus, +, ImmutableArray,ImmutableArray); + DEFINE_OPERATOR(minus, -, ImmutableArray,ImmutableArray); + DEFINE_OPERATOR(multiplies, *, ImmutableArray,ImmutableArray); + DEFINE_OPERATOR(multiplies, *, is_in_ZRC,ImmutableArray); + DEFINE_OPERATOR(multiplies, *, ImmutableArray,is_in_ZRC); + DEFINE_OPERATOR(divides, /, ImmutableArray,ImmutableArray); + DEFINE_OPERATOR(divides, /, is_in_ZRC,ImmutableArray); + DEFINE_OPERATOR(divides, /, ImmutableArray,is_in_ZRC); #undef DEFINE_OPERATOR // the unary is special - template typename std::enable_if::value, array_unary_m_expr>::type - operator - (A1 const & a1) { return {a1};} - - template array - make_array( Expr const & e) { return array(e);} + template + typename std::enable_if< + ImmutableArray::value, + array_unary_m_expr::type > + >::type + operator - (A1 && a1) { return {std::forward(a1)};} }}//namespace triqs::arrays #endif diff --git a/triqs/arrays/expression_template/matrix_algebra.hpp b/triqs/arrays/expression_template/matrix_algebra.hpp index c8f58208..6b043902 100644 --- a/triqs/arrays/expression_template/matrix_algebra.hpp +++ b/triqs/arrays/expression_template/matrix_algebra.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -62,37 +62,34 @@ namespace triqs { namespace arrays { blas::gemv(1.0,m,v,0.0,R); return R; } - + // expression template - template + template struct matrix_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { - typedef typename keeper_type::type L_t; - typedef typename keeper_type::type R_t; + typedef typename std::remove_reference::type L_t; + typedef typename std::remove_reference::type R_t; static_assert( get_rank::value==0 || get_rank::value==0 || get_rank::value == get_rank::value, "rank mismatch in matrix operations"); - typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; + typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; typedef typename std::remove_cv< typename std::remove_reference::type>::type>::type domain_type; - L_t l; R_t r; + L l; R r; template matrix_expr(LL && l_, RR && r_) : l(std::forward(l_)), r(std::forward(r_)) {} domain_type domain() const { return combine_domain()(l,r); } - - //template value_type operator[](KeyType && key) const { return operation()(l[std::forward(key)] , r[std::forward(key)]);} - template value_type operator()(Args && ... args) const { return operation()(l(std::forward(args)...) , r(std::forward(args)...));} - friend std::ostream &operator <<(std::ostream &sout, matrix_expr const &expr){return sout << "("<::name << " "< value_type operator()(Args && ... args) const { return utility::operation()(l(std::forward(args)...) , r(std::forward(args)...));} + friend std::ostream &operator <<(std::ostream &sout, matrix_expr const &expr){return sout << "("<::name << " "< // a special case : the unary operator ! struct matrix_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { - typedef typename keeper_type::type L_t; + typedef typename std::remove_reference::type L_t; typedef typename L_t::value_type value_type; typedef typename L_t::domain_type domain_type; - L_t l; + + L l; template matrix_unary_m_expr(LL && l_) : l(std::forward(l_)) {} domain_type domain() const { return l.domain(); } - - //template value_type operator[](KeyType&& key) const {return -l[key];} template value_type operator()(Args && ... args) const { return -l(std::forward(args)...);} friend std::ostream &operator <<(std::ostream &sout, matrix_unary_m_expr const &expr){return sout << '-'<\ - typename std::enable_if::value && TRAIT2 ::value, matrix_expr>::type\ - operator OP (A1 const & a1, A2 const & a2) { return matrix_expr(a1,a2);} + typename std::enable_if::value && TRAIT2 ::value, \ + matrix_expr::type, typename node_t::type>>::type\ + operator OP (A1 && a1, A2 && a2) { return {std::forward(a1),std::forward(a2)};} DEFINE_OPERATOR(plus, +, ImmutableMatrix,ImmutableMatrix); DEFINE_OPERATOR(minus, -, ImmutableMatrix,ImmutableMatrix); @@ -115,8 +113,9 @@ namespace triqs { namespace arrays { // the addition/substraction of diagonal matrix is special : all scalar are diagonal matrices here... #define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \ template\ - typename std::enable_if::value && TRAIT2 ::value, matrix_expr>::type\ - operator OP (A1 const & a1, A2 const & a2) { return matrix_expr(a1,a2);} + typename std::enable_if::value && TRAIT2 ::value, \ + matrix_expr::type, typename node_t::type>>::type\ + operator OP (A1 && a1, A2 && a2) { return {std::forward(a1),std::forward(a2)};} DEFINE_OPERATOR(plus, +, ImmutableMatrix,is_in_ZRC); DEFINE_OPERATOR(plus, +, is_in_ZRC,ImmutableMatrix); @@ -125,13 +124,27 @@ namespace triqs { namespace arrays { #undef DEFINE_OPERATOR // the unary is special - template typename std::enable_if::value, matrix_unary_m_expr>::type - operator - (A1 const & a1) { return {a1};} + template + typename std::enable_if< + ImmutableMatrix::value, + matrix_unary_m_expr::type > + >::type + operator - (A1 && a1) { return {std::forward(a1)};} + template::value>::type > + struct _a_div_matrix { + template auto operator() (A && a, M && m) DECL_AND_RETURN ( std::forward(a) * inverse(std::forward(m))); + }; + //typedef decltype ( std::declval::type>() * inverse(std::declval::type>() )) type; + template // anything / matrix ---> anything * inverse(matrix) - typename boost::lazy_enable_if< ImmutableMatrix, type_of_mult > >::type - operator/ (A const & a, M const & m) { return a * inverse(m);} + //typename boost::lazy_enable_if< ImmutableMatrix, std::result_of<_a_div_matrix(A,M)__type_of_mult_expr_matrix >::type + //typename std::result_of<_a_div_matrix(A,M)>::type + //operator/ (A && a, M && m) {return _a_div_matrix()( std::forward(a), std::forward(m));} + auto operator/ (A && a, M && m) DECL_AND_RETURN( _a_div_matrix()( std::forward(a), std::forward(m))); + // -> typename std::enable_if< ImmutableMatrix::value, decltype(std::forward(a) * inverse(std::forward(m)))>::type + //{ return std::forward(a) * inverse(std::forward(m));} }}//namespace triqs::arrays #endif diff --git a/triqs/arrays/expression_template/tools.hpp b/triqs/arrays/expression_template/tools.hpp index 86675b9b..eef32809 100644 --- a/triqs/arrays/expression_template/tools.hpp +++ b/triqs/arrays/expression_template/tools.hpp @@ -21,8 +21,11 @@ #ifndef TRIQS_ARRAYS_EXPRESSION_TOOLS_H #define TRIQS_ARRAYS_EXPRESSION_TOOLS_H #include +#include namespace triqs { namespace arrays { + using utility::is_in_ZRC; + /* namespace tags { struct plus{}; struct minus{}; struct multiplies{}; struct divides{}; } // The basic operations put in a template.... @@ -48,35 +51,44 @@ namespace triqs { namespace arrays { template struct is_in_ZRC : boost::is_arithmetic {}; template<> struct is_in_ZRC : mpl::true_ {}; template struct is_in_ZRC > : mpl::true_ {}; +*/ // Wrapping the scalar in a little struct to recognize it // In a matrix expression, the evaluation differ (it is a scalar matrix ....) - template struct scalar_wrap; + template struct _scalar_wrap; // First the scalar for an array expression ... - template struct scalar_wrap { - typedef S value_type; - S s; scalar_wrap(S const &s_):s(s_){} - template value_type operator[](KeyType&& key) const {return s;} - template inline value_type operator()(Args && ... args) const { return s;} - friend std::ostream &operator <<(std::ostream &sout, scalar_wrap const &expr){return sout << expr.s; } + template struct _scalar_wrap { + typedef typename std::remove_reference::type value_type; + S s; + template _scalar_wrap(T && x):s(std::forward(x)){} + template S operator[](KeyType&& key) const {return s;} + template inline S operator()(Args && ... args) const { return s;} + friend std::ostream &operator <<(std::ostream &sout, _scalar_wrap const &expr){return sout << expr.s; } }; // Second the scalar for a matrix expression ... - template struct scalar_wrap { - typedef S value_type; - S s; scalar_wrap(S const &s_):s(s_){} - template value_type operator[](KeyType&& key) const {return (key[0]==key[1] ? s : S());} - template inline value_type operator()(A1 && a1, A2 && a2) const { return (a1==a2 ? s : S());} - friend std::ostream &operator <<(std::ostream &sout, scalar_wrap const &expr){return sout << expr.s; } + template struct _scalar_wrap { + static_assert(!std::is_same::value, "lll"); + typedef typename std::remove_reference::type value_type; + S s; value_type zero; + template _scalar_wrap(T && x):s(std::forward(x)), zero{} {} + template S operator[](KeyType&& key) const {return (key[0]==key[1] ? s : S());} + template inline S operator()(A1 const & a1, A2 const & a2) const { return (a1==a2 ? s : zero);} + friend std::ostream &operator <<(std::ostream &sout, _scalar_wrap const &expr){return sout << expr.s; } }; + // type of the node + template struct node_t : + std::conditional::value, _scalar_wrap::type,IsMatrix>, typename utility::remove_rvalue_ref::type> {}; + //std::conditional::value, _scalar_wrap::type,IsMatrix>, typename utility::remove_rvalue_ref::type> {}; + // get the rank of something .... template struct get_rank { static constexpr int value = std::remove_reference::type::domain_type::rank;}; - template struct get_rank> { static constexpr int value =0;}; + template struct get_rank<_scalar_wrap> { static constexpr int value =0;}; // - template struct keeper_type : boost::mpl::if_, scalar_wrap, typename view_type_if_exists_else_type::type> {}; + //template struct keeper_type : boost::mpl::if_, _scalar_wrap, typename view_type_if_exists_else_type::type> {}; // Combine the two domains of LHS and RHS : need to specialize where there is a scalar struct combine_domain { @@ -85,16 +97,10 @@ namespace triqs { namespace arrays { if (l.domain().lengths() != r.domain().lengths()) TRIQS_RUNTIME_ERROR << "Domain size mismatch : "<< l.domain().lengths()<<" vs" < auto operator() (scalar_wrap const & w, R const & r) const -> decltype(r.domain()) { return r.domain();} - template auto operator() (L const & l, scalar_wrap const & w) const -> decltype(l.domain()) { return l.domain();} + template auto operator() (_scalar_wrap const & w, R const & r) const -> decltype(r.domain()) { return r.domain();} + template auto operator() (L const & l, _scalar_wrap const & w) const -> decltype(l.domain()) { return l.domain();} }; - template - struct type_of_mult{ - //typedef typename boost::remove_reference::type T1; typedef typename boost::remove_reference::type T2; - //typedef BOOST_TYPEOF_TPL( (*(pseudo_default_construct())) * (*(pseudo_default_construct()))) type; - typedef decltype ( std::declval::type>() * std::declval::type>() ) type; - }; - + }}//namespace triqs::arrays #endif diff --git a/triqs/arrays/expression_template/vector_algebra.hpp b/triqs/arrays/expression_template/vector_algebra.hpp index 235bf291..0bfc7821 100644 --- a/triqs/arrays/expression_template/vector_algebra.hpp +++ b/triqs/arrays/expression_template/vector_algebra.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -26,49 +26,54 @@ namespace triqs { namespace arrays { template struct vector_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableVector) { - typedef typename keeper_type::type L_t; - typedef typename keeper_type::type R_t; - L_t l; R_t r; - template vector_expr(LL && l_, RR && r_) : l(std::forward(l_)), r(std::forward(r_)) {} + typedef typename std::remove_reference::type L_t; + typedef typename std::remove_reference::type R_t; static_assert( get_rank::value==0 || get_rank::value==0 || get_rank::value == get_rank::value, "rank mismatch in array operations"); - typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; + typedef typename std::result_of(typename L_t::value_type,typename R_t::value_type)>::type value_type; typedef typename std::remove_reference::type>::type domain_type; + + L l; R r; + template vector_expr(LL && l_, RR && r_) : l(std::forward(l_)), r(std::forward(r_)) {} domain_type domain() const { return combine_domain()(l,r); } - mini_vector shape() const { return this->domain().lengths();} + //mini_vector shape() const { return this->domain().lengths();} size_t size() const { return this->domain().lengths()[0];} - //template value_type operator[](KeyType && key) const { return operation()(l[std::forward(key)] , r[std::forward(key)]);} - template value_type operator()(Args && ... args) const { return operation()(l(std::forward(args)...) , r(std::forward(args)...));} - friend std::ostream &operator <<(std::ostream &sout, vector_expr const &expr){return sout << "("<::name << " "< value_type operator[](KeyType && key) const { return utility::operation()(l[std::forward(key)] , r[std::forward(key)]);} + template value_type operator()(Args && ... args) const { return utility::operation()(l(std::forward(args)...) , r(std::forward(args)...));} + friend std::ostream &operator <<(std::ostream &sout, vector_expr const &expr){return sout << "("<::name << " "< make_vector(vector_expr const & e) { return e;} }; template // a special case : the unary operator ! struct vector_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableVector) { - typedef typename keeper_type::type L_t; - L_t l; - template vector_unary_m_expr(LL && l_) : l(std::forward(l_)) {} + typedef typename std::remove_reference::type L_t; typedef typename L_t::value_type value_type; typedef typename L_t::domain_type domain_type; + L l; + template vector_unary_m_expr(LL && l_) : l(std::forward(l_)) {} + domain_type domain() const { return l.domain(); } - mini_vector shape() const { return this->domain().lengths();} + //mini_vector shape() const { return this->domain().lengths();} size_t size() const { return this->domain().lengths()[0];} - //template value_type operator[](KeyType&& key) const {return -l[key];} + template value_type operator[](KeyType&& key) const {return -l[std::forward(key)];} template value_type operator()(Args && ... args) const { return -l(std::forward(args)...);} friend std::ostream &operator <<(std::ostream &sout, vector_unary_m_expr const &expr){return sout << '-'< make_vector(vector_unary_m_expr const & e) { return e;} }; // Now we can define all the C++ operators ... #define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \ template\ - typename std::enable_if::value && TRAIT2 ::value, vector_expr>::type\ - operator OP (A1 const & a1, A2 const & a2) { return vector_expr(a1,a2);} + typename std::enable_if::value && TRAIT2 ::value,\ + vector_expr::type, typename node_t::type>>::type\ + operator OP (A1 && a1, A2 && a2) { return {std::forward(a1),std::forward(a2)};} #define DELETE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \ template\ - typename std::enable_if::value && TRAIT2 ::value, vector_expr>::type\ - operator OP (A1 const & a1, A2 const & a2) = delete; + typename std::enable_if::value && TRAIT2 ::value, vector_expr>::type\ + operator OP (A1 && a1, A2 && a2) = delete; DEFINE_OPERATOR(plus, +, ImmutableVector,ImmutableVector); DEFINE_OPERATOR(minus, -, ImmutableVector,ImmutableVector); DELETE_OPERATOR(multiplies, *, ImmutableVector,ImmutableVector); @@ -79,11 +84,12 @@ namespace triqs { namespace arrays { #undef DEFINE_OPERATOR // the unary is special - template typename std::enable_if::value, vector_unary_m_expr>::type - operator - (A1 const & a1) { return {a1};} - - template vector - make_vector( Expr const & e) { return vector(e);} + template + typename std::enable_if< + ImmutableVector::value, + vector_unary_m_expr::type > + >::type + operator - (A1 && a1) { return {std::forward(a1)};} }}//namespace triqs::arrays #endif diff --git a/triqs/arrays/h5/simple_read_write.hpp b/triqs/arrays/h5/simple_read_write.hpp index f1d0be3a..5d84e255 100644 --- a/triqs/arrays/h5/simple_read_write.hpp +++ b/triqs/arrays/h5/simple_read_write.hpp @@ -116,11 +116,12 @@ namespace triqs { namespace arrays { mini_vector dims_out; //int ndims = dataspace.getSimpleExtentDims( &dims_out[0], NULL); dataspace.getSimpleExtentDims( &dims_out[0], NULL); - mini_vector d2; for (size_t u=0; u d2; + for (size_t u=0; u(A).view(),false);} //if (C_reorder) { read_array(g,name, make_cache(A).view(),false);} - ds.read( get_array_data_ptr(A), h5::data_type_memory(), data_space(A) , dataspace ); + else ds.read( get_array_data_ptr(A), h5::data_type_memory(), data_space(A) , dataspace ); } TRIQS_ARRAYS_H5_CATCH_EXCEPTION; } diff --git a/triqs/arrays/impl/indexmap_storage_pair.hpp b/triqs/arrays/impl/indexmap_storage_pair.hpp index 3a2dcdcf..ea6aabcf 100644 --- a/triqs/arrays/impl/indexmap_storage_pair.hpp +++ b/triqs/arrays/impl/indexmap_storage_pair.hpp @@ -319,7 +319,8 @@ namespace triqs { namespace arrays { // ------------------------------- resize -------------------------------------------- // void resize (domain_type const & d) { - this->indexmap_ = IndexMapType(d);// build a new one with the lengths of IND + this->indexmap_ = IndexMapType(d,this->indexmap_.memory_indices_layout()); + // build a new one with the lengths of IND BUT THE SAME layout ! // optimisation. Construct a storage only if the new index is not compatible (size mismatch). if (this->storage_.size() != this->indexmap_.domain().number_of_elements()) this->storage_ = StorageType(this->indexmap_.domain().number_of_elements(), typename flags::init_tag::type() ); diff --git a/triqs/arrays/linalg/det_and_inverse.hpp b/triqs/arrays/linalg/det_and_inverse.hpp index 062cc0e0..b4f41475 100644 --- a/triqs/arrays/linalg/det_and_inverse.hpp +++ b/triqs/arrays/linalg/det_and_inverse.hpp @@ -43,12 +43,13 @@ namespace triqs { namespace arrays { * It keeps view of the object A if it a matrix, a copy if it is a formal expression. */ template struct inverse_lazy; + template struct inverse_lazy_impl;//debug ony /// template struct determinant_lazy; /// Lazy inversion - template inverse_lazy inverse (A const & a) { return inverse_lazy(a); } + template inverse_lazy::type> inverse (A && a) { return {std::forward(a)}; } /// Lazy computation of det template determinant_lazy determinant (A const & a) { return determinant_lazy(a); } @@ -111,58 +112,67 @@ namespace triqs { namespace arrays { // an implementation class to gather the common part to matrix and expression.... template struct inverse_lazy_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { - typedef typename std::remove_const::type value_type; - typedef typename A::domain_type domain_type; - typedef typename const_view_type_if_exists_else_type::type A_type; - const A_type a; - inverse_lazy_impl(A const & a_):a (a_) { + public: + + typedef typename std::remove_reference::type A_t; + typedef typename std::remove_const::type value_type; + typedef typename A_t::domain_type domain_type; + A a; + template inverse_lazy_impl(AA && a_):a (std::forward(a_)) { if (first_dim(a) != second_dim(a)) TRIQS_RUNTIME_ERROR<< "Inverse : matrix is not square but of size "<< first_dim(a)<<" x "<< second_dim(a); } - //typename A::shape_type shape() const { return a.shape();} domain_type domain() const { return a.domain(); } - template value_type operator() (K0 const & k0, K1 const & k1) const { activate(); return _id->M(k0,k1); } + template value_type operator() (K0 const & k0, K1 const & k1) const { activate(); return _id->M(k0,k1); } friend std::ostream & operator<<(std::ostream & out,inverse_lazy_impl const&x){return out<<"inverse("< M_type; typedef matrix_view M_view_type; M_type M; - internal_data(inverse_lazy_impl const & P):M(P.a){det_and_inverse_worker worker(M); worker.inverse();} + internal_data(inverse_lazy_impl const & P):M(P.a){ + det_and_inverse_worker worker(M); + worker.inverse(); + } }; friend struct internal_data; mutable std::shared_ptr _id; - void activate() const { if (!_id) _id= std::make_shared(*this);} + void activate() const { if (!_id) _id= std::make_shared(*this);} }; // The general case template - struct inverse_lazy >::type > : inverse_lazy_impl { - inverse_lazy(A const & a):inverse_lazy_impl(a) { } + struct inverse_lazy::type>) > : inverse_lazy_impl { + template inverse_lazy(AA && a_):inverse_lazy_impl(std::forward(a_)) {} }; // for matrix and matrix_views, we have more optimisation possible .... template - struct inverse_lazy >::type >:inverse_lazy_impl{ - inverse_lazy(A const & a):inverse_lazy_impl(a) { } + struct inverse_lazy::type>) >:inverse_lazy_impl{ + template inverse_lazy(AA && a_):inverse_lazy_impl(std::forward(a_)) {} template // Optimized implementation of = friend void triqs_arrays_assign_delegation (MT & lhs, inverse_lazy const & rhs) { + //std::cerr << " DELEGATING1"<< lhs <::value, "Internal error"); + //std::cerr << " DELEGATING"<< lhs << std::endl; if ((lhs.indexmap().memory_indices_layout() !=rhs.a.indexmap().memory_indices_layout())|| (lhs.data_start() != rhs.a.data_start()) || !(has_contiguous_data(lhs))) { rhs.activate(); lhs = rhs._id->M;} else {// if M = inverse(M) with the SAME object, then we do not need to copy the data + //std::cerr << " DELEGATING"<< lhs << std::endl; blas_lapack_tools::reflexive_qcache C(lhs);// a reflexive cache will use a temporary "regrouping" copy if and only if needed det_and_inverse_worker W(C());// the worker to make the inversion of the lhs... W.inverse(); // worker is working ... + //std::cerr << " DELEGATING"<< lhs << std::endl; } + //std::cerr << " okok "<< lhs << std::endl; } friend std::ostream & operator<<(std::ostream & out,inverse_lazy const&x){return out<<"inverse("< struct determinant_lazy { // : { Tag::expression_terminal, Tag::scalar_expression_terminal { + template struct determinant_lazy { typedef typename A::value_type value_type; typedef typename const_view_type_if_exists_else_type::type A_type; A_type a;