mirror of
https://github.com/triqs/dft_tools
synced 2024-12-25 05:43:40 +01:00
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).
This commit is contained in:
Olivier Parcollet
Michel Ferrero
parent
59288e597f
commit
4d0bb56790
@ -72,15 +72,15 @@ int main(int argc, char **argv) {
|
||||
std::cout<<" matrix( Af * (Bf + Cf) )"<<matrix<double>(Af*(Bf+ Cf))<<std::endl;
|
||||
|
||||
TEST (A);
|
||||
TEST (tqa::make_matrix( 2*A ));
|
||||
TEST (make_matrix( 2*A ));
|
||||
TEST (A+ 2);
|
||||
TEST (tqa::make_matrix(A+2 ));
|
||||
TEST (make_matrix(A+2 ));
|
||||
TEST (make_matrix(1 + A ));
|
||||
|
||||
// test the vector ops : scalar * vector, vector + vector, ...
|
||||
tqa::vector<double> V(3); V(0) = 1; V(1)= 2; V(2) = 3;
|
||||
tqa::vector<double> 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);
|
||||
|
||||
@ -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]]
|
||||
|
||||
@ -30,14 +30,15 @@ using std::cout; using std::endl;
|
||||
using namespace triqs::arrays;
|
||||
|
||||
template<typename Expr >
|
||||
matrix_view <typename Expr::value_type>
|
||||
eval_as_matrix( Expr const & e) { return matrix<typename Expr::value_type>(e);}
|
||||
matrix <typename Expr::value_type>
|
||||
eval_as_matrix( Expr const & e) { return e;}
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
|
||||
try {
|
||||
|
||||
triqs::arrays::matrix<double> W(3,3,FORTRAN_LAYOUT),Wi(3,3,FORTRAN_LAYOUT),Wkeep(3,3,FORTRAN_LAYOUT),A(FORTRAN_LAYOUT);
|
||||
//triqs::arrays::matrix<double> 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<double > (W*Wi)<<std::endl<<std::endl;
|
||||
|
||||
W= inverse(W);
|
||||
std::cout<< " invert of W= "<< W<<std::endl<<std::endl;
|
||||
|
||||
|
||||
A= inverse(W);
|
||||
std::cerr<< " A = inv W= "<< A<<std::endl<<std::endl;
|
||||
for (int i =0; i<3; ++i)
|
||||
@ -94,9 +96,7 @@ int main(int argc, char **argv) {
|
||||
|
||||
|
||||
}
|
||||
catch (std::string ERR) { std::cout<<"ERROR : "<< ERR<<std::endl;}
|
||||
|
||||
return 0;
|
||||
TRIQS_CATCH_AND_ABORT;
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -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
|
||||
@ -23,63 +23,76 @@
|
||||
#include "./tools.hpp"
|
||||
namespace triqs { namespace arrays {
|
||||
|
||||
// a trait to compute the return type of the operator()(0,...0) const for anything with the ImmutableCuboidArray concept
|
||||
template<typename A, int R = get_rank<A>::value> struct get_call_const_return_type;
|
||||
template<typename A> struct get_call_const_return_type<A,0> { typedef decltype(std::declval<A const>()()) type; };
|
||||
template<typename A> struct get_call_const_return_type<A,1> { typedef decltype(std::declval<A const>()(0)) type; };
|
||||
template<typename A> struct get_call_const_return_type<A,2> { typedef decltype(std::declval<A const>()(0,0)) type; };
|
||||
template<typename A> struct get_call_const_return_type<A,3> { typedef decltype(std::declval<A const>()(0,0,0)) type; };
|
||||
template<typename A> struct get_call_const_return_type<A,4> { typedef decltype(std::declval<A const>()(0,0,0,0)) type; };
|
||||
template<typename A> struct get_call_const_return_type<A,5> { typedef decltype(std::declval<A const>()(0,0,0,0,0)) type; };
|
||||
|
||||
template<typename Tag, typename L, typename R>
|
||||
struct array_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) {
|
||||
typedef typename keeper_type<L>::type L_t;
|
||||
typedef typename keeper_type<R>::type R_t;
|
||||
typedef typename std::remove_reference<L>::type L_t;
|
||||
typedef typename std::remove_reference<R>::type R_t;
|
||||
static_assert( get_rank<R_t>::value==0 || get_rank<L_t>::value==0 || get_rank<L_t>::value == get_rank<R_t>::value, "rank mismatch in array operations");
|
||||
typedef typename std::result_of<operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::result_of<utility::operation<Tag>(typename get_call_const_return_type<L_t>::type,typename get_call_const_return_type<R_t>::type)>::type value_type;
|
||||
//typedef typename std::result_of<utility::operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::remove_reference<typename std::result_of<combine_domain(L_t,R_t)>::type>::type domain_type;
|
||||
|
||||
L_t l; R_t r;
|
||||
|
||||
L l; R r;
|
||||
template<typename LL, typename RR> array_expr(LL && l_, RR && r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
|
||||
|
||||
|
||||
domain_type domain() const { return combine_domain()(l,r); }
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType && key) const { return operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]);}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
friend std::ostream &operator <<(std::ostream &sout, array_expr const &expr){return sout << "("<<expr.l << " "<<operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
//template<typename ... Args> auto operator()(Args && ... args) const DECL_AND_RETURN( utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...)));
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
|
||||
friend std::ostream &operator <<(std::ostream &sout, array_expr const &expr){return sout << "("<<expr.l << " "<<utility::operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
friend array<value_type, domain_type::rank> make_array(array_expr const & e) { return e;}
|
||||
};
|
||||
|
||||
// a special case : the unary operator !
|
||||
template<typename L>
|
||||
struct array_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) {
|
||||
typedef typename std::remove_reference<L>::type L_t;
|
||||
typedef typename L_t::value_type value_type;
|
||||
typedef typename L_t::domain_type domain_type;
|
||||
|
||||
L l;
|
||||
template<typename LL> array_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
|
||||
|
||||
domain_type domain() const { return l.domain(); }
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return -l(std::forward<Args>(args)...);}
|
||||
|
||||
friend std::ostream &operator <<(std::ostream &sout, array_unary_m_expr const &expr){return sout << '-'<<expr.l; }
|
||||
friend array<value_type, domain_type::rank> make_array(array_unary_m_expr const & e) { return e;}
|
||||
};
|
||||
|
||||
template<typename L> struct array_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) {
|
||||
typedef typename keeper_type<L>::type L_t;
|
||||
typedef typename L_t::value_type value_type;
|
||||
typedef typename L_t::domain_type domain_type;
|
||||
|
||||
L_t l;
|
||||
template<typename LL> array_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
|
||||
|
||||
domain_type domain() const { return l.domain(); }
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType&& key) const {return -l[key];}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return -l(std::forward<Args>(args)...);}
|
||||
friend std::ostream &operator <<(std::ostream &sout, array_unary_m_expr const &expr){return sout << '-'<<expr.l; }
|
||||
};
|
||||
|
||||
// Now we can define all the C++ operators ...
|
||||
#define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \
|
||||
template<typename A1, typename A2>\
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, array_expr<tags::TAG, A1,A2>>::type\
|
||||
operator OP (A1 const & a1, A2 const & a2) { return array_expr<tags::TAG, A1,A2>(a1,a2);}
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
|
||||
array_expr<utility::tags::TAG, typename node_t<A1,false>::type, typename node_t<A2,false>::type>>::type\
|
||||
operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(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 A1> typename std::enable_if<ImmutableArray<A1>::value, array_unary_m_expr<A1>>::type
|
||||
operator - (A1 const & a1) { return {a1};}
|
||||
|
||||
template<typename Expr > array<typename Expr::value_type, Expr::domain_type::rank>
|
||||
make_array( Expr const & e) { return array<typename Expr::value_type, Expr::domain_type::rank>(e);}
|
||||
template<typename A1>
|
||||
typename std::enable_if<
|
||||
ImmutableArray<A1>::value,
|
||||
array_unary_m_expr<typename node_t<A1,false>::type >
|
||||
>::type
|
||||
operator - (A1 && a1) { return {std::forward<A1>(a1)};}
|
||||
|
||||
}}//namespace triqs::arrays
|
||||
#endif
|
||||
|
||||
@ -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<typename Tag, typename L, typename R, bool scalar_are_diagonal_matrices= false>
|
||||
template<typename Tag, typename L, typename R>
|
||||
struct matrix_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) {
|
||||
typedef typename keeper_type<L,scalar_are_diagonal_matrices>::type L_t;
|
||||
typedef typename keeper_type<R,scalar_are_diagonal_matrices>::type R_t;
|
||||
typedef typename std::remove_reference<L>::type L_t;
|
||||
typedef typename std::remove_reference<R>::type R_t;
|
||||
static_assert( get_rank<R_t>::value==0 || get_rank<L_t>::value==0 || get_rank<L_t>::value == get_rank<R_t>::value, "rank mismatch in matrix operations");
|
||||
typedef typename std::result_of<operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::result_of<utility::operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::remove_cv< typename std::remove_reference<typename std::result_of<combine_domain(L_t,R_t)>::type>::type>::type domain_type;
|
||||
|
||||
L_t l; R_t r;
|
||||
L l; R r;
|
||||
template<typename LL, typename RR> matrix_expr(LL && l_, RR && r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
|
||||
|
||||
domain_type domain() const { return combine_domain()(l,r); }
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType && key) const { return operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]);}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
friend std::ostream &operator <<(std::ostream &sout, matrix_expr const &expr){return sout << "("<<expr.l << " "<<operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
friend std::ostream &operator <<(std::ostream &sout, matrix_expr const &expr){return sout << "("<<expr.l << " "<<utility::operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
};
|
||||
|
||||
template<typename L> // a special case : the unary operator !
|
||||
struct matrix_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) {
|
||||
typedef typename keeper_type<L>::type L_t;
|
||||
typedef typename std::remove_reference<L>::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<typename LL> matrix_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
|
||||
|
||||
domain_type domain() const { return l.domain(); }
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType&& key) const {return -l[key];}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return -l(std::forward<Args>(args)...);}
|
||||
friend std::ostream &operator <<(std::ostream &sout, matrix_unary_m_expr const &expr){return sout << '-'<<expr.l; }
|
||||
};
|
||||
@ -100,8 +97,9 @@ namespace triqs { namespace arrays {
|
||||
// Now we can define all the C++ operators ...
|
||||
#define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \
|
||||
template<typename A1, typename A2>\
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, matrix_expr<tags::TAG, A1,A2>>::type\
|
||||
operator OP (A1 const & a1, A2 const & a2) { return matrix_expr<tags::TAG, A1,A2>(a1,a2);}
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
|
||||
matrix_expr<utility::tags::TAG, typename node_t<A1,false>::type, typename node_t<A2,false>::type>>::type\
|
||||
operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(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 A1, typename A2>\
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, matrix_expr<tags::TAG, A1,A2,true>>::type\
|
||||
operator OP (A1 const & a1, A2 const & a2) { return matrix_expr<tags::TAG, A1,A2,true>(a1,a2);}
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
|
||||
matrix_expr<utility::tags::TAG, typename node_t<A1,true>::type, typename node_t<A2,true>::type>>::type\
|
||||
operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(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 A1> typename std::enable_if<ImmutableMatrix<A1>::value, matrix_unary_m_expr<A1>>::type
|
||||
operator - (A1 const & a1) { return {a1};}
|
||||
template<typename A1>
|
||||
typename std::enable_if<
|
||||
ImmutableMatrix<A1>::value,
|
||||
matrix_unary_m_expr<typename node_t<A1,true>::type >
|
||||
>::type
|
||||
operator - (A1 && a1) { return {std::forward<A1>(a1)};}
|
||||
|
||||
|
||||
template<typename M1, typename Enable = typename std::enable_if<ImmutableMatrix<M1>::value>::type >
|
||||
struct _a_div_matrix {
|
||||
template<typename A, typename M> auto operator() (A && a, M && m) DECL_AND_RETURN ( std::forward<A>(a) * inverse(std::forward<M>(m)));
|
||||
};
|
||||
//typedef decltype ( std::declval<typename std::remove_reference<A>::type>() * inverse(std::declval<typename std::remove_reference<M>::type>() )) type;
|
||||
|
||||
template<typename A, typename M> // anything / matrix ---> anything * inverse(matrix)
|
||||
typename boost::lazy_enable_if< ImmutableMatrix<M>, type_of_mult<A, inverse_lazy <M> > >::type
|
||||
operator/ (A const & a, M const & m) { return a * inverse(m);}
|
||||
//typename boost::lazy_enable_if< ImmutableMatrix<M>, std::result_of<_a_div_matrix(A,M)__type_of_mult_expr_matrix<A,M> >::type
|
||||
//typename std::result_of<_a_div_matrix(A,M)>::type
|
||||
//operator/ (A && a, M && m) {return _a_div_matrix<M>()( std::forward<A>(a), std::forward<M>(m));}
|
||||
auto operator/ (A && a, M && m) DECL_AND_RETURN( _a_div_matrix<M>()( std::forward<A>(a), std::forward<M>(m)));
|
||||
// -> typename std::enable_if< ImmutableMatrix<M>::value, decltype(std::forward<A>(a) * inverse(std::forward<M>(m)))>::type
|
||||
//{ return std::forward<A>(a) * inverse(std::forward<M>(m));}
|
||||
|
||||
}}//namespace triqs::arrays
|
||||
#endif
|
||||
|
||||
@ -21,8 +21,11 @@
|
||||
#ifndef TRIQS_ARRAYS_EXPRESSION_TOOLS_H
|
||||
#define TRIQS_ARRAYS_EXPRESSION_TOOLS_H
|
||||
#include <type_traits>
|
||||
#include <triqs/utility/expression_template_tools.hpp>
|
||||
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<typename T> struct is_in_ZRC : boost::is_arithmetic<T> {};
|
||||
template<> struct is_in_ZRC<bool> : mpl::true_ {};
|
||||
template<typename T> struct is_in_ZRC<std::complex<T> > : 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<typename S, bool IsMatrix> struct scalar_wrap;
|
||||
template<typename S, bool IsMatrix> struct _scalar_wrap;
|
||||
|
||||
// First the scalar for an array expression ...
|
||||
template<typename S> struct scalar_wrap<S,false> {
|
||||
typedef S value_type;
|
||||
S s; scalar_wrap(S const &s_):s(s_){}
|
||||
template<typename KeyType> value_type operator[](KeyType&& key) const {return s;}
|
||||
template<typename ... Args> 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<typename S> struct _scalar_wrap<S,false> {
|
||||
typedef typename std::remove_reference<S>::type value_type;
|
||||
S s;
|
||||
template<typename T> _scalar_wrap(T && x):s(std::forward<T>(x)){}
|
||||
template<typename KeyType> S operator[](KeyType&& key) const {return s;}
|
||||
template<typename ... Args> 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<typename S> struct scalar_wrap<S,true> {
|
||||
typedef S value_type;
|
||||
S s; scalar_wrap(S const &s_):s(s_){}
|
||||
template<typename KeyType> value_type operator[](KeyType&& key) const {return (key[0]==key[1] ? s : S());}
|
||||
template<typename A1, typename A2> 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<typename S> struct _scalar_wrap<S,true> {
|
||||
static_assert(!std::is_same<double &, S>::value, "lll");
|
||||
typedef typename std::remove_reference<S>::type value_type;
|
||||
S s; value_type zero;
|
||||
template<typename T> _scalar_wrap(T && x):s(std::forward<T>(x)), zero{} {}
|
||||
template<typename KeyType> S operator[](KeyType&& key) const {return (key[0]==key[1] ? s : S());}
|
||||
template<typename A1, typename A2> 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<typename T, bool IsMatrix> struct node_t :
|
||||
std::conditional<utility::is_in_ZRC<T>::value, _scalar_wrap<typename std::remove_reference<T>::type,IsMatrix>, typename utility::remove_rvalue_ref<T>::type> {};
|
||||
//std::conditional<utility::is_in_ZRC<T>::value, _scalar_wrap<typename std::add_const<T>::type,IsMatrix>, typename utility::remove_rvalue_ref<T>::type> {};
|
||||
|
||||
// get the rank of something ....
|
||||
template<typename T> struct get_rank { static constexpr int value = std::remove_reference<T>::type::domain_type::rank;};
|
||||
template<typename S, bool IsMatrix> struct get_rank<scalar_wrap<S,IsMatrix>> { static constexpr int value =0;};
|
||||
template<typename S, bool IsMatrix> struct get_rank<_scalar_wrap<S,IsMatrix>> { static constexpr int value =0;};
|
||||
|
||||
//
|
||||
template<typename T, bool IsMatrix=false> struct keeper_type : boost::mpl::if_<is_in_ZRC<T>, scalar_wrap<T,IsMatrix>, typename view_type_if_exists_else_type<T>::type> {};
|
||||
//template<typename T, bool IsMatrix=false> struct keeper_type : boost::mpl::if_<is_in_ZRC<T>, _scalar_wrap<T,IsMatrix>, typename view_type_if_exists_else_type<T>::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" <<r.domain().lengths();
|
||||
return l.domain();
|
||||
}
|
||||
template<typename S, bool IsMatrix, typename R> auto operator() (scalar_wrap<S,IsMatrix> const & w, R const & r) const -> decltype(r.domain()) { return r.domain();}
|
||||
template<typename S, bool IsMatrix, typename L> auto operator() (L const & l, scalar_wrap<S,IsMatrix> const & w) const -> decltype(l.domain()) { return l.domain();}
|
||||
template<typename S, bool IsMatrix, typename R> auto operator() (_scalar_wrap<S,IsMatrix> const & w, R const & r) const -> decltype(r.domain()) { return r.domain();}
|
||||
template<typename S, bool IsMatrix, typename L> auto operator() (L const & l, _scalar_wrap<S,IsMatrix> const & w) const -> decltype(l.domain()) { return l.domain();}
|
||||
};
|
||||
|
||||
template<typename A, typename B>
|
||||
struct type_of_mult{
|
||||
//typedef typename boost::remove_reference<A>::type T1; typedef typename boost::remove_reference<B>::type T2;
|
||||
//typedef BOOST_TYPEOF_TPL( (*(pseudo_default_construct<T1>())) * (*(pseudo_default_construct<T2>()))) type;
|
||||
typedef decltype ( std::declval<typename std::remove_reference<A>::type>() * std::declval<typename std::remove_reference<B>::type>() ) type;
|
||||
};
|
||||
|
||||
|
||||
}}//namespace triqs::arrays
|
||||
#endif
|
||||
|
||||
@ -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<typename Tag, typename L, typename R>
|
||||
struct vector_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableVector) {
|
||||
typedef typename keeper_type<L>::type L_t;
|
||||
typedef typename keeper_type<R>::type R_t;
|
||||
L_t l; R_t r;
|
||||
template<typename LL, typename RR> vector_expr(LL && l_, RR && r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
|
||||
typedef typename std::remove_reference<L>::type L_t;
|
||||
typedef typename std::remove_reference<R>::type R_t;
|
||||
static_assert( get_rank<R_t>::value==0 || get_rank<L_t>::value==0 || get_rank<L_t>::value == get_rank<R_t>::value, "rank mismatch in array operations");
|
||||
typedef typename std::result_of<operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::result_of<utility::operation<Tag>(typename L_t::value_type,typename R_t::value_type)>::type value_type;
|
||||
typedef typename std::remove_reference<typename std::result_of<combine_domain(L_t,R_t)>::type>::type domain_type;
|
||||
|
||||
L l; R r;
|
||||
template<typename LL, typename RR> vector_expr(LL && l_, RR && r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
|
||||
|
||||
domain_type domain() const { return combine_domain()(l,r); }
|
||||
mini_vector<size_t,1> shape() const { return this->domain().lengths();}
|
||||
//mini_vector<size_t,1> shape() const { return this->domain().lengths();}
|
||||
size_t size() const { return this->domain().lengths()[0];}
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType && key) const { return operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]);}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
friend std::ostream &operator <<(std::ostream &sout, vector_expr const &expr){return sout << "("<<expr.l << " "<<operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
template<typename KeyType> value_type operator[](KeyType && key) const { return utility::operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]);}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...));}
|
||||
friend std::ostream &operator <<(std::ostream &sout, vector_expr const &expr){return sout << "("<<expr.l << " "<<utility::operation<Tag>::name << " "<<expr.r<<")" ; }
|
||||
friend arrays::vector<value_type> make_vector(vector_expr const & e) { return e;}
|
||||
};
|
||||
|
||||
template<typename L> // a special case : the unary operator !
|
||||
struct vector_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableVector) {
|
||||
typedef typename keeper_type<L>::type L_t;
|
||||
L_t l;
|
||||
template<typename LL> vector_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
|
||||
typedef typename std::remove_reference<L>::type L_t;
|
||||
typedef typename L_t::value_type value_type;
|
||||
typedef typename L_t::domain_type domain_type;
|
||||
|
||||
L l;
|
||||
template<typename LL> vector_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
|
||||
|
||||
domain_type domain() const { return l.domain(); }
|
||||
mini_vector<size_t,1> shape() const { return this->domain().lengths();}
|
||||
//mini_vector<size_t,1> shape() const { return this->domain().lengths();}
|
||||
size_t size() const { return this->domain().lengths()[0];}
|
||||
|
||||
//template<typename KeyType> value_type operator[](KeyType&& key) const {return -l[key];}
|
||||
template<typename KeyType> value_type operator[](KeyType&& key) const {return -l[std::forward<KeyType>(key)];}
|
||||
template<typename ... Args> value_type operator()(Args && ... args) const { return -l(std::forward<Args>(args)...);}
|
||||
friend std::ostream &operator <<(std::ostream &sout, vector_unary_m_expr const &expr){return sout << '-'<<expr.l; }
|
||||
friend arrays::vector<value_type> 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 A1, typename A2>\
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, vector_expr<tags::TAG, A1,A2>>::type\
|
||||
operator OP (A1 const & a1, A2 const & a2) { return vector_expr<tags::TAG, A1,A2>(a1,a2);}
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value,\
|
||||
vector_expr<utility::tags::TAG, typename node_t<A1,false>::type, typename node_t<A2,false>::type>>::type\
|
||||
operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(a2)};}
|
||||
#define DELETE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \
|
||||
template<typename A1, typename A2>\
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, vector_expr<tags::TAG, A1,A2>>::type\
|
||||
operator OP (A1 const & a1, A2 const & a2) = delete;
|
||||
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, vector_expr<utility::tags::TAG, A1,A2>>::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 A1> typename std::enable_if<ImmutableVector<A1>::value, vector_unary_m_expr<A1>>::type
|
||||
operator - (A1 const & a1) { return {a1};}
|
||||
|
||||
template<typename Expr > vector<typename Expr::value_type>
|
||||
make_vector( Expr const & e) { return vector<typename Expr::value_type>(e);}
|
||||
template<typename A1>
|
||||
typename std::enable_if<
|
||||
ImmutableVector<A1>::value,
|
||||
vector_unary_m_expr<typename node_t<A1,true>::type >
|
||||
>::type
|
||||
operator - (A1 && a1) { return {std::forward<A1>(a1)};}
|
||||
|
||||
}}//namespace triqs::arrays
|
||||
#endif
|
||||
|
||||
@ -116,11 +116,12 @@ namespace triqs { namespace arrays {
|
||||
mini_vector<hsize_t,Rank> dims_out;
|
||||
//int ndims = dataspace.getSimpleExtentDims( &dims_out[0], NULL);
|
||||
dataspace.getSimpleExtentDims( &dims_out[0], NULL);
|
||||
mini_vector<size_t,ArrayType::rank > d2; for (size_t u=0; u<ArrayType::rank ; ++u) d2[u] = dims_out[u];
|
||||
resize_or_check(A, d2 );
|
||||
mini_vector<size_t,ArrayType::rank > d2;
|
||||
for (size_t u=0; u<ArrayType::rank ; ++u) d2[u] = dims_out[u];
|
||||
resize_or_check(A, d2);
|
||||
if (C_reorder) { read_array(g,name, cache<ArrayType,typename ArrayType::regular_type>(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<typename ArrayType::value_type>(), data_space(A) , dataspace );
|
||||
else ds.read( get_array_data_ptr(A), h5::data_type_memory<typename ArrayType::value_type>(), data_space(A) , dataspace );
|
||||
}
|
||||
TRIQS_ARRAYS_H5_CATCH_EXCEPTION;
|
||||
}
|
||||
|
||||
@ -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<OptionFlags>::type() );
|
||||
|
||||
@ -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<typename A, class Enable = void> struct inverse_lazy;
|
||||
template<typename A> struct inverse_lazy_impl;//debug ony
|
||||
|
||||
///
|
||||
template<typename A> struct determinant_lazy;
|
||||
|
||||
/// Lazy inversion
|
||||
template<class A> inverse_lazy<A> inverse (A const & a) { return inverse_lazy<A>(a); }
|
||||
template<class A> inverse_lazy<typename utility::remove_rvalue_ref<A>::type> inverse (A && a) { return {std::forward<A>(a)}; }
|
||||
|
||||
/// Lazy computation of det
|
||||
template<typename A> determinant_lazy<A> determinant (A const & a) { return determinant_lazy<A>(a); }
|
||||
@ -111,58 +112,67 @@ namespace triqs { namespace arrays {
|
||||
|
||||
// an implementation class to gather the common part to matrix and expression....
|
||||
template<typename A> struct inverse_lazy_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) {
|
||||
typedef typename std::remove_const<typename A::value_type>::type value_type;
|
||||
typedef typename A::domain_type domain_type;
|
||||
typedef typename const_view_type_if_exists_else_type<A>::type A_type;
|
||||
const A_type a;
|
||||
inverse_lazy_impl(A const & a_):a (a_) {
|
||||
public:
|
||||
|
||||
typedef typename std::remove_reference<A>::type A_t;
|
||||
typedef typename std::remove_const<typename A_t::value_type>::type value_type;
|
||||
typedef typename A_t::domain_type domain_type;
|
||||
A a;
|
||||
template<typename AA> inverse_lazy_impl(AA && a_):a (std::forward<AA>(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<typename K0, typename K1> value_type operator() (K0 const & k0, K1 const & k1) const { activate(); return _id->M(k0,k1); }
|
||||
template<typename K0, typename K1> 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("<<x.a<<")";}
|
||||
protected:
|
||||
struct internal_data { // implementing the pattern LazyPreCompute
|
||||
//typedef typename A_type::regular_type M_type;
|
||||
typedef matrix<value_type> M_type;
|
||||
typedef matrix_view<value_type> M_view_type;
|
||||
M_type M;
|
||||
internal_data(inverse_lazy_impl const & P):M(P.a){det_and_inverse_worker<M_view_type> worker(M); worker.inverse();}
|
||||
internal_data(inverse_lazy_impl const & P):M(P.a){
|
||||
det_and_inverse_worker<M_view_type> worker(M);
|
||||
worker.inverse();
|
||||
}
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
};
|
||||
|
||||
// The general case
|
||||
template<typename A>
|
||||
struct inverse_lazy<A,typename boost::disable_if< is_matrix_or_view<A> >::type > : inverse_lazy_impl<A> {
|
||||
inverse_lazy(A const & a):inverse_lazy_impl<A>(a) { }
|
||||
struct inverse_lazy<A, DISABLE_IF(is_matrix_or_view<typename std::remove_reference<A>::type>) > : inverse_lazy_impl<A> {
|
||||
template<typename AA> inverse_lazy(AA && a_):inverse_lazy_impl<A>(std::forward<AA>(a_)) {}
|
||||
};
|
||||
|
||||
// for matrix and matrix_views, we have more optimisation possible ....
|
||||
template<typename A>
|
||||
struct inverse_lazy<A,typename boost::enable_if< is_matrix_or_view<A> >::type >:inverse_lazy_impl<A>{
|
||||
inverse_lazy(A const & a):inverse_lazy_impl<A>(a) { }
|
||||
struct inverse_lazy<A, ENABLE_IF(is_matrix_or_view<typename std::remove_reference<A>::type>) >:inverse_lazy_impl<A>{
|
||||
template<typename AA> inverse_lazy(AA && a_):inverse_lazy_impl<A>(std::forward<AA>(a_)) {}
|
||||
|
||||
template<typename MT> // Optimized implementation of =
|
||||
friend void triqs_arrays_assign_delegation (MT & lhs, inverse_lazy const & rhs) {
|
||||
//std::cerr << " DELEGATING1"<< lhs <<std::endl;
|
||||
std::cerr << " DELEGATING2"<< rhs << rhs.a<<std::endl;
|
||||
static_assert(is_matrix_or_view<MT>::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<MT> C(lhs);// a reflexive cache will use a temporary "regrouping" copy if and only if needed
|
||||
det_and_inverse_worker<typename MT::view_type> 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("<<x.a<<")";}
|
||||
};
|
||||
|
||||
//------------------- det ----------------------------------------
|
||||
|
||||
template<typename A> struct determinant_lazy { // : { Tag::expression_terminal, Tag::scalar_expression_terminal {
|
||||
template<typename A> struct determinant_lazy {
|
||||
typedef typename A::value_type value_type;
|
||||
typedef typename const_view_type_if_exists_else_type<A>::type A_type;
|
||||
A_type a;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user