mirror of
https://github.com/triqs/dft_tools
synced 2024-12-25 13:53:40 +01:00
arrays: new ref counting and weak views
- improve the mem_block and shared_block. - the reference counting is now done in the mem_block and shared_block, removing the need of shared_ptr. - speed tests shows that shared_ptr is very slow (due to thread safety?) the new version is much better, though not perfect. - Hence introducing weak views. - also : -- clean the guard mechanism for python (to allow returning from python without any python ref left). -- clean code, add documentation for mem_block -- remove nan init, which was not working, and corresponding test -- serialisation of view still unchanged (need to forbid serialization of view ??). - tests ok, incl. valgrind tests.
This commit is contained in:
Olivier Parcollet
parent
6bf0d0f6b5
commit
fc2a620eae
77
test/speed/array_view.cpp
Normal file
77
test/speed/array_view.cpp
Normal file
@ -0,0 +1,77 @@
|
||||
//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
|
||||
#include <triqs/arrays.hpp>
|
||||
#include <triqs/arrays/matrix_view_proxy.hpp>
|
||||
|
||||
using namespace triqs::arrays;
|
||||
|
||||
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
|
||||
|
||||
const int nl_interne = 100000;
|
||||
const int N = 1000;
|
||||
|
||||
typedef double VALUE_TYPE;
|
||||
//typedef int VALUE_TYPE;
|
||||
inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8.0);}
|
||||
//inline VALUE_TYPE fnt(size_t i) { return i;} //*(i+2.0)*(i-8);}
|
||||
//inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8);}
|
||||
|
||||
struct with_view {
|
||||
void operator()() {
|
||||
|
||||
array<double,3> A(N,2,2);
|
||||
A() =0;
|
||||
|
||||
for (int u =0; u<nl_interne; ++u)
|
||||
for (int i =0; i<N-1; ++i) A(i,range(), range()) (0,0) = fnt(i);
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
struct with_proxy {
|
||||
void operator()() {
|
||||
|
||||
array<double,3> A(N,2,2);
|
||||
A() =0;
|
||||
|
||||
for (int u =0; u<nl_interne; ++u)
|
||||
for (int i =0; i<N-1; ++i) matrix_view_proxy<array<double,3>, 0>(A,i) (0,0) = fnt(i);
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
struct direct {
|
||||
void operator()() {
|
||||
|
||||
array<double,3> A(N,2,2);
|
||||
A() =0;
|
||||
|
||||
for (int u =0; u<nl_interne; ++u)
|
||||
for (int i =0; i<N-1; ++i) A(i,0,0) = fnt(i);
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
struct direct2 {
|
||||
void operator()() {
|
||||
|
||||
array<double,1> A(N);
|
||||
A() =0;
|
||||
|
||||
for (int u =0; u<nl_interne; ++u)
|
||||
for (int i =0; i<N-1; ++i) A(i) = fnt(i);
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
#include "./speed_tester.hpp"
|
||||
int main() {
|
||||
speed_tester<direct> (50);
|
||||
speed_tester<with_proxy> (50);
|
||||
speed_tester<with_view> (50);
|
||||
speed_tester<direct2> (50);
|
||||
//speed_tester<with_slices> (5000);
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -22,9 +22,12 @@
|
||||
#ifndef COMMON_TEST_ARRAY_H
|
||||
#define COMMON_TEST_ARRAY_H
|
||||
#include<triqs/utility/first_include.hpp>
|
||||
#include<sstream>
|
||||
|
||||
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl;
|
||||
|
||||
#define TEST_ASSERT(X) if(!bool(X)) { std::stringstream fs; fs<< "Failed Assertion line "<< __LINE__<< " of file "<< __FILE__<<" :\n"<<BOOST_PP_STRINGIZE((X)); std::cout << fs.str()<< std::endl ; TRIQS_RUNTIME_ERROR<< fs.str();}
|
||||
|
||||
#define AS_STRING(X) AS_STRING2(X)
|
||||
#define AS_STRING2(X) #X
|
||||
|
||||
|
||||
@ -1,43 +0,0 @@
|
||||
/*******************************************************************************
|
||||
*
|
||||
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||
*
|
||||
* Copyright (C) 2011 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
|
||||
* Foundation, either version 3 of the License, or (at your option) any later
|
||||
* version.
|
||||
*
|
||||
* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License along with
|
||||
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
******************************************************************************/
|
||||
#include "./common.hpp"
|
||||
#include "./src/array.hpp"
|
||||
#include "./src/matrix.hpp"
|
||||
#include <iostream>
|
||||
|
||||
using std::cout; using std::endl;
|
||||
using namespace triqs::arrays;
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
|
||||
|
||||
|
||||
array<long,2, NanInit > A (2,2);
|
||||
array<double,2, NanInit > B (2,2);
|
||||
array<std::complex<double> ,2, NanInit > C (2,2);
|
||||
|
||||
std::cerr << A << std::endl ;
|
||||
std::cerr << B << std::endl ;
|
||||
std::cerr << C << std::endl ;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
@ -51,16 +51,14 @@ int main(int argc, char **argv) {
|
||||
triqs::arrays::vector<double> V(3), W(4);;
|
||||
V() = 3; W()=4;
|
||||
|
||||
auto VV = V(tqa::range (0,2));
|
||||
auto VW = W(tqa::range (0,2));
|
||||
|
||||
std::cout << "V = "<< V << " W = " << W<< " V view "<< VV<< " W view "<< VW<< std::endl;
|
||||
std::cout << "V = "<< V << " W = " << W<< std::endl;
|
||||
std::swap(V,W);
|
||||
std::cout << "V = "<< V << " W = " << W<< " V view "<< VV<< " W view "<< VW<< std::endl;
|
||||
std::cout << "V = "<< V << " W = " << W<< std::endl;
|
||||
|
||||
// This does NOT compile, the std::swap specialization is deleted.
|
||||
auto VV = V(tqa::range (0,2));
|
||||
auto VW = W(tqa::range (0,2));
|
||||
//std::swap(VV,VW);
|
||||
//std::cout << "V = "<< V << " W = " << W<< " V view "<< VV<< " W view "<< VW<< std::endl;
|
||||
|
||||
}
|
||||
std::cout << "With deep_swap"<< std::endl;
|
||||
|
||||
@ -3,8 +3,8 @@ V = [3,3,3] W = [4,4,4,4] V view [3,3] W view [4,4]
|
||||
V = [4,4,4,4] W = [3,3,3] V view [3,3] W view [4,4]
|
||||
V = [4,4,4,4] W = [3,3,3] V view [4,4] W view [3,3]
|
||||
With std::swap
|
||||
V = [3,3,3] W = [4,4,4,4] V view [3,3] W view [4,4]
|
||||
V = [4,4,4,4] W = [3,3,3] V view [3,3] W view [4,4]
|
||||
V = [3,3,3] W = [4,4,4,4]
|
||||
V = [4,4,4,4] W = [3,3,3]
|
||||
With deep_swap
|
||||
V = [3,3,3] W = [4,4,4] V view [3,3] W view [4,4]
|
||||
V = [4,4,4] W = [3,3,3] V view [4,4] W view [3,3]
|
||||
|
||||
@ -48,7 +48,7 @@ int main(int argc, char **argv) {
|
||||
for (int j=0; j<3; ++j)
|
||||
A(i,j) = 10*i+ j;
|
||||
|
||||
array_view<long const ,2 > AA (A);
|
||||
array_view<long,2 > AA (A);
|
||||
|
||||
std::cout<<"A is "<<A<<std::endl;
|
||||
std::cout<<"another view : "<<AA<<std::endl;
|
||||
|
||||
@ -28,22 +28,23 @@
|
||||
#include "impl/flags.hpp"
|
||||
namespace triqs { namespace arrays {
|
||||
|
||||
template <typename ValueType, int Rank, ull_t Opt=0, ull_t TraversalOrder= 0> class array_view;
|
||||
template <typename ValueType, int Rank, ull_t Opt=0, ull_t TraversalOrder= 0, bool Borrowed=false > class array_view;
|
||||
template <typename ValueType, int Rank, ull_t Opt=0, ull_t TraversalOrder= 0> class array;
|
||||
|
||||
// ---------------------- array_view --------------------------------
|
||||
|
||||
#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<Rank,Opt,TraversalOrder>, \
|
||||
storages::shared_block<ValueType>, Opt, TraversalOrder, Tag::array_view >
|
||||
storages::shared_block<ValueType, Borrowed>, Opt, TraversalOrder, Tag::array_view >
|
||||
|
||||
template <typename ValueType, int Rank, ull_t Opt, ull_t TraversalOrder>
|
||||
template <typename ValueType, int Rank, ull_t Opt, ull_t TraversalOrder, bool Borrowed>
|
||||
class array_view : Tag::array_view, TRIQS_MODEL_CONCEPT(MutableCuboidArray), public IMPL_TYPE {
|
||||
static_assert( Rank>0, " Rank must be >0");
|
||||
public:
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
typedef typename IMPL_TYPE::storage_type storage_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder> view_type;
|
||||
typedef array<ValueType,Rank,Opt,TraversalOrder> non_view_type;
|
||||
typedef array <ValueType,Rank,Opt,TraversalOrder> non_view_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder> view_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
/// Build from an IndexMap and a storage
|
||||
@ -80,19 +81,27 @@ namespace triqs { namespace arrays {
|
||||
// Move assignment not defined : will use the copy = since view must copy data
|
||||
|
||||
TRIQS_DEFINE_COMPOUND_OPERATORS(array_view);
|
||||
// to forbid serialization of views...
|
||||
//template<class Archive> void serialize(Archive & ar, const unsigned int version) = delete;
|
||||
|
||||
};
|
||||
|
||||
#undef IMPL_TYPE
|
||||
//------------------------------- array ---------------------------------------------------
|
||||
|
||||
#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<Rank,Opt,TraversalOrder>, \
|
||||
storages::shared_block<ValueType>, Opt, TraversalOrder, Tag::array_view >
|
||||
|
||||
template <typename ValueType, int Rank, ull_t Opt, ull_t TraversalOrder>
|
||||
class array: Tag::array, TRIQS_MODEL_CONCEPT(MutableCuboidArray), public IMPL_TYPE {
|
||||
public:
|
||||
typedef typename IMPL_TYPE::value_type value_type;
|
||||
typedef typename IMPL_TYPE::storage_type storage_type;
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder> view_type;
|
||||
typedef array<ValueType,Rank,Opt,TraversalOrder> non_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
typedef array <ValueType,Rank,Opt,TraversalOrder> non_view_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder> view_type;
|
||||
typedef array_view<ValueType,Rank,Opt,TraversalOrder,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
/// Empty array.
|
||||
explicit array(memory_layout<Rank> ml = memory_layout<Rank>(IMPL_TYPE::indexmap_type::traversal_order)) :IMPL_TYPE(ml){}
|
||||
@ -190,7 +199,8 @@ namespace triqs { namespace arrays {
|
||||
//----------------------------------------------------------------------------------
|
||||
|
||||
// how to build the view type ....
|
||||
template < class V, int R, ull_t Opt, ull_t TraversalOrder > struct ViewFactory< V, R, Opt, TraversalOrder,Tag::array_view > { typedef array_view<V,R,Opt,TraversalOrder> type; };
|
||||
template < class V, int R, ull_t Opt, ull_t TraversalOrder, bool Borrowed >
|
||||
struct ISPViewType< V, R, Opt, TraversalOrder,Tag::array_view, Borrowed> { typedef array_view<V,R,Opt,TraversalOrder, Borrowed> type; };
|
||||
|
||||
}}//namespace triqs::arrays
|
||||
|
||||
|
||||
@ -79,10 +79,10 @@ namespace triqs { namespace arrays { namespace blas_lapack_tools {
|
||||
// first the function to take the decision
|
||||
// should be copy to call blas/lapack ? only if we have a view and if the min_stride
|
||||
// of the matrix is not 1, otherwise we use the LD parameter of blas/lapack
|
||||
template <typename T> constexpr bool copy_is_needed (T const & A) { return false;}
|
||||
template <typename T> constexpr TYPE_ENABLE_IFC(bool, !is_matrix_view<T>::value) copy_is_needed (T const & A) { return false;}
|
||||
|
||||
template <typename T, ull_t Opt, ull_t To>
|
||||
bool copy_is_needed (matrix_view<T,Opt,To> const & A) {
|
||||
template <typename MatrixView>
|
||||
TYPE_ENABLE_IF(bool, is_matrix_view<MatrixView>) copy_is_needed (MatrixView const & A) {
|
||||
auto min_stride = A.indexmap().strides()[A.memory_layout_is_fortran() ? 0 : 1];
|
||||
return min_stride !=1;
|
||||
}
|
||||
|
||||
@ -65,8 +65,8 @@ namespace triqs { namespace arrays {
|
||||
}
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable boost::shared_ptr<internal_data> _id;
|
||||
internal_data & id() const { if (!_id) _id= boost::make_shared<internal_data>(*this,ml); return *_id;}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
internal_data & id() const { if (!_id) _id= std::make_shared<internal_data>(*this,ml); return *_id;}
|
||||
|
||||
//exposed_view_type & view2 () { if (need_copy) return id().view; else return keeper;}
|
||||
//exposed_view_type const & view2 () const { if (need_copy) return id().view; else return keeper;}
|
||||
|
||||
@ -33,8 +33,6 @@
|
||||
#include <sstream>
|
||||
|
||||
#include <boost/utility/enable_if.hpp>
|
||||
#include <boost/shared_ptr.hpp>
|
||||
#include <boost/make_shared.hpp>
|
||||
#include <boost/mpl/or.hpp>
|
||||
#include <boost/mpl/and.hpp>
|
||||
#include <boost/mpl/not.hpp>
|
||||
@ -57,12 +55,17 @@ namespace boost { namespace serialization { class access;}}
|
||||
|
||||
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
|
||||
#define TRIQS_ARRAYS_CHECK_IM_STORAGE_COMPAT
|
||||
#define TRIQS_ARRAYS_ENFORCE_INIT_NAN_INF
|
||||
|
||||
#else
|
||||
#define TRIQS_ARRAYS_DEBUG_CHECK(Cond,Error)
|
||||
#endif
|
||||
|
||||
// In particular, gcc 4.6, 4.7 do not have "rvalue for *this", which forbid the implementation
|
||||
// of weak references. In that we, we will revert to strong view instead.
|
||||
#ifdef TRIQS_COMPILER_IS_OBSOLETE
|
||||
#define TRIQS_ARRAYS_SLICE_DEFAUT_IS_SHARED
|
||||
#endif
|
||||
|
||||
namespace triqs {
|
||||
|
||||
typedef unsigned long long ull_t;
|
||||
|
||||
@ -37,7 +37,7 @@ namespace triqs { namespace arrays {
|
||||
constexpr ull_t BoundCheck = 1ull << 0;
|
||||
constexpr ull_t TraversalOrderC = 1ull << 1;
|
||||
constexpr ull_t TraversalOrderFortran = 1ull << 2;
|
||||
constexpr ull_t NanInit = 1ull << 3;
|
||||
//constexpr ull_t NanInit = 1ull << 3;
|
||||
constexpr ull_t DefaultInit = 1ull << 4;
|
||||
|
||||
#define BOUND_CHECK triqs::arrays::BoundCheck
|
||||
@ -48,12 +48,12 @@ namespace triqs { namespace arrays {
|
||||
// NB : flags MUST be insensitive to slicing ...
|
||||
// i.e. when I slice, the flags does not change.
|
||||
|
||||
namespace flags {
|
||||
namespace flags {
|
||||
|
||||
#ifndef TRIQS_WORKAROUND_INTEL_COMPILER_BUGS
|
||||
|
||||
constexpr ull_t get(ull_t f, ull_t a) { return ((f & (1ull<<a)) >> a);}
|
||||
constexpr ull_t get2(ull_t f, ull_t a) { return (f & (((1ull<<a) + (1ull<< (a+1ull))) >> a) );}
|
||||
constexpr ull_t get2(ull_t f, ull_t a) { return ((f & ((1ull<<a) + (1ull<< (a+1ull)))) >> a);}
|
||||
|
||||
#ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
|
||||
constexpr bool bound_check (ull_t f) { return true;}
|
||||
@ -71,14 +71,14 @@ namespace triqs { namespace arrays {
|
||||
|
||||
template<ull_t F> struct init_tag1;
|
||||
template<> struct init_tag1<0> { typedef Tag::no_init type;};
|
||||
template<> struct init_tag1<1> { typedef Tag::nan_inf_init type;};
|
||||
//template<> struct init_tag1<1> { typedef Tag::nan_inf_init type;};
|
||||
template<> struct init_tag1<2> { typedef Tag::default_init type;};
|
||||
|
||||
// for the init_tag, we pass the *whole* option flag.
|
||||
template<ull_t F> struct init_tag : init_tag1 < init_mode(F)> {};
|
||||
|
||||
template<ull_t F, ull_t To> struct assert_make_sense {
|
||||
static constexpr bool bc = bound_check(F);;
|
||||
static constexpr bool bc = bound_check(F);
|
||||
static constexpr bool is_c = traversal_order_c(F);
|
||||
static constexpr bool is_f = traversal_order_fortran(F);
|
||||
static constexpr bool inm = check_init_mode(F);
|
||||
@ -90,8 +90,8 @@ namespace triqs { namespace arrays {
|
||||
#else
|
||||
|
||||
#define TRIQS_FLAGS_GET(f,a) ((f & (1ull<<a)) >> a)
|
||||
#define TRIQS_FLAGS_GET2(f,a) (f & (((1ull<<a) + (1ull<< (a+1ull))) >> a))
|
||||
|
||||
#define TRIQS_FLAGS_GET2(f,a) ((f & ((1ull<<a) + (1ull<< (a+1ull)))) >> a)
|
||||
|
||||
#ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
|
||||
template< ull_t F> struct bound_check_trait { static constexpr bool value = true;};
|
||||
#else
|
||||
@ -105,7 +105,7 @@ namespace triqs { namespace arrays {
|
||||
|
||||
template<ull_t F> struct init_tag1;
|
||||
template<> struct init_tag1<0> { typedef Tag::no_init type;};
|
||||
template<> struct init_tag1<1> { typedef Tag::nan_inf_init type;};
|
||||
//template<> struct init_tag1<1> { typedef Tag::nan_inf_init type;};
|
||||
template<> struct init_tag1<2> { typedef Tag::default_init type;};
|
||||
template<ull_t F> struct init_tag : init_tag1 < init_mode_tr<F>::value > {};
|
||||
|
||||
|
||||
@ -40,7 +40,7 @@ namespace triqs { namespace arrays {
|
||||
|
||||
template <bool Const, typename IndexMapIterator, typename StorageType > class iterator_adapter;
|
||||
|
||||
template < class V, int R, ull_t OptionFlags, ull_t TraversalOrder, class ViewTag > struct ViewFactory;
|
||||
template < class V, int R, ull_t OptionFlags, ull_t TraversalOrder, class ViewTag, bool Borrowed > struct ISPViewType;
|
||||
|
||||
template <typename IndexMapType, typename StorageType, ull_t OptionFlags, ull_t TraversalOrder, typename ViewTag >
|
||||
class indexmap_storage_pair : Tag::indexmap_storage_pair, TRIQS_MODEL_CONCEPT(MutableArray) {
|
||||
@ -48,7 +48,7 @@ namespace triqs { namespace arrays {
|
||||
public :
|
||||
typedef typename StorageType::value_type value_type;
|
||||
static_assert(std::is_constructible<value_type>::value, "array/array_view and const operate only on values");
|
||||
//static_assert(!std::is_const<value_type>::value, "no const type");
|
||||
static_assert(!std::is_const<value_type>::value, "no const type");
|
||||
typedef StorageType storage_type;
|
||||
typedef IndexMapType indexmap_type;
|
||||
static constexpr unsigned int rank = IndexMapType::domain_type::rank;
|
||||
@ -63,6 +63,7 @@ namespace triqs { namespace arrays {
|
||||
|
||||
indexmap_storage_pair (const indexmap_type & IM, const storage_type & ST):
|
||||
indexmap_(IM),storage_(ST){
|
||||
//std::cerr << " construct ISP "<< storage_type::is_weak<< std::endl;
|
||||
#ifdef TRIQS_ARRAYS_CHECK_IM_STORAGE_COMPAT
|
||||
if (ST.size() != IM.domain().number_of_elements())
|
||||
TRIQS_RUNTIME_ERROR<<"index_storage_pair construction : storage and indices are not compatible";
|
||||
@ -71,12 +72,21 @@ namespace triqs { namespace arrays {
|
||||
|
||||
indexmap_storage_pair (indexmap_type && IM, storage_type && ST):
|
||||
indexmap_(std::move(IM)),storage_(std::move(ST)){
|
||||
//std::cerr << " construct ISP && IM, && ST "<< storage_type::is_weak<< std::endl;
|
||||
#ifdef TRIQS_ARRAYS_CHECK_IM_STORAGE_COMPAT
|
||||
if (ST.size() != IM.domain().number_of_elements())
|
||||
TRIQS_RUNTIME_ERROR<<"index_storage_pair construction : storage and indices are not compatible";
|
||||
#endif
|
||||
}
|
||||
|
||||
indexmap_storage_pair (indexmap_type const & IM, storage_type && ST):
|
||||
indexmap_(IM),storage_(std::move(ST)){
|
||||
//std::cerr << " construct ISP && ST "<< storage_type::is_weak<< std::endl;
|
||||
#ifdef TRIQS_ARRAYS_CHECK_IM_STORAGE_COMPAT
|
||||
if (ST.size() != IM.domain().number_of_elements())
|
||||
TRIQS_RUNTIME_ERROR<<"index_storage_pair construction : storage and indices are not compatible";
|
||||
#endif
|
||||
}
|
||||
|
||||
/// The storage is allocated from the size of IM.
|
||||
indexmap_storage_pair (const indexmap_type & IM): indexmap_(IM),storage_(){
|
||||
@ -107,7 +117,7 @@ namespace triqs { namespace arrays {
|
||||
<<"\nfrom the python object \n"<< numpy_interface::object_to_string(X)
|
||||
<<"\nThe error was :\n "<<s.what();
|
||||
}
|
||||
//std::cout << " Leave IPS ref count = "<< X->ob_refcnt << std::endl;
|
||||
//std::cerr << " Leave IPS ref count = "<< X->ob_refcnt << std::endl;
|
||||
}
|
||||
#endif
|
||||
// ------------------------------- swap --------------------------------------------
|
||||
@ -172,6 +182,9 @@ namespace triqs { namespace arrays {
|
||||
|
||||
// ------------------------------- operator () --------------------------------------------
|
||||
|
||||
typedef typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag,false >::type view_type;
|
||||
typedef typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag,true >::type weak_view_type;
|
||||
|
||||
// Evaluation and slices
|
||||
template<typename... Args>
|
||||
typename std::enable_if<
|
||||
@ -185,32 +198,82 @@ namespace triqs { namespace arrays {
|
||||
, value_type const &>::type
|
||||
operator()(Args const & ... args) const { return storage_[indexmap_(args...)]; }
|
||||
|
||||
template<bool is_const, typename ... Args> struct result_of_call_as_view {
|
||||
template<bool is_const, bool ForceBorrowed, typename ... Args> struct result_of_call_as_view {
|
||||
typedef typename indexmaps::slicer<indexmap_type,Args...>::r_type IM2;
|
||||
typedef typename std::conditional<is_const, typename std::add_const<value_type>::type, value_type>::type V2;
|
||||
//typedef typename std::conditional<is_const, typename std::add_const<value_type>::type, value_type>::type V2;
|
||||
typedef value_type V2;
|
||||
static_assert(IM2::domain_type::rank !=0, "Internal error");
|
||||
typedef typename ViewFactory<V2,IM2::domain_type::rank, OptionFlags, IM2::traversal_order_in_template, ViewTag >::type type;
|
||||
typedef typename ISPViewType<V2,IM2::domain_type::rank, OptionFlags, IM2::traversal_order_in_template, ViewTag, ForceBorrowed || StorageType::is_weak >::type type;
|
||||
};
|
||||
|
||||
// change here only for debug
|
||||
//#define TRIQS_ARRAYS_SLICE_DEFAUT_IS_SHARED
|
||||
#ifndef TRIQS_ARRAYS_SLICE_DEFAUT_IS_SHARED
|
||||
template<typename... Args> // non const version
|
||||
typename boost::lazy_enable_if_c<
|
||||
(!clef::is_any_lazy<Args...>::value) && (indexmaps::slicer<indexmap_type,Args...>::r_type::domain_type::rank!=0)
|
||||
, result_of_call_as_view<false,Args...>
|
||||
, result_of_call_as_view<false,true,Args...>
|
||||
>::type // enable_if
|
||||
operator()(Args const & ... args) {
|
||||
return typename result_of_call_as_view<false,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
operator()(Args const & ... args) & {
|
||||
return typename result_of_call_as_view<false,true,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
|
||||
template<typename... Args> // const version
|
||||
typename boost::lazy_enable_if_c<
|
||||
(!clef::is_any_lazy<Args...>::value) && (indexmaps::slicer<indexmap_type,Args...>::r_type::domain_type::rank!=0)
|
||||
, result_of_call_as_view<true,Args...>
|
||||
, result_of_call_as_view<true,true,Args...>
|
||||
>::type // enable_if
|
||||
operator()(Args const & ... args) const & {
|
||||
return typename result_of_call_as_view<true,true,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
|
||||
template<typename... Args> // rvalue version : same value of weak as this
|
||||
typename boost::lazy_enable_if_c<
|
||||
(!clef::is_any_lazy<Args...>::value) && (indexmaps::slicer<indexmap_type,Args...>::r_type::domain_type::rank!=0)
|
||||
, result_of_call_as_view<true,false,Args...>
|
||||
>::type // enable_if
|
||||
operator()(Args const & ... args) && {
|
||||
//std::cerr << "slicing a temporary"<< this->storage().is_weak<< result_of_call_as_view<true,true,Args...>::type::storage_type::is_weak << std::endl;
|
||||
return typename result_of_call_as_view<true,false,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), std::move(storage()));
|
||||
}
|
||||
|
||||
/// Equivalent to make_view
|
||||
//typename ISPViewType<typename std::add_const<value_type>::type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag >::type
|
||||
typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag, true >::type
|
||||
operator()() const & { return *this; }
|
||||
|
||||
typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag, true >::type
|
||||
operator()() & { return *this; }
|
||||
|
||||
typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag, StorageType::is_weak >::type
|
||||
operator()() && { return *this; }
|
||||
|
||||
#else
|
||||
|
||||
template<typename... Args> // non const version
|
||||
typename boost::lazy_enable_if_c<
|
||||
(!clef::is_any_lazy<Args...>::value) && (indexmaps::slicer<indexmap_type,Args...>::r_type::domain_type::rank!=0)
|
||||
, result_of_call_as_view<false,false,Args...>
|
||||
>::type // enable_if
|
||||
operator()(Args const & ... args) {
|
||||
return typename result_of_call_as_view<false,false,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
|
||||
template<typename... Args> // const version
|
||||
typename boost::lazy_enable_if_c<
|
||||
(!clef::is_any_lazy<Args...>::value) && (indexmaps::slicer<indexmap_type,Args...>::r_type::domain_type::rank!=0)
|
||||
, result_of_call_as_view<true,false,Args...>
|
||||
>::type // enable_if
|
||||
operator()(Args const & ... args) const {
|
||||
return typename result_of_call_as_view<true,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
return typename result_of_call_as_view<true,false,Args...>::type ( indexmaps::slicer<indexmap_type,Args...>::invoke(indexmap_,args...), storage()); }
|
||||
|
||||
typedef typename ViewFactory<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag >::type view_type;
|
||||
|
||||
// Interaction with the CLEF library : calling with any clef expression as argument build a new clef expression
|
||||
/// Equivalent to make_view
|
||||
//typename ISPViewType<typename std::add_const<value_type>::type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag >::type
|
||||
typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag, false >::type
|
||||
operator()() const { return *this; }
|
||||
|
||||
typename ISPViewType<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag, false >::type
|
||||
operator()() { return *this; }
|
||||
|
||||
#endif
|
||||
// Interaction with the CLEF library : calling with any clef expression as argument build a new clef expression
|
||||
// NB : this is ok because indexmap_storage_pair has a shallow copy constructor ...
|
||||
// so A(i_) if A is an array will NOT copy the data....
|
||||
template< typename... Args>
|
||||
@ -223,12 +286,6 @@ namespace triqs { namespace arrays {
|
||||
|
||||
template<typename Fnt> friend void triqs_clef_auto_assign (indexmap_storage_pair & x, Fnt f) { assign_foreach(x,f);}
|
||||
|
||||
/// Equivalent to make_view
|
||||
typename ViewFactory<typename std::add_const<value_type>::type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag >::type
|
||||
operator()() const { return *this; }
|
||||
typename ViewFactory<value_type,domain_type::rank, OptionFlags, TraversalOrder, ViewTag >::type
|
||||
operator()() { return *this; }
|
||||
|
||||
// ------------------------------- Iterators --------------------------------------------
|
||||
typedef iterator_adapter<true,typename IndexMapType::iterator, StorageType> const_iterator;
|
||||
typedef iterator_adapter<false,typename IndexMapType::iterator, StorageType> iterator;
|
||||
|
||||
@ -35,6 +35,7 @@ namespace triqs { namespace arrays {
|
||||
r(0) = A(1)* B(2) - B(1) * A(2);
|
||||
r(1) = - A(0)* B(2) + B(0) * A(2);
|
||||
r(2) = A(0)*B(1) - B(0) * A(1);
|
||||
std::cout << "in cross product "<< A << B << r << std::endl;
|
||||
return r;
|
||||
}
|
||||
|
||||
|
||||
@ -133,8 +133,8 @@ namespace triqs { namespace arrays {
|
||||
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 boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
};
|
||||
|
||||
// The general case
|
||||
@ -179,8 +179,8 @@ namespace triqs { namespace arrays {
|
||||
internal_data(determinant_lazy const & P):M(P.a){det_and_inverse_worker<A_type> worker(M); det = worker.det();}
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
};
|
||||
|
||||
}} // namespace triqs::arrays
|
||||
|
||||
@ -80,7 +80,7 @@ namespace triqs { namespace arrays { namespace linalg {
|
||||
protected:
|
||||
eigenelements_worker_base ( matrix_view <T,Opt> the_matrix) : eigenelements_worker_base <T,Opt,false> (the_matrix) {this->compz='V'; }
|
||||
public:
|
||||
matrix<T,Opt> vectors() const {
|
||||
typename matrix_view<T,Opt>::non_view_type vectors() const {
|
||||
if (!this->has_run) TRIQS_RUNTIME_ERROR<<"eigenelements_worker has not been invoked !";
|
||||
return this->mat;
|
||||
}
|
||||
@ -152,8 +152,10 @@ namespace triqs { namespace arrays { namespace linalg {
|
||||
* if true : a copy is made, M is preserved, but of course it is slower...
|
||||
*/
|
||||
template<typename MatrixViewType >
|
||||
std::pair<array<double,1>, typename MatrixViewType::non_view_type> eigenelements( MatrixViewType M, bool take_copy =false) {
|
||||
eigenelements_worker<MatrixViewType, true> W(take_copy ? make_clone(M)() : M()); W.invoke(); return std::make_pair(W.values(),W.vectors());
|
||||
std::pair<array<double,1>, typename MatrixViewType::non_view_type> eigenelements( MatrixViewType const & M, bool take_copy =false) {
|
||||
eigenelements_worker<typename MatrixViewType::view_type, true> W(take_copy ? MatrixViewType(make_clone(M)) : M);
|
||||
W.invoke();
|
||||
return std::make_pair(W.values(),W.vectors());
|
||||
}
|
||||
|
||||
//--------------------------------
|
||||
@ -167,8 +169,8 @@ namespace triqs { namespace arrays { namespace linalg {
|
||||
* if true : a copy is made, M is preserved, but of course it is slower...
|
||||
*/
|
||||
template<typename MatrixViewType >
|
||||
triqs::arrays::vector_view <double> eigenvalues( MatrixViewType M, bool take_copy = false) {
|
||||
eigenelements_worker<MatrixViewType,false> W(take_copy ? make_clone(M)() : M()); W.invoke(); return W.values();
|
||||
triqs::arrays::vector_view <double> eigenvalues( MatrixViewType const & M, bool take_copy = false) {
|
||||
eigenelements_worker<MatrixViewType,false> W(take_copy ? MatrixViewType(make_clone(M)) : M); W.invoke(); return W.values();
|
||||
}
|
||||
|
||||
}}} // namespace triqs::arrays::linalg
|
||||
|
||||
@ -111,7 +111,7 @@ namespace triqs { namespace arrays {
|
||||
|
||||
// an implementation class to gather the common part to matrix and expression....
|
||||
template<typename A> struct inverse_lazy_impl : TRIQS_MODEL_CONCEPT(ImmutableMatrix) {
|
||||
typedef typename boost::remove_const<typename A::value_type>::type value_type;
|
||||
typedef typename std::remove_const<typename A::value_type>::type value_type;
|
||||
typedef typename A::domain_type domain_type;
|
||||
typedef typename domain_type::index_value_type index_value_type;
|
||||
typedef typename const_view_type_if_exists_else_type<A>::type A_type;
|
||||
@ -133,8 +133,8 @@ namespace triqs { namespace arrays {
|
||||
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 boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
};
|
||||
|
||||
// The general case
|
||||
@ -179,8 +179,8 @@ namespace triqs { namespace arrays {
|
||||
internal_data(determinant_lazy const & P):M(P.a){det_and_inverse_worker<A_type> worker(M); det = worker.det();}
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
};
|
||||
|
||||
}} // namespace triqs::arrays
|
||||
|
||||
@ -58,8 +58,8 @@ namespace triqs { namespace arrays {
|
||||
internal_data(mat_vec_mul_lazy const & P): R(P.M.dim0()) { blas::gemv(1,P.M,P.V,0,R); }
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
|
||||
public:
|
||||
mat_vec_mul_lazy( MT const & M_, VT const & V_):M(M_),V(V_){
|
||||
|
||||
@ -52,8 +52,8 @@ namespace triqs { namespace arrays {
|
||||
internal_data(matmul_lazy const & P): R( P.a.dim0(), P.b.dim1()) { blas::gemm(1.0,P.a, P.b, 0.0, R); }
|
||||
};
|
||||
friend struct internal_data;
|
||||
mutable boost::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= boost::make_shared<internal_data>(*this);}
|
||||
mutable std::shared_ptr<internal_data> _id;
|
||||
void activate() const { if (!_id) _id= std::make_shared<internal_data>(*this);}
|
||||
|
||||
public:
|
||||
matmul_lazy( A const & a_, B const & b_):a(a_),b(b_){
|
||||
@ -83,7 +83,7 @@ namespace triqs { namespace arrays {
|
||||
|
||||
private:
|
||||
template<typename LHS> void assign_comp_impl (LHS & lhs, double S) const {
|
||||
static_assert((is_matrix_or_view<LHS>::value), "LHS is not a matrix");
|
||||
static_assert((is_matrix_or_view<LHS>::value), "LHS is not a matrix");
|
||||
if (lhs.dim0() != dim0())
|
||||
TRIQS_RUNTIME_ERROR<< "Matmul : +=/-= operator : first dimension mismatch in A*B "<< lhs.dim0()<<" vs "<< dim0();
|
||||
if (lhs.dim1() != dim1())
|
||||
|
||||
@ -27,7 +27,7 @@
|
||||
#include "vector.hpp"
|
||||
namespace triqs { namespace arrays {
|
||||
|
||||
template <typename ValueType, ull_t Opt=0, ull_t TraversalOrder= 0> class matrix_view;
|
||||
template <typename ValueType, ull_t Opt=0, ull_t TraversalOrder= 0, bool Borrowed = false> class matrix_view;
|
||||
template <typename ValueType, ull_t Opt=0, ull_t TraversalOrder= 0> class matrix;
|
||||
|
||||
// ---------------------- matrix --------------------------------
|
||||
@ -46,13 +46,14 @@ namespace triqs { namespace arrays {
|
||||
bool memory_layout_is_fortran() const { return this->indexmap().strides()[0] < this->indexmap().strides()[1]; }
|
||||
|
||||
#define IMPL_TYPE indexmap_storage_pair < indexmaps::cuboid::map<2,Opt,TraversalOrder>, \
|
||||
storages::shared_block<ValueType>, Opt, TraversalOrder, Tag::matrix_view >
|
||||
storages::shared_block<ValueType,Borrowed>, Opt, TraversalOrder, Tag::matrix_view >
|
||||
|
||||
template <typename ValueType, ull_t Opt, ull_t TraversalOrder >
|
||||
template <typename ValueType, ull_t Opt, ull_t TraversalOrder, bool Borrowed>
|
||||
class matrix_view : Tag::matrix_view, TRIQS_MODEL_CONCEPT(MutableMatrix), public IMPL_TYPE {
|
||||
public :
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder> view_type;
|
||||
typedef matrix<ValueType,Opt,TraversalOrder> non_view_type;
|
||||
typedef matrix <ValueType,Opt,TraversalOrder> non_view_type;
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder> view_type;
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
@ -96,11 +97,15 @@ namespace triqs { namespace arrays {
|
||||
|
||||
//---------------------------------------------------------------------
|
||||
// this traits is used by indexmap_storage_pair, when slicing to find the correct view type.
|
||||
template < class V, int R, ull_t OptionFlags, ull_t To > struct ViewFactory< V, R, OptionFlags, To, Tag::matrix_view> {
|
||||
typedef typename std::conditional <R == 1, vector_view<V,OptionFlags >, matrix_view<V,OptionFlags, To > >::type type;
|
||||
template < class V, int R, ull_t OptionFlags, ull_t To, bool Borrowed >
|
||||
struct ISPViewType< V, R, OptionFlags, To, Tag::matrix_view, Borrowed> {
|
||||
typedef typename std::conditional <R == 1, vector_view<V,OptionFlags,Borrowed >, matrix_view<V,OptionFlags, To, Borrowed > >::type type;
|
||||
};
|
||||
#undef IMPL_TYPE
|
||||
|
||||
// ---------------------- matrix --------------------------------
|
||||
#define IMPL_TYPE indexmap_storage_pair < indexmaps::cuboid::map<2,Opt,TraversalOrder>, \
|
||||
storages::shared_block<ValueType>, Opt, TraversalOrder, Tag::matrix_view >
|
||||
|
||||
template <typename ValueType, ull_t Opt, ull_t TraversalOrder >
|
||||
class matrix: Tag::matrix, TRIQS_MODEL_CONCEPT(MutableMatrix), public IMPL_TYPE {
|
||||
@ -108,8 +113,9 @@ namespace triqs { namespace arrays {
|
||||
typedef typename IMPL_TYPE::value_type value_type;
|
||||
typedef typename IMPL_TYPE::storage_type storage_type;
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder> view_type;
|
||||
typedef matrix<ValueType,Opt,TraversalOrder> non_view_type;
|
||||
typedef matrix <ValueType,Opt,TraversalOrder> non_view_type;
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder> view_type;
|
||||
typedef matrix_view<ValueType,Opt,TraversalOrder,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
/// Empty matrix.
|
||||
|
||||
@ -52,12 +52,12 @@ namespace triqs { namespace arrays { namespace numpy_interface {
|
||||
}
|
||||
if (!PyArray_Check(res)) TRIQS_RUNTIME_ERROR<<" array_view_from_numpy : internal error : the python object is not a numpy";
|
||||
PyArrayObject * arr = (PyArrayObject *)(res);
|
||||
//PyArray_SetBaseObject(arr, A.storage().new_ref_to_guard());
|
||||
//PyArray_SetBaseObject(arr, A.storage().new_python_ref());
|
||||
#ifdef TRIQS_NUMPY_VERSION_LT_17
|
||||
arr->base = A.storage().new_ref_to_guard();
|
||||
arr->base = A.storage().new_python_ref();
|
||||
assert( arr->flags == (arr->flags & ~NPY_OWNDATA));
|
||||
#else
|
||||
int r = PyArray_SetBaseObject(arr,A.storage().new_ref_to_guard());
|
||||
int r = PyArray_SetBaseObject(arr,A.storage().new_python_ref());
|
||||
if (r!=0) TRIQS_RUNTIME_ERROR << "Internal Error setting the guard in numpy !!!!";
|
||||
assert( PyArray_FLAGS(arr) == (PyArray_FLAGS(arr) & ~NPY_ARRAY_OWNDATA));
|
||||
#endif
|
||||
|
||||
@ -1,68 +0,0 @@
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||
*
|
||||
* Copyright (C) 2011 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
|
||||
* Foundation, either version 3 of the License, or (at your option) any later
|
||||
* version.
|
||||
*
|
||||
* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
|
||||
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
||||
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
|
||||
* details.
|
||||
*
|
||||
* You should have received a copy of the GNU General Public License along with
|
||||
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
******************************************************************************/
|
||||
|
||||
#ifndef TRIQS_ARRAYS_BASIC_BLOCK_H
|
||||
#define TRIQS_ARRAYS_BASIC_BLOCK_H
|
||||
#include <boost/noncopyable.hpp>
|
||||
|
||||
// A minimal basic block : replace by a vector<> or is it quicker ?
|
||||
namespace triqs { namespace arrays { namespace storages { namespace details {
|
||||
|
||||
template<typename ValueType >
|
||||
struct basic_block : boost::noncopyable {
|
||||
typedef typename boost::remove_const<ValueType>::type non_const_value_type;
|
||||
size_t size_;
|
||||
non_const_value_type * restrict p;
|
||||
|
||||
basic_block():size_(0),p(NULL){}
|
||||
|
||||
basic_block (size_t s):size_(s),p(new non_const_value_type[s]) { }
|
||||
|
||||
~basic_block(){if (p) {delete[] p; } }
|
||||
|
||||
void operator=(const basic_block & X) {
|
||||
assert(size_==X.size_);assert(p); assert(X.p);
|
||||
memcpy (p,X.p,size_ * sizeof(ValueType));
|
||||
}
|
||||
|
||||
non_const_value_type & operator[](size_t i) {return p[i];}
|
||||
|
||||
size_t size() const {return size_;}
|
||||
|
||||
template<class Archive>
|
||||
void save(Archive & ar, const unsigned int version) const {
|
||||
ar << boost::serialization::make_nvp("size",size_);
|
||||
for (size_t i=0; i<size_; ++i) ar << boost::serialization::make_nvp("data",p[i]);
|
||||
}
|
||||
|
||||
template<class Archive>
|
||||
void load(Archive & ar, const unsigned int version) {
|
||||
ar >> size_;
|
||||
assert (p==NULL);
|
||||
p = new non_const_value_type[size_];
|
||||
for (size_t i=0; i<size_; ++i) ar >> p[i];
|
||||
}
|
||||
BOOST_SERIALIZATION_SPLIT_MEMBER();
|
||||
};
|
||||
}}}}//namespace triqs::arrays ...
|
||||
#endif
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
|
||||
/*******************************************************************************
|
||||
*
|
||||
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||
@ -19,137 +18,263 @@
|
||||
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
|
||||
*
|
||||
******************************************************************************/
|
||||
#ifndef TRIQS_MEM_BLOCK_H
|
||||
#ifndef TRIQS_MEM_BLOCK_H
|
||||
#define TRIQS_MEM_BLOCK_H
|
||||
#include "../../utility/exceptions.hpp"
|
||||
#include <Python.h>
|
||||
#include <numpy/arrayobject.h>
|
||||
#include "./memcopy.hpp"
|
||||
#include <triqs/utility/macros.hpp>
|
||||
|
||||
//#define TRIQS_ARRAYS_DEBUG_TRACE_MEM
|
||||
#ifdef TRIQS_ARRAYS_DEBUG_TRACE_MEM
|
||||
#include <iostream>
|
||||
#define TRACE_MEM_DEBUG(X) std::cerr<< X << std::endl;
|
||||
#define TRACE_MEM_DEBUG(X) std::cerr<< X << std::endl;
|
||||
#else
|
||||
#define TRACE_MEM_DEBUG(X)
|
||||
#define TRACE_MEM_DEBUG(X)
|
||||
#endif
|
||||
|
||||
namespace triqs { namespace arrays { namespace storages { namespace details {
|
||||
|
||||
template<typename ValueType> class mem_block {
|
||||
// if no python, same code but remove python parts...
|
||||
#ifndef TRIQS_WITH_PYTHON_SUPPORT
|
||||
typedef void PyObject;
|
||||
#define Py_DECREF(P)
|
||||
#endif
|
||||
|
||||
typedef typename boost::remove_const<ValueType>::type non_const_value_type;
|
||||
size_t size_;
|
||||
non_const_value_type * restrict p;
|
||||
PyObject * py_obj;
|
||||
static void import_numpy_array() { if (_import_array()!=0) TRIQS_RUNTIME_ERROR <<"Internal Error in importing numpy";}
|
||||
namespace triqs { namespace arrays { namespace storages { //namespace details {
|
||||
|
||||
public :
|
||||
template<typename ValueType> struct mem_block; // forward
|
||||
|
||||
mem_block():size_(0),p(NULL),py_obj(NULL){}
|
||||
// debug only, to check also weak refs. This will slow down a bit critical loops..
|
||||
// I am not sure this is really useful ...
|
||||
//#define TRIQS_ARRAYS_CHECK_WEAK_REFS
|
||||
|
||||
mem_block (size_t s):size_(s),py_obj(NULL){
|
||||
//p = new non_const_value_type[s]; // check speed penalty for try ??
|
||||
try { p = new non_const_value_type[s];}
|
||||
catch (std::bad_alloc& ba) { TRIQS_RUNTIME_ERROR<< "Memory allocation error : bad_alloc : "<< ba.what();}
|
||||
}
|
||||
// the ref counting functions : weak means : no use of the reference counting system
|
||||
template<bool Weak, typename ValueType> ENABLE_IFC(Weak) inc_ref(mem_block<ValueType> * const m) {
|
||||
#ifdef TRIQS_ARRAYS_CHECK_WEAK_REFS
|
||||
m->weak_ref_count++;
|
||||
#endif
|
||||
}
|
||||
template<bool Weak, typename ValueType> DISABLE_IFC(Weak) inc_ref(mem_block<ValueType> * const m) { m->ref_count++;}
|
||||
|
||||
mem_block (PyObject * obj, bool borrowed ) {
|
||||
TRACE_MEM_DEBUG(" construct memblock from pyobject"<<obj<< " # ref ="<<obj->ob_refcnt<<" borrowed = "<< borrowed);
|
||||
assert(obj); import_numpy_array();
|
||||
if (borrowed) Py_INCREF(obj);
|
||||
if (!PyArray_Check(obj)) TRIQS_RUNTIME_ERROR<<"Internal error : mem_block construct from pyo : obj is not an array";
|
||||
PyArrayObject * arr = (PyArrayObject *)(obj);
|
||||
size_ = PyArray_SIZE(arr);
|
||||
this->p = (non_const_value_type*)PyArray_DATA(arr);
|
||||
py_obj = obj;
|
||||
}
|
||||
// force inline on gcc ?
|
||||
//template<typename ValueType> void dec_ref(mem_block<ValueType> * const m) __attribute__((always_inline));
|
||||
|
||||
~mem_block(){ // delete memory manually iif py_obj is not set. Otherwise the python interpreter will do that for us.
|
||||
TRACE_MEM_DEBUG("deleting mem block p ="<<p<< " py_obj = "<< py_obj << " ref of py obj if exists"<<(py_obj ? py_obj->ob_refcnt: -1));
|
||||
if (py_obj==NULL) { if (p) { delete[] p;} }
|
||||
else Py_DECREF(py_obj);
|
||||
}
|
||||
template<bool Weak, typename ValueType> ENABLE_IFC(Weak) dec_ref(mem_block<ValueType> * const m) {
|
||||
#ifdef TRIQS_ARRAYS_CHECK_WEAK_REFS
|
||||
m->weak_ref_count--;
|
||||
#endif
|
||||
}
|
||||
|
||||
void operator=(const mem_block & X) {
|
||||
//assert( py_obj==NULL);
|
||||
assert(size_==X.size_);assert(p); assert(X.p);
|
||||
storages::memcopy (p, X.p, size_);
|
||||
}
|
||||
template<bool Weak, typename ValueType> DISABLE_IFC(Weak) dec_ref(mem_block<ValueType> * const m) {
|
||||
m->ref_count--; if (m->ref_count ==0) {
|
||||
#ifdef TRIQS_ARRAYS_CHECK_WEAK_REFS
|
||||
std::cout << " detroying "<< m->weak_ref_count <<std::endl;
|
||||
if (m->weak_ref_count !=0) TRIQS_RUNTIME_ERROR << "Deleting an memory block of an array with still "<< m->weak_ref_count<< " weak references";
|
||||
#endif
|
||||
delete m;
|
||||
}
|
||||
}
|
||||
|
||||
mem_block ( mem_block const & X): size_(X.size()), py_obj(NULL) {
|
||||
//p = new non_const_value_type[s]; // check speed penalty for try ??
|
||||
try { p = new non_const_value_type[X.size()];}
|
||||
catch (std::bad_alloc& ba) { TRIQS_RUNTIME_ERROR<< "Memory allocation error : bad_alloc : "<< ba.what();}
|
||||
if ((X.py_obj==NULL) || (PyCObject_Check(X.py_obj))) { (*this) = X; }
|
||||
else {
|
||||
// else make a new copy of the numpy ...
|
||||
import_numpy_array();
|
||||
//assert(PyArray_Check(X.py_obj));
|
||||
if (!is_scalar_or_pod<ValueType>::value) TRIQS_RUNTIME_ERROR << "Internal Error : memcpy on non-scalar";
|
||||
/**
|
||||
* This is a block of memory (pointer p and size size_).
|
||||
* INTERNAL USE only, by shared_block only
|
||||
*
|
||||
* The memory can be :
|
||||
*
|
||||
* - allocated (and deleted in C++)
|
||||
* - owned by a numpy python object (py_numpy)
|
||||
*
|
||||
* The block contains its own reference system, to avoid the use of shared_ptr in shared_block
|
||||
* (which was very slow in critical codes).
|
||||
*
|
||||
* In addition, the python guard system can return to python an array allocated in C++.
|
||||
*
|
||||
* The memory block has 4 possible states (and only 4) :
|
||||
*
|
||||
* * State 0) : p ==nullptr && py_numpy == nullptr && py_guard == nullptr
|
||||
* Default state. Also obtained when the object has been moved from
|
||||
*
|
||||
* * State 1) : p !=nullptr && py_numpy == nullptr && py_guard == nullptr
|
||||
* Memory block allocated and used in C++.
|
||||
*
|
||||
* * State 2) : p !=nullptr && py_numpy != nullptr && py_guard == nullptr
|
||||
* Memory block was allocated by python, and is used in C++.
|
||||
* The block keeps a python ref (owned!) in py_numpy which is released at destruction
|
||||
*
|
||||
* * State 3) : p !=nullptr && py_numpy == nullptr && py_guard != nullptr
|
||||
* Memory block allocated by C++, but passed to python.
|
||||
* The guard is an owned reference to a python object, which itself owns a c++ reference to the mem_block
|
||||
* When the python is done with this object, it releases this reference and the mem_block
|
||||
* can be destroyed (when nobody else uses it).
|
||||
* This guarantees that :
|
||||
* * as long as python does not clean the guard (which is then later attached to a numpy array, cf numpy interface)
|
||||
* the block will not be deleted.
|
||||
* * when python is done with this numpy, hence the guard, the c++ reference is dec_refed and
|
||||
* the usage can continue normally in c++ (without *any* python ref contrary to a previous design).
|
||||
*
|
||||
* * Invariants :
|
||||
* * py_numpy == nullptr || py_guard == nullptr :
|
||||
* * ref_count >=1.
|
||||
*/
|
||||
template<typename ValueType> struct mem_block {
|
||||
|
||||
size_t size_; // size of the block
|
||||
ValueType * restrict p; // the memory block. Owned by this, except when py_numpy is not null
|
||||
size_t ref_count; // number of refs. : >=1
|
||||
size_t weak_ref_count; // number of refs. : >=1
|
||||
PyObject * py_numpy; // if not null, an owned reference to a numpy which is the data of this block
|
||||
PyObject * py_guard; // if not null, a BORROWED reference to the guard. If null, the guard does not exist
|
||||
static_assert(!std::is_const<ValueType>::value, "internal error");
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
static void import_numpy_array() { if (_import_array()!=0) TRIQS_RUNTIME_ERROR <<"Internal Error in importing numpy";}
|
||||
#endif
|
||||
|
||||
//Construct to state 0
|
||||
mem_block():size_(0),p(nullptr),py_numpy(nullptr), py_guard(nullptr), ref_count(1), weak_ref_count(0){}
|
||||
|
||||
// construct to state 1 with a given size.
|
||||
mem_block (size_t s):size_(s),py_numpy(nullptr), py_guard(nullptr){
|
||||
try { p = new ValueType[s];}
|
||||
catch (std::bad_alloc& ba) { TRIQS_RUNTIME_ERROR<< "Memory allocation error : bad_alloc : "<< ba.what();}
|
||||
ref_count=1;
|
||||
weak_ref_count =0;
|
||||
}
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
// construct to state 2. python_object_is_borrowed : whether the python ref is borrowed
|
||||
mem_block (PyObject * obj, bool python_object_is_borrowed) {
|
||||
TRACE_MEM_DEBUG(" construct memblock from pyobject"<<obj<< " # ref ="<<obj->ob_refcnt<<" python_object_is_borrowed = "<< python_object_is_borrowed);
|
||||
assert(obj); import_numpy_array();
|
||||
if (python_object_is_borrowed) Py_INCREF(obj);
|
||||
if (!PyArray_Check(obj)) TRIQS_RUNTIME_ERROR<<"Internal error : mem_block construct from pyo : obj is not an array";
|
||||
PyArrayObject * arr = (PyArrayObject *)(obj);
|
||||
size_ = PyArray_SIZE(arr);
|
||||
p = (ValueType*)PyArray_DATA(arr);
|
||||
py_numpy = obj;
|
||||
py_guard = nullptr;
|
||||
ref_count=1;
|
||||
weak_ref_count =0;
|
||||
}
|
||||
#endif
|
||||
|
||||
// destructor : release memory only in state 1. This should NEVER be called in state 3 (first need to get back to 1).
|
||||
~mem_block(){ // delete memory manually iif py_obj is not set. Otherwise the python interpreter will do that for us.
|
||||
TRACE_MEM_DEBUG("deleting mem block p ="<<p<< " py_obj = "<< py_obj << " ref of py obj if exists"<<(py_obj ? py_obj->ob_refcnt: -1));
|
||||
assert(ref_count<=1); assert(py_guard==nullptr);// state 3 forbidden
|
||||
if (py_numpy) Py_DECREF(py_numpy); // state 1
|
||||
else { if (p) delete[] p; } // state 2 or state 0
|
||||
}
|
||||
|
||||
// can not be copied or moved.
|
||||
mem_block & operator=(mem_block const & X) = delete;
|
||||
mem_block & operator=(mem_block && X) = delete;
|
||||
mem_block(mem_block && X) noexcept {
|
||||
size_ = X.size_; p = X.p; ref_count = X.ref_count; weak_ref_count = X.weak_ref_count; py_numpy=X.py_numpy; py_guard = X.py_guard;
|
||||
X.p =nullptr; X.py_numpy= nullptr; X.py_guard = nullptr; // state 0, ready to destruct
|
||||
}
|
||||
|
||||
// deep copy of data from another block.
|
||||
// MUST be of the same size...
|
||||
void copy_from(const mem_block & X) {
|
||||
assert(size_==X.size_);assert(p); assert(X.p);
|
||||
storages::memcopy (p, X.p, size_);
|
||||
}
|
||||
|
||||
// copy construct into state 1, always.
|
||||
// This is a choice, even if X is state 2 (a numpy).
|
||||
// We copy a numpy into a regular C++ array, which can then be used at max speed.
|
||||
mem_block (mem_block const & X): size_(X.size()), py_numpy(nullptr), py_guard(nullptr) {
|
||||
try { p = new ValueType[X.size()];}
|
||||
catch (std::bad_alloc& ba) { TRIQS_RUNTIME_ERROR<< "Memory allocation error : bad_alloc : "<< ba.what();}
|
||||
ref_count=1;
|
||||
weak_ref_count =0;
|
||||
// now we copy the data
|
||||
#ifndef TRIQS_WITH_PYTHON_SUPPORT
|
||||
copy_from(X);
|
||||
#else
|
||||
// if X is in state 1 or 3
|
||||
if (X.py_numpy==nullptr) { copy_from(X); }
|
||||
else { // X was in state 2
|
||||
// else make a new copy of the numpy ...
|
||||
import_numpy_array();
|
||||
if (!is_scalar_or_pod<ValueType>::value) TRIQS_RUNTIME_ERROR << "Internal Error : memcpy on non-scalar";
|
||||
#ifdef TRIQS_NUMPY_VERSION_LT_17
|
||||
if ( ( PyArray_ISFORTRAN(X.py_obj)) || (PyArray_ISCONTIGUOUS(X.py_obj))) {
|
||||
memcpy (p,PyArray_DATA(X.py_obj),size_ * sizeof(ValueType));
|
||||
}
|
||||
PyObject * arr3 = X.py_numpy;
|
||||
#else
|
||||
// STRANGE : uncommenting this leads to a segfault on mac ???
|
||||
// TO BE INVESTIGATED, IT IS NOT NORMAL
|
||||
//if (!PyArray_Check(X.py_obj)) TRIQS_RUNTIME_ERROR<<"Internal error : is not an array";
|
||||
PyArrayObject * arr3 = (PyArrayObject *)(X.py_obj);
|
||||
if ( ( PyArray_ISFORTRAN(arr3)) || (PyArray_ISCONTIGUOUS(arr3))) {
|
||||
memcpy (p,PyArray_DATA(arr3),size_ * sizeof(ValueType));
|
||||
}
|
||||
// STRANGE : uncommenting this leads to a segfault on mac ???
|
||||
// TO BE INVESTIGATED, IT IS NOT NORMAL
|
||||
//if (!PyArray_Check(X.py_numpy)) TRIQS_RUNTIME_ERROR<<"Internal error : is not an array";
|
||||
PyArrayObject * arr3 = (PyArrayObject *)(X.py_numpy);
|
||||
#endif
|
||||
else { // if the X.py_obj is not contiguous, first let numpy copy it properly
|
||||
PyObject * na = PyObject_CallMethod(X.py_obj,(char *)"copy",NULL);
|
||||
assert(na);
|
||||
// if we can make a memcpy, do it.
|
||||
if ( ( PyArray_ISFORTRAN(arr3)) || (PyArray_ISCONTIGUOUS(arr3))) {
|
||||
memcpy (p,PyArray_DATA(arr3),size_ * sizeof(ValueType));
|
||||
}
|
||||
else { // if the X.py_numpy is not contiguous, first let numpy copy it properly, then memcpy
|
||||
PyObject * na = PyObject_CallMethod(X.py_numpy,(char *)"copy",nullptr);
|
||||
assert(na);
|
||||
#ifdef TRIQS_NUMPY_VERSION_LT_17
|
||||
assert( ( PyArray_ISFORTRAN(na)) || (PyArray_ISCONTIGUOUS(na)));
|
||||
memcpy (p,PyArray_DATA(na),size_ * sizeof(ValueType));
|
||||
PyObject * arr = na;
|
||||
#else
|
||||
if (!PyArray_Check(na)) TRIQS_RUNTIME_ERROR<<"Internal error : is not an array";
|
||||
PyArrayObject * arr = (PyArrayObject *)(na);
|
||||
assert( ( PyArray_ISFORTRAN(arr)) || (PyArray_ISCONTIGUOUS(arr)));
|
||||
memcpy (p,PyArray_DATA(arr),size_ * sizeof(ValueType));
|
||||
if (!PyArray_Check(na)) TRIQS_RUNTIME_ERROR<<"Internal error : is not an array";
|
||||
PyArrayObject * arr = (PyArrayObject *)(na);
|
||||
#endif
|
||||
Py_DECREF(na);
|
||||
}
|
||||
assert( ( PyArray_ISFORTRAN(arr)) || (PyArray_ISCONTIGUOUS(arr)));
|
||||
memcpy (p,PyArray_DATA(arr),size_ * sizeof(ValueType));
|
||||
Py_DECREF(na);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
static void delete_pointeur( void *ptr ) {
|
||||
TRACE_MEM_DEBUG("deleting data block"<<(non_const_value_type*) ptr);
|
||||
delete [] ( (non_const_value_type*) ptr) ;
|
||||
}
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
// Precondition : state 3, postcondition: state 1
|
||||
static void delete_python_guard(void *ptr ) {
|
||||
mem_block * m = static_cast<mem_block*>(ptr);
|
||||
TRACE_MEM_DEBUG("deleting data block"<<m);
|
||||
assert(m->ref_count>0); assert(m->py_guard !=nullptr);assert(m->py_numpy ==nullptr);
|
||||
m->py_guard=nullptr;
|
||||
dec_ref<false>(m);// release the reference which was owned by the guard
|
||||
}
|
||||
|
||||
PyObject * new_ref_to_guard() {
|
||||
if (py_obj==NULL) {
|
||||
TRACE_MEM_DEBUG(" activating python guard for C++ block"<<p<< " py_obj = "<< py_obj);
|
||||
py_obj = PyCObject_FromVoidPtr( (void*) p, &mem_block<ValueType>::delete_pointeur);
|
||||
}
|
||||
Py_INCREF(py_obj);
|
||||
return py_obj;
|
||||
}
|
||||
// returns a NEW python ref either to the numpy or to the guard.
|
||||
// if the guard does not yet exists, create it.
|
||||
PyObject * new_python_ref() {
|
||||
// if we just have a borrowed numpy, just return a new ref to it
|
||||
if (py_numpy) { Py_INCREF(py_numpy); return py_numpy;}
|
||||
// if the guard is already set, then return it, otherwise create it...
|
||||
if (py_guard) { Py_INCREF(py_guard); return py_guard;}
|
||||
else {
|
||||
TRACE_MEM_DEBUG(" activating python guard for C++ block"<<p<< " py_guard = "<< py_guard);
|
||||
inc_ref<false>(this); // the guard owns a C++ reference !!
|
||||
py_guard = PyCObject_FromVoidPtr( static_cast<void *>(this), &mem_block<ValueType>::delete_python_guard);
|
||||
}
|
||||
return py_guard;
|
||||
}
|
||||
#endif
|
||||
|
||||
non_const_value_type & operator[](size_t i) {return p[i];}
|
||||
ValueType & operator[](size_t i) {return p[i];}
|
||||
|
||||
size_t size() const {return size_;}
|
||||
size_t size() const {return size_;}
|
||||
|
||||
template<class Archive>
|
||||
void save(Archive & ar, const unsigned int version) const {
|
||||
ar << boost::serialization::make_nvp("size",size_);
|
||||
for (size_t i=0; i<size_; ++i) ar << boost::serialization::make_nvp("data",p[i]);
|
||||
}
|
||||
template<class Archive>
|
||||
void save(Archive & ar, const unsigned int version) const {
|
||||
ar << boost::serialization::make_nvp("size",size_);
|
||||
for (size_t i=0; i<size_; ++i) ar << boost::serialization::make_nvp("data",p[i]);
|
||||
}
|
||||
|
||||
template<class Archive>
|
||||
void load(Archive & ar, const unsigned int version) {
|
||||
ar >> size_;
|
||||
assert (p==NULL);
|
||||
p = new non_const_value_type[size_];
|
||||
for (size_t i=0; i<size_; ++i) ar >> p[i];
|
||||
}
|
||||
BOOST_SERIALIZATION_SPLIT_MEMBER();
|
||||
template<class Archive>
|
||||
void load(Archive & ar, const unsigned int version) {
|
||||
ar >> size_;
|
||||
assert (p==nullptr);
|
||||
p = new ValueType[size_];
|
||||
for (size_t i=0; i<size_; ++i) ar >> p[i];
|
||||
}
|
||||
BOOST_SERIALIZATION_SPLIT_MEMBER();
|
||||
};
|
||||
}}}}//namespace triqs::arrays
|
||||
|
||||
}}}//namespace triqs::arrays
|
||||
|
||||
#undef TRACE_MEM_DEBUG
|
||||
#endif
|
||||
|
||||
@ -22,128 +22,89 @@
|
||||
#define TRIQS_STORAGE_SHARED_POINTER_H
|
||||
#include <string.h>
|
||||
#include <limits>
|
||||
#include <boost/shared_ptr.hpp>
|
||||
#include <boost/type_traits/is_const.hpp>
|
||||
#include <boost/type_traits/add_const.hpp>
|
||||
#include <boost/type_traits/remove_const.hpp>
|
||||
#include <boost/type_traits/is_same.hpp>
|
||||
#include <boost/serialization/shared_ptr.hpp>
|
||||
#include "../impl/make_const.hpp"
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
#include "./mem_block.hpp"
|
||||
#else
|
||||
#include <vector>
|
||||
//#include "basic_block.hpp"
|
||||
|
||||
namespace triqs { namespace arrays { namespace storages {
|
||||
|
||||
/*
|
||||
* This is a shared memory block, similar to a shared_ptr<mem_block>
|
||||
* but a lot faster and adapted to model the Storage concept.
|
||||
*/
|
||||
template<typename ValueType, bool Weak=false>
|
||||
class shared_block {
|
||||
static_assert( !std::is_const<ValueType>::value, "oops : internal error : shared_block should never be instanciated with a const");
|
||||
|
||||
mem_block<ValueType> * sptr;
|
||||
ValueType * restrict data_; // for optimization on some compilers. ?? obsolete : to be removed ?
|
||||
size_t s;
|
||||
|
||||
void construct_delegate(size_t size) {
|
||||
s = size;
|
||||
if (size ==0) { sptr = nullptr; data_=nullptr; }
|
||||
else { sptr = new mem_block<ValueType>(size); data_ = sptr->p;}
|
||||
}
|
||||
|
||||
public:
|
||||
typedef ValueType value_type;
|
||||
static constexpr bool is_weak = Weak;
|
||||
|
||||
/// Construct a new block of memory of given size
|
||||
explicit shared_block(size_t size, Tag::no_init) { construct_delegate(size); }
|
||||
|
||||
explicit shared_block(size_t size, Tag::default_init) {// C++11 : delegate to previous constructor when gcc 4.6 support is out.
|
||||
construct_delegate(size);
|
||||
const auto s = this->size(); for (size_t u=0; u<s; ++u) data_[u] = ValueType();
|
||||
}
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
explicit shared_block(PyObject * arr, bool weak): sptr(new mem_block<ValueType>(arr,weak)) { data_ = sptr->p; s= sptr->size(); }
|
||||
#endif
|
||||
|
||||
namespace triqs { namespace arrays {
|
||||
namespace Tag {struct storage{}; struct shared_block:storage{}; }
|
||||
explicit shared_block() { sptr =nullptr; data_=nullptr; s=0; }
|
||||
|
||||
namespace storages {
|
||||
shared_block(shared_block const & X) noexcept { sptr =X.sptr; data_ = X.data_; s= X.s; if (sptr) inc_ref<Weak>(sptr); }
|
||||
shared_block(shared_block<ValueType,!Weak> const & X) noexcept { sptr =X.sptr; data_ = X.data_; s= X.s; if (sptr) inc_ref<Weak>(sptr); }
|
||||
|
||||
template <class V, class Opt> struct __init_value;
|
||||
template <class V> struct __init_value<V, Tag::no_init> { static void invoke (V * data, size_t s) {}};
|
||||
template <class V> struct __init_value< V, Tag::default_init> {
|
||||
static void invoke (V * restrict data, size_t s) { if (data!=NULL) for (size_t u=0; u<s; ++u) data[u] = V(); }
|
||||
};
|
||||
template <class V> struct __init_value< V, Tag::nan_inf_init> {
|
||||
static void invoke (V * restrict data, size_t s) {
|
||||
static_assert ( ( std::numeric_limits<V>::has_quiet_NaN || std::numeric_limits<V>::is_integer), "type has no Nan and is not integer. This init is impossible");
|
||||
if (data==NULL) return;
|
||||
if (std::numeric_limits<V>::has_quiet_NaN) for (size_t u=0; u<s; ++u) data[u] = std::numeric_limits<V>::quiet_NaN();
|
||||
if (std::numeric_limits<V>::is_integer) for (size_t u=0; u<s; ++u) data[u] = std::numeric_limits<V>::max();
|
||||
shared_block(shared_block && X) noexcept { sptr =X.sptr; data_ = X.data_; s= X.s; X.sptr=nullptr; }
|
||||
|
||||
~shared_block() { if (sptr) dec_ref<Weak>(sptr); }
|
||||
|
||||
shared_block & operator=(shared_block const & X) noexcept {
|
||||
if (sptr) dec_ref<Weak>(sptr); sptr =X.sptr; if (sptr) inc_ref<Weak>(sptr);
|
||||
data_ = X.data_; s = X.s;
|
||||
return *this;
|
||||
}
|
||||
|
||||
shared_block & operator=(shared_block && X) noexcept {
|
||||
if (sptr) dec_ref<Weak>(sptr); sptr =X.sptr; X.sptr=nullptr;
|
||||
data_ = X.data_; s = X.s;
|
||||
return *this;
|
||||
}
|
||||
|
||||
/// True copy of the data
|
||||
shared_block clone() const {
|
||||
shared_block res;
|
||||
if (!empty()) { res.sptr = new mem_block<ValueType> (*sptr); res.data_ = res.sptr->p; res.s = s;}
|
||||
return res;
|
||||
}
|
||||
|
||||
value_type & operator[](size_t i) const { return data_[i];}
|
||||
bool empty() const {return (sptr==nullptr);}
|
||||
size_t size() const {return s;}
|
||||
//size_t size() const {return (empty () ? 0 : sptr->size());}
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
PyObject * new_python_ref() const {return sptr->new_python_ref();}
|
||||
#endif
|
||||
|
||||
private:
|
||||
friend class shared_block<ValueType,!Weak>;
|
||||
friend class boost::serialization::access;
|
||||
template<class Archive> void serialize(Archive & ar, const unsigned int version) {
|
||||
ar & boost::serialization::make_nvp("ptr",sptr); data_ = (sptr ? sptr->p : nullptr); s = (sptr ? sptr->size() : 0);
|
||||
}
|
||||
};
|
||||
|
||||
template <class V> struct __init_value< std::complex<V>, Tag::nan_inf_init> {
|
||||
static void invoke (std::complex<V> * restrict data, size_t s) {
|
||||
static_assert ( ( std::numeric_limits<V>::has_quiet_NaN || std::numeric_limits<V>::is_integer), "type has no Nan and is not integer. This init is impossible");
|
||||
if (data==NULL) return;
|
||||
if (std::numeric_limits<V>::has_quiet_NaN) for (size_t u=0; u<s; ++u) data[u] = std::complex<V>(std::numeric_limits<V>::quiet_NaN(),std::numeric_limits<V>::quiet_NaN()) ;
|
||||
if (std::numeric_limits<V>::is_integer) for (size_t u=0; u<s; ++u) data[u] = std::complex<V>(std::numeric_limits<V>::max(),std::numeric_limits<V>::max());
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/* Storage as a shared_ptr to a basic_block
|
||||
* The shared pointer guarantees that the data will not be destroyed during the life of the array.
|
||||
* Impl: we are not using shared_array directly because of serialization support for shared_ptr */
|
||||
template<typename ValueType >
|
||||
class shared_block : Tag::shared_block {
|
||||
typedef typename boost::add_const<ValueType>::type const_value_type;
|
||||
typedef typename boost::remove_const<ValueType>::type non_const_value_type;
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
typedef details::mem_block<non_const_value_type> block_type;
|
||||
#else
|
||||
// typedef details::basic_block<non_const_value_type> block_type;
|
||||
typedef std::vector<non_const_value_type> block_type;
|
||||
#endif
|
||||
|
||||
public:
|
||||
typedef triqs::arrays::Tag::shared_block tag;
|
||||
typedef ValueType value_type;
|
||||
typedef shared_block<value_type> clone_type;
|
||||
typedef shared_block<const_value_type> const_clone_type;
|
||||
|
||||
/// Construct a new block of memory of given size
|
||||
template<typename InitOpt>
|
||||
explicit shared_block(size_t size, InitOpt opt ): sptr(size ? new block_type(size) : NULL) {
|
||||
init_data();
|
||||
__init_value<value_type,InitOpt>::invoke (this-> data_, this->size());
|
||||
}
|
||||
|
||||
explicit shared_block(): sptr() { init_data(); }
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
/// Construct from a numpy object
|
||||
explicit shared_block(PyObject * arr, bool borrowed): sptr(new block_type(arr,borrowed)) { init_data();}
|
||||
#endif
|
||||
|
||||
/// Shallow copy
|
||||
shared_block(const shared_block<const_value_type> & X): sptr(X.sptr) { init_data(); }
|
||||
|
||||
/// Shallow copy
|
||||
shared_block(const shared_block<non_const_value_type> & X): sptr(X.sptr) { init_data(); }
|
||||
|
||||
void operator=(const shared_block & X) { sptr=X.sptr; init_data(); }
|
||||
|
||||
/// True copy of the data
|
||||
clone_type clone() const {
|
||||
if (empty()) return clone_type ();
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
clone_type res; res.sptr = boost::make_shared<block_type > (*sptr); res.init_data();
|
||||
#else
|
||||
clone_type res(this->size(), Tag::no_init() ); (*res.sptr) = (*sptr);
|
||||
#endif
|
||||
return res;
|
||||
}
|
||||
|
||||
// Make a clone forced to have const value_type
|
||||
const_clone_type const_clone() const {return clone();}
|
||||
|
||||
value_type & operator[](size_t i) const { return data_[i];}
|
||||
bool empty() const {return (sptr.get()==NULL);}
|
||||
size_t size() const {return (empty () ? 0 : sptr.get()->size());}
|
||||
|
||||
#ifdef TRIQS_WITH_PYTHON_SUPPORT
|
||||
PyObject * new_ref_to_guard() const {return sptr->new_ref_to_guard();}
|
||||
#endif
|
||||
|
||||
private:
|
||||
boost::shared_ptr<block_type > sptr;
|
||||
value_type * restrict data_; // for optimization on some compilers.
|
||||
void init_data(){ data_ = (sptr ? &((*sptr)[0]) : NULL); }
|
||||
friend class shared_block <non_const_value_type>; friend class shared_block <const_value_type>;
|
||||
friend class boost::serialization::access;
|
||||
template<class Archive> void serialize(Archive & ar, const unsigned int version) { ar & boost::serialization::make_nvp("ptr",sptr); init_data(); }
|
||||
};
|
||||
}
|
||||
|
||||
namespace details {
|
||||
template<bool Const, typename T> struct make_const_type<Const,storages::shared_block<T> > {
|
||||
typedef storages::shared_block<typename make_const_type<Const,T>::type> type;
|
||||
};
|
||||
}
|
||||
}}//namespace triqs::arrays
|
||||
}}}//namespace
|
||||
#endif
|
||||
|
||||
|
||||
@ -29,19 +29,20 @@
|
||||
|
||||
namespace triqs { namespace arrays {
|
||||
|
||||
template <typename ValueType, ull_t Opt=0> class vector_view;
|
||||
template <typename ValueType, ull_t Opt=0, bool Borrowed =false> class vector_view;
|
||||
template <typename ValueType, ull_t Opt=0> class vector;
|
||||
|
||||
// ---------------------- vector_view --------------------------------
|
||||
|
||||
#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<1,Opt,0> , storages::shared_block<ValueType>, Opt, 0, Tag::vector_view >
|
||||
#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<1,Opt,0> , storages::shared_block<ValueType,Borrowed>, Opt, 0, Tag::vector_view >
|
||||
|
||||
/** */
|
||||
template <typename ValueType, ull_t Opt >
|
||||
template <typename ValueType, ull_t Opt, bool Borrowed>
|
||||
class vector_view : Tag::vector_view, TRIQS_MODEL_CONCEPT(MutableVector), public IMPL_TYPE {
|
||||
public :
|
||||
typedef vector_view<ValueType,Opt> view_type;
|
||||
typedef vector<ValueType,Opt> non_view_type;
|
||||
typedef vector <ValueType,Opt> non_view_type;
|
||||
typedef vector_view<ValueType,Opt,false> view_type;
|
||||
typedef vector_view<ValueType,Opt,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
@ -93,9 +94,12 @@ namespace triqs { namespace arrays {
|
||||
//template<typename Arg> auto operator[](Arg && arg) const DECL_AND_RETURN((*this)(std::forward<Arg>(arg)));
|
||||
//template<typename Arg> auto operator[](Arg && arg) DECL_AND_RETURN((*this)(std::forward<Arg>(arg)));
|
||||
};
|
||||
#undef IMPL_TYPE
|
||||
|
||||
template < class V, int R, ull_t OptionFlags, ull_t To > struct ViewFactory< V, R,OptionFlags,To, Tag::vector_view> { typedef vector_view<V,OptionFlags> type; };
|
||||
template < class V, int R, ull_t OptionFlags, ull_t To, bool Borrowed >
|
||||
struct ISPViewType< V, R,OptionFlags,To, Tag::vector_view, Borrowed> { typedef vector_view<V,OptionFlags,Borrowed> type; };
|
||||
// ---------------------- vector--------------------------------
|
||||
#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<1,Opt,0> , storages::shared_block<ValueType>, Opt, 0, Tag::vector_view >
|
||||
|
||||
template <typename ValueType, ull_t Opt>
|
||||
class vector: Tag::vector, TRIQS_MODEL_CONCEPT(MutableVector), public IMPL_TYPE {
|
||||
@ -103,8 +107,9 @@ namespace triqs { namespace arrays {
|
||||
typedef typename IMPL_TYPE::value_type value_type;
|
||||
typedef typename IMPL_TYPE::storage_type storage_type;
|
||||
typedef typename IMPL_TYPE::indexmap_type indexmap_type;
|
||||
typedef vector_view<ValueType,Opt> view_type;
|
||||
typedef vector<ValueType,Opt> non_view_type;
|
||||
typedef vector <ValueType,Opt> non_view_type;
|
||||
typedef vector_view<ValueType,Opt> view_type;
|
||||
typedef vector_view<ValueType,Opt,true> weak_view_type;
|
||||
typedef void has_view_type_tag;
|
||||
|
||||
/// Empty vector.
|
||||
@ -277,8 +282,8 @@ namespace triqs { namespace arrays {
|
||||
}
|
||||
|
||||
// swapping 2 vector
|
||||
template <typename V, ull_t S1, ull_t S2>
|
||||
void deep_swap(vector_view <V,S1> x, vector_view<V,S2> y) {
|
||||
template <typename V, ull_t S1, ull_t S2, bool B1, bool B2>
|
||||
void deep_swap(vector_view <V,S1,B1> x, vector_view<V,S2,B2> y) {
|
||||
blas::swap(x,y);
|
||||
}
|
||||
|
||||
|
||||
@ -65,7 +65,10 @@ namespace triqs { namespace lattice_tools {
|
||||
|
||||
brillouin_zone::brillouin_zone( bravais_lattice const & bl_) : lattice_(bl_), K_reciprocal(3,3) {
|
||||
bravais_lattice::units_type Units(lattice().units());
|
||||
std::cout << Units << std::endl;
|
||||
double delta = dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range())));
|
||||
std::cout << dot(Units(0,range()), cross_product(Units(1,range()),Units(2,range())))<<std::endl;
|
||||
std::cout << cross_product(Units(1,range()),Units(2,range()))<<std::endl;
|
||||
if (abs(delta)<almost_zero) TRIQS_RUNTIME_ERROR<<"Tiling : the 3 vectors of Units are not independant";
|
||||
K_reciprocal(0,range()) = cross_product(Units(1,range()),Units(2,range())) / delta;
|
||||
K_reciprocal(1,range()) = cross_product(Units(2,range()),Units(0,range())) / delta;
|
||||
|
||||
@ -24,6 +24,7 @@
|
||||
#include <triqs/arrays/matrix.hpp>
|
||||
#include <triqs/arrays/vector.hpp>
|
||||
#include <triqs/utility/exceptions.hpp>
|
||||
#include <unordered_map>
|
||||
#include <map>
|
||||
#include <string>
|
||||
#include <unordered_map>
|
||||
|
||||
@ -22,8 +22,17 @@
|
||||
#define TRIQS_COMPILER_DETAILS_H
|
||||
|
||||
#ifdef __GNUC__
|
||||
#define restrict __restrict__
|
||||
#define restrict __restrict__
|
||||
|
||||
#define GCC_VERSION (__GNUC__ * 10000 \
|
||||
+ __GNUC_MINOR__ * 100 \
|
||||
+ __GNUC_PATCHLEVEL__)
|
||||
#if GCC_VERSION < 40801
|
||||
#define TRIQS_COMPILER_IS_OBSOLETE
|
||||
// we still accept gcc 4.6 and 4.7, but only the 4.8.1 and higher is a compliant c++11
|
||||
#endif
|
||||
#endif
|
||||
|
||||
|
||||
#ifdef __GXX_EXPERIMENTAL_CXX0X__
|
||||
#define TRIQS_USE_STATIC_ASSERT
|
||||
|
||||
Loading…
Reference in New Issue
Block a user