From 47cb8a03f7da19e85f87be1732874ca75aea92c4 Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Sun, 7 Sep 2014 14:34:10 +0200 Subject: [PATCH] [arrays] Important changes in implementation. - Simplify group_indices - Only for C ordered, remove complex compile time. - Could be generalized to non C ordered, but no need. - Fix slice for custom orders. - Generalize the group_indices for the custom order. - Add c_ordered_transposed_view (useful ?) - Improve slice, special for ellipsis (quicker). - Simplify TraversalOrder - Assignement. Specialize one case for speed. - use FORCEINLINE in foreach, according to speed test for clang - add one speed test - Modify iterators for better speed. - along the lines decided for the foreach - update doc. --- doc/reference/arrays/algebras_1.cpp | 2 +- .../arrays/containers/reg_constructors_0.cpp | 12 +- doc/reference/arrays/containers/regular.rst | 17 +- doc/reference/arrays/containers/resize_0.cpp | 2 +- doc/reference/arrays/containers/stream.rst | 12 +- .../arrays/containers/template_parameters.rst | 60 +-- .../arrays/containers/view_constructors.rst | 2 +- doc/reference/arrays/containers/views.rst | 21 +- doc/reference/arrays/foreach.rst | 11 +- test/pytriqs/arrays/f.hpp | 2 +- test/speed/array1.cpp | 8 +- test/speed/array3.cpp | 160 +++++++ test/speed/array_lazy_assign.cpp | 4 +- test/speed/array_short1.cpp | 8 +- test/speed/gf_s.cpp | 2 +- test/triqs/arrays/bound_check.cpp | 4 +- .../arrays/c_ordered_transposed_view.cpp | 36 ++ test/triqs/arrays/common.hpp | 1 + test/triqs/arrays/const_view.cpp | 2 +- test/triqs/arrays/ellipsis.cpp | 1 + test/triqs/arrays/fill_with_lazy.cpp | 2 +- test/triqs/arrays/group_indices.cpp | 138 +++--- test/triqs/arrays/group_indices.output | 46 -- test/triqs/arrays/iterators.cpp | 2 +- test/triqs/arrays/reinterpret_view.cpp | 30 +- test/triqs/arrays/slice_index_order.cpp | 135 +++--- test/triqs/arrays/views3.cpp | 70 ++- test/triqs/utility/tuple_tools.cpp | 4 +- test/triqs/utility/view_tools.cpp | 4 +- triqs/arrays.hpp | 5 + triqs/arrays/array.hpp | 131 ++--- triqs/arrays/blas_lapack/gemm.hpp | 4 +- triqs/arrays/blas_lapack/gemv.hpp | 4 +- triqs/arrays/blas_lapack/ger.hpp | 4 +- triqs/arrays/cache.hpp | 6 +- triqs/arrays/group_indices.hpp | 77 +++ triqs/arrays/impl/assignment.hpp | 14 +- triqs/arrays/impl/flags.hpp | 87 ---- triqs/arrays/impl/indexmap_storage_pair.hpp | 108 ++--- triqs/arrays/impl/print.hpp | 77 +++ triqs/arrays/indexmaps/cuboid/domain.hpp | 258 ++++------ triqs/arrays/indexmaps/cuboid/foreach.hpp | 139 +++++- .../arrays/indexmaps/cuboid/group_indices.hpp | 107 ----- triqs/arrays/indexmaps/cuboid/map.hpp | 164 ++++--- triqs/arrays/indexmaps/cuboid/mem_layout.hpp | 185 +++---- triqs/arrays/indexmaps/cuboid/slice.hpp | 176 ++++--- .../cuboid/slice_traversal_order.hpp | 48 -- triqs/arrays/linalg/det_and_inverse.hpp | 2 +- triqs/arrays/matrix.hpp | 361 +++++++------- triqs/arrays/matrix_tensor_proxy.hpp | 2 + triqs/arrays/python/numpy_extractor.cpp | 1 - triqs/arrays/storages/shared_block.hpp | 10 +- triqs/arrays/vector.hpp | 450 +++++++++--------- triqs/gfs/data_proxies.hpp | 26 +- triqs/utility/mini_vector.hpp | 148 +++++- triqs/utility/mpi1.hpp | 11 +- triqs/utility/tuple_tools.hpp | 1 + 57 files changed, 1855 insertions(+), 1549 deletions(-) create mode 100644 test/speed/array3.cpp create mode 100644 test/triqs/arrays/c_ordered_transposed_view.cpp delete mode 100644 test/triqs/arrays/group_indices.output create mode 100644 triqs/arrays/group_indices.hpp delete mode 100644 triqs/arrays/impl/flags.hpp create mode 100644 triqs/arrays/impl/print.hpp delete mode 100644 triqs/arrays/indexmaps/cuboid/group_indices.hpp delete mode 100644 triqs/arrays/indexmaps/cuboid/slice_traversal_order.hpp diff --git a/doc/reference/arrays/algebras_1.cpp b/doc/reference/arrays/algebras_1.cpp index 20ed6d3b..4b1287e5 100644 --- a/doc/reference/arrays/algebras_1.cpp +++ b/doc/reference/arrays/algebras_1.cpp @@ -2,7 +2,7 @@ using namespace triqs::arrays; int main() { array A(2, 2); - matrix M(2, 2); + matrix M(2, 2); // M + A; // --> ERROR. Rejected by the compiler. M + make_matrix_view(A); //--> OK. } diff --git a/doc/reference/arrays/containers/reg_constructors_0.cpp b/doc/reference/arrays/containers/reg_constructors_0.cpp index 68b8ee11..ad7e127f 100644 --- a/doc/reference/arrays/containers/reg_constructors_0.cpp +++ b/doc/reference/arrays/containers/reg_constructors_0.cpp @@ -17,17 +17,17 @@ int main() { // a vector of double vector V(10); - // arrays with custom TraversalOrder + // arrays with custom memory layout // C-style - array A0(2, 3, 4); - array A0b; // same type but empty + array A0(2, 3, 4); + array A0b; // same type but empty // Fortran-style - array A4(2, 3, 4); - array A1b; // same type but empty + array A4(2, 3, 4, FORTRAN_LAYOUT); + array A1b(FORTRAN_LAYOUT); // same type but empty // custom : (i,j,k) : index j is fastest, then k, then i - array A2(2, 3, 4); + array A2(2, 3, 4, triqs::arrays::make_memory_layout(1, 0, 2)); } diff --git a/doc/reference/arrays/containers/regular.rst b/doc/reference/arrays/containers/regular.rst index 5ee575c9..b4a7031a 100644 --- a/doc/reference/arrays/containers/regular.rst +++ b/doc/reference/arrays/containers/regular.rst @@ -7,15 +7,9 @@ array, matrix & vector .. code-block:: c - template class array; - template class matrix; - template class vector; - -where triqs::ull_t is the type defined by : - -.. code-block:: c - - typedef unsigned long long ull_t; + template class array; + template class matrix; + template class vector; * The library provides three basic containers: @@ -41,9 +35,8 @@ Template parameters +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ | Rank | int | rank | The rank of the array | +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ -| :ref:`OptionsFlags` | ull_t | opt_flags | Compile time options | -+-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ -| :ref:`TraversalOrder` | ull_t | | Traversal Order for all loops | +| :ref:`TraversalOrder` | ull_t | | Traversal Order for loops and | +| | | | iterators | +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ NB: Rank is only present for array, since matrix have rank 2 and vector rank 1. diff --git a/doc/reference/arrays/containers/resize_0.cpp b/doc/reference/arrays/containers/resize_0.cpp index 8f9f0703..49ee13fd 100644 --- a/doc/reference/arrays/containers/resize_0.cpp +++ b/doc/reference/arrays/containers/resize_0.cpp @@ -4,7 +4,7 @@ int main() { array A(2, 3); A.resize(make_shape(5, 5)); - matrix M; + matrix M; M.resize(3, 3); vector V; diff --git a/doc/reference/arrays/containers/stream.rst b/doc/reference/arrays/containers/stream.rst index 2f6d3d09..c2b03423 100644 --- a/doc/reference/arrays/containers/stream.rst +++ b/doc/reference/arrays/containers/stream.rst @@ -7,14 +7,14 @@ operator << **Synopsis**:: - template - std::ostream & operator << (std::ostream & out, array_view const &); + template + std::ostream & operator << (std::ostream & out, array_view const &); - template - std::ostream & operator << (std::ostream & out, matrix_view const &); + template + std::ostream & operator << (std::ostream & out, matrix_view const &); - template - std::ostream & operator << (std::ostream & out, vector_view const &); + template + std::ostream & operator << (std::ostream & out, vector_view const &); Prints the view (or the container since the view is trivially build from the container) to the stream out. diff --git a/doc/reference/arrays/containers/template_parameters.rst b/doc/reference/arrays/containers/template_parameters.rst index 44c2cb69..3dca6209 100644 --- a/doc/reference/arrays/containers/template_parameters.rst +++ b/doc/reference/arrays/containers/template_parameters.rst @@ -1,51 +1,41 @@ Template parameters of the containers and views ====================================================== -.. _arr_templ_par_opt: - -OptionFlags ----------------------------- - -* OptionFlags is a series of flags determining various options at compile time. - The possible flags are accessible via a constexpr ull_t in triqs::arrays or a macro: - - ======================== ======================================= - Macro constexpr equivalent - ======================== ======================================= - BOUND_CHECK triqs::arrays::BoundCheck - TRAVERSAL_ORDER_C triqs::arrays::TraversalOrderC - TRAVERSAL_ORDER_FORTRAN triqs::arrays::TraversalOrderFortran - DEFAULT_INIT triqs::arrays::DefaultInit - ======================== ======================================= - - - Defaults can be modified with the macros: - - * `TRIQS_ARRAYS_ENFORCE_BOUNDCHECK` : enforce BoundCheck by default (slows down codes ! Use only for debugging purposes). - .. _arr_templ_par_to: TraversalOrder ---------------------------- -* TraversalOrder is a coding of the optimal ordering of indices, given by a permutation - evaluated at **compile time**. - The traversal of the arrays (iterators, foreach loop) will be written and optimised for this +* TraversalOrder is a type encoding of the optimal traversal in memory. + Default is void, corresponding to as regular C-style ordering (slowest index first). + +* The traversal of the arrays (iterators, foreach loop) will be written and optimised for this order. - The default (0) is understood as regular C-style ordering (slowest index first). +* It is indeed sometimes necessary to know *at compile time* the traversal order to generate + an optimal code. The default (C-style) is sufficient is most cases, but not always... - Note that an array can use any index ordering in memory, and that decision is take at run time - (this is necessary for interfacing with python numpy arrays, see below). - The code will be correct for any order, but optimised for the TraversalOrder. +* More explanations here... + +* Note that this notion is completely independent of the real memory layout of the array, + which is a runtime parameter. + + The code will be correct for any order, but may be faster for the TraversalOrder. - For a few very specials operations (e.g. regrouping of indices), the indices ordering in memory and TraversalOrder - must coincide. This will be explicitely said below. By default, it is not necessary (but can be suboptimal). +* TraversalOrder is not present for vector since there is only one possibility in 1d. - The permutation P encodes the position of the index: P[0] is the fastest index, P[rank - 1] the slowest index - (see examples below). - - TraversalOrder is not present for vector since there is only one possibility in 1d. +* TODO : document the various possibility beyond C style. ++--------------------------+------------------------------------------------------------+ +| TraversalOrder | Meaning | ++==========================+============================================================+ +| void | C style traversal (lowest index first). | ++--------------------------+------------------------------------------------------------+ +| _traversal_fortran | Fortran style traversal (last index first) | ++--------------------------+------------------------------------------------------------+ +| _traversal_dynamical | Traverse in the order specified by the memory layout | ++--------------------------+------------------------------------------------------------+ +| _traversal_custom<2,1,0> | Traverse in the order specified by the permutation (2,1,0) | ++--------------------------+------------------------------------------------------------+ diff --git a/doc/reference/arrays/containers/view_constructors.rst b/doc/reference/arrays/containers/view_constructors.rst index 6b1f73bc..7197509f 100644 --- a/doc/reference/arrays/containers/view_constructors.rst +++ b/doc/reference/arrays/containers/view_constructors.rst @@ -34,7 +34,7 @@ Constructors of array_const_view --------------------------------------- +------------------------------------------------------------------------+---------------------------------------------------------------------------------------+ -| Constructors of array_const_view | Comments | +| Constructors of array_const_view | Comments | +========================================================================+=======================================================================================+ | array_const_view(array_const_view const &) | Create a new view, viewing the same data. Does **not** copy data. (copy constructor) | +------------------------------------------------------------------------+---------------------------------------------------------------------------------------+ diff --git a/doc/reference/arrays/containers/views.rst b/doc/reference/arrays/containers/views.rst index ed5dbbc2..f8e8d0fc 100644 --- a/doc/reference/arrays/containers/views.rst +++ b/doc/reference/arrays/containers/views.rst @@ -9,20 +9,13 @@ Views .. code-block:: c - template class array_view; - template class matrix_view; - template class vector_view; + template class array_view; + template class matrix_view; + template class vector_view; - template class array_const_view; - template class matrix_const_view; - template class vector_const_view; - - -where triqs::ull_t is the type defined by : - -.. code-block:: c - - typedef unsigned long long ull_t; + template class array_const_view; + template class matrix_const_view; + template class vector_const_view; * The view types of X (= array, matrix, vector) are called X_view, and X_const_view, with the same template parameters as the regular type. @@ -65,8 +58,6 @@ Template parameters +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ | Rank | int | rank | The rank of the array | +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ -| :ref:`OptionsFlags` | ull_t | opt_flags | Compile time options | -+-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ | :ref:`TraversalOrder` | ull_t | | Traversal Order for all loops | +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ diff --git a/doc/reference/arrays/foreach.rst b/doc/reference/arrays/foreach.rst index f1558dae..3714fcde 100644 --- a/doc/reference/arrays/foreach.rst +++ b/doc/reference/arrays/foreach.rst @@ -16,9 +16,8 @@ several advantages: * it is more compact, less error prone (one does not need to specify the bounds in the loop). -* the library orders the for loop in the most efficient way to traverse the memory. - Hence if one change the TraversalOrder for some arrays, there is no need to change the for loop - ordering to obtain a better performance, just recompile... +* The library orders the loops in way specified by the template parameters TraversalOrder (by default, the standard C order), + to obtain the most efficient way to traverse the memory. * it is easier to write generic code for array of several dimensions. @@ -26,7 +25,7 @@ foreach ------------ The *foreach* function call a given function *f* successively on the indices of an array *A*, -in the order specified by the TraversalOrder flag of the array. +in the order specified by the TraversalOrder of the array. * Synopsis:: @@ -44,7 +43,7 @@ in the order specified by the TraversalOrder flag of the array. for (i,j,k...) F(i,j,k...) * The for loop are automatically organised to optimize the traversal order of A - using the TraversalOrder flag of the array. + using the TraversalOrder of the array. * As a result this is always equally or more optimized than a manually written loop. @@ -80,7 +79,7 @@ Synopsis:: for (i,j,k...) A(i,j,k...) = F(i,j,k...) * The for loop are automatically organised to optimize the traversal order of A - using the TraversalOrder flag of the array. + using the TraversalOrder of the array. .. triqs_example:: ./foreach_1.cpp .. note:: diff --git a/test/pytriqs/arrays/f.hpp b/test/pytriqs/arrays/f.hpp index 99b3ecfd..1d8a710c 100644 --- a/test/pytriqs/arrays/f.hpp +++ b/test/pytriqs/arrays/f.hpp @@ -6,7 +6,7 @@ void f(array_view A) { std::cout << " A(range(),range(),0) = "<< A(range(),range(),0) <(l); + speed_tester(l); + speed_tester(l); + speed_tester(l); + speed_tester(l); + speed_tester(l); + speed_tester(l); + //speed_tester(l); + return 0; +} + diff --git a/test/speed/array_lazy_assign.cpp b/test/speed/array_lazy_assign.cpp index 22c7a5bc..d076a2d9 100644 --- a/test/speed/array_lazy_assign.cpp +++ b/test/speed/array_lazy_assign.cpp @@ -41,7 +41,7 @@ struct plain { struct lazy { void operator()() { tql::placeholder<0> i_; tql::placeholder<1> j_; - triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); + triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); for (int u =0; u<5000; ++u) A(i_,j_) << i_+ 2*j_; } @@ -49,7 +49,7 @@ struct lazy { struct foreach_lambda { void operator()() { - triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); + triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); for (int u =0; u<5000; ++u) //assign_foreach(A, [](size_t i, size_t j) { return i + 2.0*j;}); assign_foreach(A, [](long i, long j) { return i + 2.0*j;}); diff --git a/test/speed/array_short1.cpp b/test/speed/array_short1.cpp index c59461b6..6a92af7d 100644 --- a/test/speed/array_short1.cpp +++ b/test/speed/array_short1.cpp @@ -37,7 +37,7 @@ struct plain_for_no_ptr_const { struct assign_to_const { void operator()() { - triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); + triqs::arrays::array A (N1,N2,FORTRAN_LAYOUT); auto V = make_view(A); for (int u =0; u<2500000; ++u) //make_view(A) = 1867; @@ -58,8 +58,8 @@ struct plain_for_no_ptr { //typedef double value_type; //typedef triqs::arrays::matrix::indexmap_type::domain_type::index_value_type index_value_type; struct F { - triqs::arrays::matrix & a; - F(triqs::arrays::matrix & a_): a(a_){} + triqs::arrays::matrix & a; + F(triqs::arrays::matrix & a_): a(a_){} //void operator()(value_type & p, index_value_type const & key) const { p=key[0]*10 + key[1] ;} //void operator()(value_type & p, value_type const & x0, value_type const & x1) const { p=x0*10 + x1 ;} template void operator()( X0 const & x0, X1 const & x1) const { a(x0,x1) =x0*10 + x1 ;} @@ -67,7 +67,7 @@ struct F { struct with_foreach { void operator()() { - triqs::arrays::matrix A (N1,N2,FORTRAN_LAYOUT); + triqs::arrays::matrix A (N1,N2,FORTRAN_LAYOUT); // triqs::arrays::indexmaps:: for (int u =0; u<2500000; ++u) foreach(A,F(A)); // for (int u =0; u<5000; ++u) make_view(A) = 1867; diff --git a/test/speed/gf_s.cpp b/test/speed/gf_s.cpp index 7915d14d..36386237 100644 --- a/test/speed/gf_s.cpp +++ b/test/speed/gf_s.cpp @@ -1,4 +1,4 @@ -#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK //#define TRIQS_ARRAYS_CHECK_WEAK_REFS diff --git a/test/triqs/arrays/bound_check.cpp b/test/triqs/arrays/bound_check.cpp index 4c478226..7441fbaa 100644 --- a/test/triqs/arrays/bound_check.cpp +++ b/test/triqs/arrays/bound_check.cpp @@ -31,10 +31,10 @@ int main(int argc, char **argv) { - array A (2,3); + array A (2,3); array B (2,3); array C(2); - array Af (2,3, FORTRAN_LAYOUT); + array Af (2,3, FORTRAN_LAYOUT); try { diff --git a/test/triqs/arrays/c_ordered_transposed_view.cpp b/test/triqs/arrays/c_ordered_transposed_view.cpp new file mode 100644 index 00000000..640c9693 --- /dev/null +++ b/test/triqs/arrays/c_ordered_transposed_view.cpp @@ -0,0 +1,36 @@ +#include "./common.hpp" + +using namespace triqs::arrays; +using namespace triqs::clef; + +template void test(T... x) { + placeholder<0> i_; + placeholder<1> j_; + placeholder<2> k_; + placeholder<3> l_; + array a(1, 2, 3, 4, make_memory_layout(x...)); + a(i_, j_, k_, l_) << i_ + 10 * j_ + 100 * k_ + 1000 * l_; + TEST_ERR(a); + + auto v = c_ordered_transposed_view(a()); + + TEST(a.shape()); + TEST(a.indexmap().get_memory_layout()); + TEST(v.shape()); + TEST(v.indexmap().get_memory_layout()); + std::cerr << "----------------" << std::endl; +} + +int main() { + + try { + + test(0, 1, 2, 3); + test(1, 0, 2, 3); + test(2, 0, 3, 1); + } + catch (std::exception const& e) { + std::cerr << e.what() << std::endl; + } +} + diff --git a/test/triqs/arrays/common.hpp b/test/triqs/arrays/common.hpp index ac1297a1..fe7696fd 100644 --- a/test/triqs/arrays/common.hpp +++ b/test/triqs/arrays/common.hpp @@ -27,6 +27,7 @@ #include #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) < "<< (X) < void f(array_view const &a) { +template void f(array_view const &a) { std::cout << a << std::endl; } diff --git a/test/triqs/arrays/ellipsis.cpp b/test/triqs/arrays/ellipsis.cpp index 60251518..938dbc23 100644 --- a/test/triqs/arrays/ellipsis.cpp +++ b/test/triqs/arrays/ellipsis.cpp @@ -21,6 +21,7 @@ #include "./common.hpp" #include +#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK #include #include #include diff --git a/test/triqs/arrays/fill_with_lazy.cpp b/test/triqs/arrays/fill_with_lazy.cpp index ec475826..38e1479d 100644 --- a/test/triqs/arrays/fill_with_lazy.cpp +++ b/test/triqs/arrays/fill_with_lazy.cpp @@ -31,7 +31,7 @@ int main(int argc, char **argv) { A(i_,j_) << i_ + 2*j_ ; std::cout << "A = "< B(2,2); + tqa::array B(2,2); B(i_,j_) << i_ + 2*j_ ; std::cout << "B = "< -#include -#include +#include -namespace tqa=triqs::arrays; -namespace tql=triqs::clef; +using namespace triqs::arrays; +using namespace triqs::clef; -using tqa::m_index; +int main() { -template void test() { + try { + placeholder<0> i_; + placeholder<1> j_; + placeholder<2> k_; + placeholder<3> l_; - tql::placeholder<0> i_; - tql::placeholder<1> j_; - tql::placeholder<2> k_; - tql::placeholder<3> l_; + { // a simple test + array A(2, 2, 2, 2); + A(i_, j_, k_, l_) << i_ + 10 * j_ + 100 * k_ + 1000 * l_; + TEST_ERR(A); + TEST_ERR(group_indices_view(A, {0,1}, {2,3})); + } - { // a simple test - tqa::array A(2,2,2,2); - A(i_,j_,k_,l_) << i_ + 10*j_ + 100*k_ + 1000*l_; - TEST(A); - TEST( group_indices(A, m_index<0,1>(), m_index<2,3>())); + // more complex : inversing a tensor product of matrices... + matrix B(2, 2), C(3, 3), Binv, Cinv; + C(i_, j_) << 1.7 / (3.4 * i_ - 2.3 * j_ + 1); + B(i_, j_) << 2 * i_ + j_; + TEST_ERR(B); + TEST_ERR(C); + Binv = inverse(B); + Cinv = inverse(C); + TEST_ERR(Binv); + TEST_ERR(Cinv); + + std::cerr << " --- (1) ---"<< std::endl; + { + array A(2, 3, 2, 3); + A(i_, j_, k_, l_) << B(i_, k_) * C(j_, l_); + + TEST_ERR(A); + + auto M = make_matrix_view(group_indices_view(A, {0, 1}, {2, 3})); + M = inverse(M); + + // checking + array R(A.shape()); + R(i_, j_, k_, l_) << Binv(i_, k_) * Cinv(j_, l_); + + TEST_ERR(R); + TEST_ERR(A); + assert_all_close(R, A, 5.e-15); + } + + std::cerr << " --- (2) ---"<< std::endl; + { + array A(2, 3, 2, 3, make_memory_layout(1, 0, 3, 2)); + A(i_, j_, k_, l_) << B(i_, k_) * C(j_, l_); + + TEST_ERR(A); + + auto M = make_matrix_view(group_indices_view(A, {0, 1}, {2, 3})); + M = inverse(M); + + // checking + array R(A.shape()); + R(i_, j_, k_, l_) << Binv(i_, k_) * Cinv(j_, l_); + + TEST_ERR(R); + TEST_ERR(A); + assert_all_close(R, A, 5.e-15); + } + + std::cerr << " --- (3) ---"<< std::endl; + { + array A(2, 2, 3, 3, make_memory_layout(0, 2, 1, 3)); + A(i_, k_, j_, l_) << B(i_, k_) * C(j_, l_); + + TEST_ERR(A); + + auto M = make_matrix_view(group_indices_view(A, {0, 2}, {1, 3})); + M = inverse(M); + + // checking + array R(A.shape()); + R(i_, k_, j_, l_) << Binv(i_, k_) * Cinv(j_, l_); + + TEST_ERR(R); + TEST_ERR(A); + assert_all_close(R, A, 5.e-15); + } } - - { // more complex : inversing a tensor product of matrices... - tqa::matrix B(2,2), C(3,3), Binv, Cinv; - C(i_,j_) << 1.7 /( 3.4*i_ - 2.3*j_ + 1) ; - B(i_,j_) << 2*i_ + j_; - TEST(B); TEST(C); - Binv = inverse(B); - Cinv = inverse(C); - TEST(Binv); TEST(Cinv); - - tqa::array A(2,3,2,3); - A(i_,j_,k_,l_) << B(i_,k_) * C(j_,l_); - - TEST(A); - - tqa::matrix_view M = group_indices (A, m_index<0,1>() , m_index<2,3>() ); - M = inverse(M); - - // checking - tqa::array R(A.shape()); - R(i_,j_,k_,l_) << Binv(i_,k_) * Cinv(j_,l_); - - TEST(R); - TEST(A); - assert_all_close(R,A,5.e-15); - //TEST(make_array(R-A)); - //TEST( max(abs(R-A))); + catch (std::exception const& e) { + std::cerr << e.what() << std::endl; } - } -int main () { - test(); - test(); - - // template dans le namespace clef pour ADL et enable sur le concept de A ... - // an alternative syntax ? - //auto V = group_indices( A(i_,j_,k_,l_), var(i_,j_), var(k_,l_)); - // match A(...), deduce the int of hte ph, and make a finder for a ph in it - // finder::find::value --> m_index<...> -} diff --git a/test/triqs/arrays/group_indices.output b/test/triqs/arrays/group_indices.output deleted file mode 100644 index c1a989e1..00000000 --- a/test/triqs/arrays/group_indices.output +++ /dev/null @@ -1,46 +0,0 @@ -(A) ---> [0 1 10 11 100 101 110 111 1000 1001 1010 1011 1100 1101 1110 1111 ] -(group_indices(A, m_index<0,1>(), m_index<2,3>())) ---> -[[0,100,1000,1100] - [1,101,1001,1101] - [10,110,1010,1110] - [11,111,1011,1111]] -(B) ---> -[[0,1] - [2,3]] -(C) ---> -[[1.7,-1.30769,-0.472222] - [0.386364,0.809524,-8.5] - [0.217949,0.309091,0.53125]] -(Binv) ---> -[[-1.5,0.5] - [1,0]] -(Cinv) ---> -[[0.386252,0.0693273,1.45257] - [-0.259977,0.1271,1.80251] - [-0.00720283,-0.102391,0.237694]] -(A) ---> [0 3.4 0 0.772727 0 0.435897 1.7 5.1 0.386364 1.15909 0.217949 0.653846 0 -2.61538 0 1.61905 0 0.618182 -1.30769 -3.92308 0.809524 2.42857 0.309091 0.927273 0 -0.944444 0 -17 0 1.0625 -0.472222 -1.41667 -8.5 -25.5 0.53125 1.59375 ] -(R) ---> [-0.579378 0.386252 0.389966 -0.259977 0.0108043 -0.00720283 0.193126 0 -0.129989 0 -0.00360142 0 -0.103991 0.0693273 -0.19065 0.1271 0.153587 -0.102391 0.0346636 0 0.06355 0 -0.0511955 0 -2.17886 1.45257 -2.70376 1.80251 -0.35654 0.237694 0.726286 0 0.901255 0 0.118847 0 ] -(A) ---> [-0.579378 0.386252 0.389966 -0.259977 0.0108043 -0.00720283 0.193126 0 -0.129989 0 -0.00360142 0 -0.103991 0.0693273 -0.19065 0.1271 0.153587 -0.102391 0.0346636 0 0.06355 0 -0.0511955 0 -2.17886 1.45257 -2.70376 1.80251 -0.35654 0.237694 0.726286 0 0.901255 0 0.118847 0 ] -(A) ---> [0 1 10 11 100 101 110 111 1000 1001 1010 1011 1100 1101 1110 1111 ] -(group_indices(A, m_index<0,1>(), m_index<2,3>())) ---> -[[0,1000,100,1100] - [10,1010,110,1110] - [1,1001,101,1101] - [11,1011,111,1111]] -(B) ---> -[[0,1] - [2,3]] -(C) ---> -[[1.7,-1.30769,-0.472222] - [0.386364,0.809524,-8.5] - [0.217949,0.309091,0.53125]] -(Binv) ---> -[[-1.5,0.5] - [1,0]] -(Cinv) ---> -[[0.386252,0.0693273,1.45257] - [-0.259977,0.1271,1.80251] - [-0.00720283,-0.102391,0.237694]] -(A) ---> [0 3.4 0 0.772727 0 0.435897 1.7 5.1 0.386364 1.15909 0.217949 0.653846 0 -2.61538 0 1.61905 0 0.618182 -1.30769 -3.92308 0.809524 2.42857 0.309091 0.927273 0 -0.944444 0 -17 0 1.0625 -0.472222 -1.41667 -8.5 -25.5 0.53125 1.59375 ] -(R) ---> [-0.579378 0.386252 0.389966 -0.259977 0.0108043 -0.00720283 0.193126 0 -0.129989 0 -0.00360142 0 -0.103991 0.0693273 -0.19065 0.1271 0.153587 -0.102391 0.0346636 0 0.06355 0 -0.0511955 0 -2.17886 1.45257 -2.70376 1.80251 -0.35654 0.237694 0.726286 0 0.901255 0 0.118847 0 ] -(A) ---> [-0.579378 0.386252 0.389966 -0.259977 0.0108043 -0.00720283 0.193126 0 -0.129989 0 -0.00360142 0 -0.103991 0.0693273 -0.19065 0.1271 0.153587 -0.102391 0.0346636 0 0.06355 0 -0.0511955 0 -2.17886 1.45257 -2.70376 1.80251 -0.35654 0.237694 0.726286 0 0.901255 0 0.118847 0 ] diff --git a/test/triqs/arrays/iterators.cpp b/test/triqs/arrays/iterators.cpp index a61a3b7b..c5fa4d19 100644 --- a/test/triqs/arrays/iterators.cpp +++ b/test/triqs/arrays/iterators.cpp @@ -65,7 +65,7 @@ int main(int argc, char **argv) { std::cout<<" F order : traversal"< Af (2,3); + array Af (2,3, FORTRAN_LAYOUT); for (auto it = Af.begin(); it; ++it) { *it =it.indices()[0] + 10 *it.indices()[1] ; diff --git a/test/triqs/arrays/reinterpret_view.cpp b/test/triqs/arrays/reinterpret_view.cpp index 545ea067..509c4f4b 100644 --- a/test/triqs/arrays/reinterpret_view.cpp +++ b/test/triqs/arrays/reinterpret_view.cpp @@ -22,40 +22,20 @@ #include #include -namespace triqs { namespace arrays { +namespace triqs { namespace arrays { -/* template - array_view reinterpret_array_view (array const & a, I ... index) { - //static int constexpr rank = sizeof...(I); - typedef array_view return_t; - //typedef array_view return_t; - return return_t { make_shape(index...), a.storage() }; - //return return_t (typename return_t::indexmap_type (mini_vector(index...)) , a.storage()); - } -*/ - template - array_view reinterpret (array const & a, I ... index) { + template + array_view reinterpret(array const &a, I... index) { return { {make_shape(index...)}, a.storage() }; } // wrong for views - template - array_view reinterpret_array_view (array_view const & a, I ... index) { + template + array_view reinterpret_array_view (array_view const & a, I ... index) { if (!has_contiguous_data(a)) TRIQS_RUNTIME_ERROR << "reinterpretation failure : data of the view are not contiguous"; return { {make_shape(index...)}, a.storage() }; } -/* - template - array_view reinterpret_array_view_c (matrix const & a, I ... index) { - return { {make_shape(index...)}, a.storage() }; - } - -template - array_view reinterpret_array_view_c (vector const & a, I ... index) { - return { {make_shape(index...)}, a.storage() }; - } -*/ }} diff --git a/test/triqs/arrays/slice_index_order.cpp b/test/triqs/arrays/slice_index_order.cpp index f54acbbc..f90ba4c8 100644 --- a/test/triqs/arrays/slice_index_order.cpp +++ b/test/triqs/arrays/slice_index_order.cpp @@ -20,96 +20,81 @@ ******************************************************************************/ #include "./common.hpp" #include -#include -#include -#include -#include +#include using namespace triqs::arrays; -using namespace triqs::arrays::permutations; +using namespace triqs::arrays::indexmaps; -template struct P { - static constexpr ull value = F; - friend std::ostream & operator <<( std::ostream & out, P const & s) { - //out << "Permutation of size " << permutations::size(s.value) << " : "<< std::hex; - out << std::hex; - s.print(out, std::integral_constant()); - return out << std::dec; - } +template void test(ML1 ml1, ML2 ml2, Args... x) { - template void print( std::ostream & out, std::integral_constant) const { out << apply(this->value,c); print(out, std::integral_constant());} - void print( std::ostream & out, std::integral_constant) const {} -}; + cuboid::map m(ml2); + slicer, Args...> S; + auto m2 = S.invoke(m, x...); + std::cout << m.get_memory_layout() << std::endl; + std::cout << m2.get_memory_layout() << std::endl; -template void test() { - - typedef indexmaps::cuboid::slicing_TO_order::sliced_memory_order S1; - constexpr auto s1 = S1::value; - //std::cout << " sliced "<< std::hex<< S1::value << " " <() << std::endl; - std::cout << P() << std::endl; - //std::cout << "mashk "<< std::hex<(); - test< permutation(0,1,2) ,permutation(0,1,2,3), range, range,int, range>(); - test< permutations::identity(4) ,permutation(0,1,2,3), range, range, range, range>(); - //test< permutation(0,1,2,3) ,permutation(0,1,2,3), range, range, range, range>(); - test< 0 ,permutation(0,1,2,3), int, int, int, int> (); + std::cout << " F order " << std::endl; + test(make_memory_layout(1, 0), make_memory_layout(3, 2, 1, 0), 0, range(), 0, range()); + test(make_memory_layout(2, 1, 0), make_memory_layout(3, 2, 1, 0), 0, range(), range(), range()); + test(make_memory_layout(2, 1, 0), make_memory_layout(3, 2, 1, 0), range(), 0, range(), range()); + test(make_memory_layout(2, 1, 0), make_memory_layout(3, 2, 1, 0), range(), range(), 0, range()); + test(make_memory_layout(2, 1, 0), make_memory_layout(3, 2, 1, 0), range(), range(), range(), 0); + test(make_memory_layout(3, 2, 1, 0), make_memory_layout(3, 2, 1, 0), range(), range(), range(), range()); + test(make_memory_layout(), make_memory_layout(3, 2, 1, 0), 0, 0, 0, 0); - std::cout << " c order " << std::endl ; - test< permutation(1,0) ,permutation(3,2,1,0), int, range,int, range>(); - test< permutation(2,1,0) ,permutation(3,2,1,0), int, range, range, range>(); - test< permutation(2,1,0) ,permutation(3,2,1,0), range,int, range, range>(); - test< permutation(2,1,0) ,permutation(3,2,1,0), range, range,int, range>(); - test< permutation(2,1,0) ,permutation(3,2,1,0), range, range, range, int>(); - test< permutation(3,2,1,0) ,permutation(3,2,1,0), range, range, range, range>(); - test< 0 ,permutation(3,2,1,0), int, int, int, int> (); + test(make_memory_layout(0), make_memory_layout(0, 1), 0, range()); - test< permutation(0), permutation(0,1), int, range>(); + std::cout << " custom order " << std::endl; + test(make_memory_layout(1, 0), make_memory_layout(0, 3, 1, 2), 0, range(), 0, range()); + test(make_memory_layout(2, 0, 1), make_memory_layout(0, 3, 1, 2), 0, range(), range(), range()); + test(make_memory_layout(0, 2, 1), make_memory_layout(0, 3, 1, 2), range(), 0, range(), range()); + test(make_memory_layout(0, 2, 1), make_memory_layout(0, 3, 1, 2), range(), range(), 0, range()); + test(make_memory_layout(0, 1, 2), make_memory_layout(0, 3, 1, 2), range(), range(), range(), 0); + test(make_memory_layout(0, 3, 1, 2), make_memory_layout(0, 3, 1, 2), range(), range(), range(), range()); + test(make_memory_layout(), make_memory_layout(0, 3, 1, 2), 0, 0, 0, 0); - std::cout << " custom order " << std::endl ; - test< permutation(1,0) ,permutation(0,3,1,2), int, range,int, range>(); - test< permutation(2,0,1) ,permutation(0,3,1,2), int, range, range, range>(); - test< permutation(0,2,1) ,permutation(0,3,1,2), range,int, range, range>(); - test< permutation(0,2,1) ,permutation(0,3,1,2), range, range,int, range>(); - test< permutation(0,1,2) ,permutation(0,3,1,2), range, range, range, int>(); - test< permutation(0,3,1,2) ,permutation(0,3,1,2), range, range, range, range>(); - test< 0 ,permutation(0,3,1,2), int, int, int, int> (); + std::cout << " ----------- custom order ------------- " << std::endl; - std::cout << " ----------- custom order ------------- " << std::endl ; + std::cout << " ---- 0 0 " << std::endl; + test(make_memory_layout(2, 0, 3, 1), make_memory_layout(2, 0, 3, 1), range(), range(), range(), range()); - std::cout << " ---- 0 int "<< std::endl ; - test< permutation(2,0,3,1) ,permutation(2,0,3,1), range, range, range, range>(); + std::cout << " ---- 1 0 " << std::endl; + test(make_memory_layout(1, 2, 0), make_memory_layout(2, 0, 3, 1), 0, range(), range(), range()); + test(make_memory_layout(1, 0, 2), make_memory_layout(2, 0, 3, 1), range(), 0, range(), range()); + test(make_memory_layout(0, 2, 1), make_memory_layout(2, 0, 3, 1), range(), range(), 0, range()); + test(make_memory_layout(2, 0, 1), make_memory_layout(2, 0, 3, 1), range(), range(), range(), 0); - std::cout << " ---- 1 int "<< std::endl ; - test< permutation(1,2,0) ,permutation(2,0,3,1), int, range, range, range>(); - test< permutation(1,0,2) ,permutation(2,0,3,1), range,int, range, range>(); - test< permutation(0,2,1) ,permutation(2,0,3,1), range, range,int, range>(); - test< permutation(2,0,1) ,permutation(2,0,3,1), range, range, range, int>(); + std::cout << " ---- 2 0 " << std::endl; + test(make_memory_layout(0, 1), make_memory_layout(2, 0, 3, 1), 0, 0, range(), range()); + test(make_memory_layout(1, 0), make_memory_layout(2, 0, 3, 1), 0, range(), 0, range()); + test(make_memory_layout(1, 0), make_memory_layout(2, 0, 3, 1), 0, range(), range(), 0); + test(make_memory_layout(0, 1), make_memory_layout(2, 0, 3, 1), range(), 0, 0, range()); + test(make_memory_layout(1, 0), make_memory_layout(2, 0, 3, 1), range(), 0, range(), 0); + test(make_memory_layout(0, 1), make_memory_layout(2, 0, 3, 1), range(), range(), 0, 0); - std::cout << " ---- 2 int "<< std::endl ; - test< permutation(0,1) ,permutation(2,0,3,1), int, int, range, range>(); - test< permutation(1,0) ,permutation(2,0,3,1), int, range,int, range>(); - test< permutation(1,0) ,permutation(2,0,3,1), int, range,range, int>(); - test< permutation(0,1) ,permutation(2,0,3,1), range,int, int, range>(); - test< permutation(1,0) ,permutation(2,0,3,1), range,int, range, int>(); - test< permutation(0,1) ,permutation(2,0,3,1), range, range,int, int>(); - - std::cout << " ---- 3 int "<< std::endl ; - test< permutation(0) ,permutation(2,0,3,1), range, int, int, int>(); - test< permutation(0) ,permutation(2,0,3,1), int,range, int, int>(); - test< permutation(0) ,permutation(2,0,3,1), int, int,range, int>(); - test< permutation(0) ,permutation(2,0,3,1), int, int, int, range>(); - - std::cout << " ---- 4 int "<< std::endl ; - test< 0 ,permutation(2,0,3,1), int, int, int, int> (); - std::cout << " OK "<< std::endl ; + std::cout << " ---- 3 0 " << std::endl; + test(make_memory_layout(0), make_memory_layout(2, 0, 3, 1), range(), 0, 0, 0); + test(make_memory_layout(0), make_memory_layout(2, 0, 3, 1), 0, range(), 0, 0); + test(make_memory_layout(0), make_memory_layout(2, 0, 3, 1), 0, 0, range(), 0); + test(make_memory_layout(0), make_memory_layout(2, 0, 3, 1), 0, 0, 0, range()); + std::cout << " ---- 4 0 " << std::endl; + test(make_memory_layout(), make_memory_layout(2, 0, 3, 1), 0, 0, 0, 0); + std::cout << " OK " << std::endl; } diff --git a/test/triqs/arrays/views3.cpp b/test/triqs/arrays/views3.cpp index e50fa07e..5bcfc908 100644 --- a/test/triqs/arrays/views3.cpp +++ b/test/triqs/arrays/views3.cpp @@ -24,67 +24,61 @@ #include #include -using std::cout; using std::endl; +using std::cout; +using std::endl; using namespace triqs::arrays; /// -template -void print_in_indexmap_order (std::ostream & out, IM const & im, A const & a) { - out<<"["; - //for (typename IM::iterator it(im); it; ++it) out< void print_in_indexmap_order(std::ostream &out, IM const &im, A const &a) { + out << "["; for (typename IM::iterator it(im); it; ++it) { auto i = it.indices(); - out<< a(i[0],i[1],i[2])<<" "; + out << a(i[0], i[1], i[2]) << " "; } - out<<"]"; + out << "]"; } -template -void f(AA & A) { +template void f(AA &A) { - for (int i =0; i<2; ++i) - for (int j=0; j<3; ++j) - for (int k=0; k<4; ++k) - A(i,j,k) = 100*(i+1)+ 10*(j+1) + (k+1); + for (int i = 0; i < 2; ++i) + for (int j = 0; j < 3; ++j) + for (int k = 0; k < 4; ++k) A(i, j, k) = 100 * (i + 1) + 10 * (j + 1) + (k + 1); - print_in_indexmap_order(std::cout,A.indexmap(),A);std::cout< > > A0 (2,3,4); - //array > > A1 (2,3,4); - //array > > A2 (2,3,4); - // permutation in triqs::arrays - array A0 (2,3,4); - array A1 (2,3,4); - array A2 (2,3,4); - array A3 (2,3,4); - array A4 (2,3,4,FORTRAN_LAYOUT); + // permutation in triqs::arrays + array A0(2, 3, 4); + array A1(2, 3, 4, make_memory_layout(2, 1, 0)); + array A2(2, 3, 4, make_memory_layout(2, 0, 1)); + array A3(2, 3, 4, make_memory_layout(0, 1, 2)); + array A4(2, 3, 4, FORTRAN_LAYOUT); - f(A0); - f(A1); - f(A2); - f(A3); - f(A4); + f(A0); + f(A1); + f(A2); + f(A3); + f(A4); } - catch (const char * ERR) { + catch (const char *ERR) { - std::cout<<"Error is "<(t); std::cout << m<< std::endl ; } - +*/ { // filter std::cout << " ----- filter ----"<< std::endl ; auto t= std::make_tuple(0,1,2,3,4,"=5"); diff --git a/test/triqs/utility/view_tools.cpp b/test/triqs/utility/view_tools.cpp index bb7749d3..f4170950 100644 --- a/test/triqs/utility/view_tools.cpp +++ b/test/triqs/utility/view_tools.cpp @@ -40,8 +40,8 @@ int main(int argc, char **argv) { check> (); check>(); - check>(); - check>(); + check>(); + check>(); return 0; } diff --git a/triqs/arrays.hpp b/triqs/arrays.hpp index 54c2033c..e77500cd 100644 --- a/triqs/arrays.hpp +++ b/triqs/arrays.hpp @@ -21,6 +21,8 @@ #ifndef TRIQS_ARRAYS_ALL_H #define TRIQS_ARRAYS_ALL_H +#define TRIQS_ARRAYS_INCLUDED + // The basic classes #include #include @@ -36,6 +38,9 @@ // HDF5 interface #include +// Regrouping indices +#include + // Linear algebra ?? Keep here ? //#include diff --git a/triqs/arrays/array.hpp b/triqs/arrays/array.hpp index e108fb57..a5e6957b 100644 --- a/triqs/arrays/array.hpp +++ b/triqs/arrays/array.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011-2013 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,44 +18,44 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_ARRAY_H -#define TRIQS_ARRAYS_ARRAY_H +#pragma once #include #include "indexmaps/cuboid/map.hpp" #include "indexmaps/cuboid/slice.hpp" #include "impl/indexmap_storage_pair.hpp" #include "impl/assignment.hpp" -#include "impl/flags.hpp" namespace triqs { namespace arrays { - template class array_view; - template class array; + template class array_view; + template class array; // ---------------------- array_view -------------------------------- -#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map, \ - storages::shared_block, Opt, TraversalOrder, IsConst, true, Tag::array_view > +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, TraversalOrder, IsConst, true, \ + Tag::array_view> - template + template class array_view: Tag::array_view, TRIQS_CONCEPT_TAG_NAME(MutableArray), 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 regular_type; - typedef array_view view_type; - typedef array_view const_view_type; - typedef array_view weak_view_type; - + + public: + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using storage_type = typename IMPL_TYPE::storage_type; + using regular_type = array; + using view_type = array_view; + using const_view_type = array_view; + using weak_view_type = array_view; + /// Build from an IndexMap and a storage template array_view(indexmap_type const& Ind, S const& Mem) : IMPL_TYPE(Ind, Mem) {} /// Copy constructor - array_view(array_view const & X): IMPL_TYPE(X.indexmap(),X.storage()) {} + array_view(array_view const& X) : IMPL_TYPE(X.indexmap(), X.storage()) {} /// Build from anything that has an indexmap and a storage compatible with this class - template array_view(const ISP & X): IMPL_TYPE(X.indexmap(),X.storage()) { - // to be activated + template array_view(const ISP& X) : IMPL_TYPE(X.indexmap(), X.storage()) { + // to be activated static_assert(IsConst || (!ISP::is_const), "Cannot construct a non const view from a const one !"); } @@ -79,7 +79,7 @@ namespace triqs { namespace arrays { } // rebind the other view, iif this is const, and the other is not. - template ENABLE_IFC(C) rebind(array_view const& X) { + template ENABLE_IFC(C) rebind(array_view const& X) { this->indexmap_ = X.indexmap_; this->storage_ = X.storage_; } @@ -96,45 +96,52 @@ namespace triqs { namespace arrays { // to forbid serialization of views... //template void serialize(Archive & ar, const unsigned int version) = delete; - template friend array_view transposed_view(array_view const& a, INT... is) { - return array_view{transpose(a.indexmap_, is...), a.storage_}; + template friend view_type transposed_view(array_view const& a, INT... is) { + return transposed_view(a, mini_vector{is...}); + } + + friend view_type transposed_view(array_view const& a, mini_vector const& perm) { + return {transpose(a.indexmap_, perm), a.storage_}; + } + + friend view_type c_ordered_transposed_view(array_view const& a) { + return transposed_view(a, a.indexmap().get_memory_layout().get_memory_positions()); } }; - #undef IMPL_TYPE - template - using array_const_view = array_view; + template + using array_const_view = array_view; //------------------------------- array --------------------------------------------------- -#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map, \ - storages::shared_block, Opt, TraversalOrder, false, false, Tag::array_view > +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, TraversalOrder, false, false, Tag::array_view> - template - class array: Tag::array, TRIQS_CONCEPT_TAG_NAME(MutableArray), 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 regular_type; - typedef array_view view_type; - typedef array_view const_view_type; - typedef array_view weak_view_type; + template + class array : Tag::array, TRIQS_CONCEPT_TAG_NAME(MutableArray), public IMPL_TYPE { + public: + using value_type = typename IMPL_TYPE::value_type; + using storage_type = typename IMPL_TYPE::storage_type; + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using regular_type = array; + using view_type = array_view; + using const_view_type = array_view; + using weak_view_type = array_view; - /// Empty array. - explicit array(memory_layout ml = memory_layout(IMPL_TYPE::indexmap_type::traversal_order)) :IMPL_TYPE(ml){} + /// Empty array. + explicit array(memory_layout ml = memory_layout{}) : IMPL_TYPE(ml) {} - /// From a domain - explicit array( typename indexmap_type::domain_type const & dom, memory_layout ml = memory_layout(IMPL_TYPE::indexmap_type::traversal_order)): - IMPL_TYPE(indexmap_type(dom,ml)){} + /// From a domain + explicit array(typename indexmap_type::domain_type const& dom, memory_layout ml = memory_layout{}) + : IMPL_TYPE(indexmap_type(dom, ml)) {} #ifdef TRIQS_DOXYGEN /// Construction from the dimensions. NB : the number of parameters must be exactly rank (checked at compile time). array (size_t I_1, .... , size_t I_rank); #else #define IMPL(z, NN, unused) \ - explicit array (BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(NN), size_t I_),memory_layout ml = memory_layout(IMPL_TYPE::indexmap_type::traversal_order)): \ + explicit array (BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(NN), size_t I_),memory_layout ml = memory_layout{}): \ IMPL_TYPE(indexmap_type(mini_vector(BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(NN), I_)),ml)) {\ static_assert(IMPL_TYPE::rank-1==NN,"array : incorrect number of variables in constructor");} BOOST_PP_REPEAT(ARRAY_NRANK_MAX , IMPL, nil) @@ -148,18 +155,19 @@ namespace triqs { namespace arrays { explicit array(array && X) { this->swap_me(X); } // from a temporary storage and an indexmap. Used for reshaping a temporary array - explicit array(typename indexmap_type::domain_type const & dom, storage_type && sto, - memory_layout ml = memory_layout(IMPL_TYPE::indexmap_type::traversal_order)): - IMPL_TYPE(indexmap_type(dom,ml), std::move(sto)){} + explicit array(typename indexmap_type::domain_type const& dom, storage_type&& sto, memory_layout ml = memory_layout{}) + : IMPL_TYPE(indexmap_type(dom, ml), std::move(sto)) {} /** * Build a new array from X.domain() and fill it with by evaluating X. X can be : * - another type of array, array_view, matrix,.... (any pair) * - a expression : e.g. array A = B+ 2*C; */ - template - array(const T & X, TYPE_ENABLE_IF(memory_layout, ImmutableCuboidArray) ml = memory_layout(IMPL_TYPE::indexmap_type::traversal_order)): - IMPL_TYPE(indexmap_type(X.domain(),ml)) { triqs_arrays_assign_delegation(*this,X); } + template + array(const T& X, TYPE_ENABLE_IF(memory_layout, ImmutableCuboidArray) ml = memory_layout{}) + : IMPL_TYPE(indexmap_type(X.domain(), ml)) { + triqs_arrays_assign_delegation(*this, X); + } #ifdef TRIQS_WITH_PYTHON_SUPPORT ///Build from a numpy.array X (or any object from which numpy can make a numpy.array). Makes a copy. @@ -169,14 +177,14 @@ namespace triqs { namespace arrays { // build from a init_list template array(std::initializer_list const & l, typename std::enable_if::type * dummy =0): - IMPL_TYPE(indexmap_type(mini_vector(l.size()),memory_layout(IMPL_TYPE::indexmap_type::traversal_order))) { + IMPL_TYPE(indexmap_type(mini_vector(l.size()),memory_layout())) { size_t i=0; for (auto const & x : l) (*this)(i++) = x; } template array (std::initializer_list> const & l, typename std::enable_if::type * dummy =0): - IMPL_TYPE(memory_layout(IMPL_TYPE::indexmap_type::traversal_order)) { + IMPL_TYPE(memory_layout()) { size_t i=0,j=0; int s=-1; for (auto const & l1 : l) { if (s==-1) s= l1.size(); else if (s != l1.size()) TRIQS_RUNTIME_ERROR << "initializer list not rectangular !";} IMPL_TYPE::resize(typename IMPL_TYPE::domain_type (mini_vector(l.size(),s))); @@ -217,10 +225,10 @@ namespace triqs { namespace arrays { TRIQS_DEFINE_COMPOUND_OPERATORS(array); template friend const_view_type transposed_view(array const& a, INT... is) { - return view_type{transpose(a.indexmap_, is...), a.storage_}; + return transposed_view(a(), is...); }; template friend view_type transposed_view(array& a, INT... is) { - return view_type{transpose(a.indexmap_, is...), a.storage_}; + return transposed_view(a(), is...); }; };//array class @@ -230,19 +238,18 @@ namespace triqs { namespace arrays { //---------------------------------------------------------------------------------- // how to build the view type .... - template < class V, int R, ull_t Opt, ull_t TraversalOrder, bool Borrowed, bool IsConst > - struct ISPViewType< V, R, Opt, TraversalOrder,Tag::array_view, Borrowed,IsConst> { typedef array_view type; }; - + template + struct ISPViewType { + using type = array_view; + }; }}//namespace triqs::arrays // The std::swap is WRONG for a view because of the copy/move semantics of view. // Use swap instead (the correct one, found by ADL). -namespace std { - template - void swap( triqs::arrays::array_view & a , triqs::arrays::array_view & b)= delete; +namespace std { +template +void swap(triqs::arrays::array_view& a, triqs::arrays::array_view& b) = delete; } #include "./expression_template/array_algebra.hpp" -#endif - diff --git a/triqs/arrays/blas_lapack/gemm.hpp b/triqs/arrays/blas_lapack/gemm.hpp index f63dfd9b..1d854a9d 100644 --- a/triqs/arrays/blas_lapack/gemm.hpp +++ b/triqs/arrays/blas_lapack/gemm.hpp @@ -132,8 +132,8 @@ namespace triqs { namespace arrays { namespace blas { // to allow gemm (alpha, a, b, beta, M(..., ...)) i.e. a temporary view, which is not matched by previos templates // which require an lvalue. This is the only version which takes an && as last argument // indeed, in the routine, c is a *lvalue*, since it has a name, and hence we call *other* overload of the function - template - void gemm (A alpha, MT1 const & a, MT2 const & b, B beta, matrix_view && c) { gemm(alpha,a,b,beta,c);} + template + void gemm (A alpha, MT1 const & a, MT2 const & b, B beta, matrix_view && c) { gemm(alpha,a,b,beta,c);} }}}// namespace diff --git a/triqs/arrays/blas_lapack/gemv.hpp b/triqs/arrays/blas_lapack/gemv.hpp index f898efb7..d5f5f5d4 100644 --- a/triqs/arrays/blas_lapack/gemv.hpp +++ b/triqs/arrays/blas_lapack/gemv.hpp @@ -98,8 +98,8 @@ namespace triqs { namespace arrays { namespace blas { // to allow gem (alpha, a, b, beta, M(..., ...)) i.e. a temporary view, which is not matched by previos templates // which require an lvalue. - template - void gemv (A alpha, MT const & a, VT const & b, B beta, vector_view && c) { gemv(alpha,a,b,beta,c);} + template + void gemv (A alpha, MT const & a, VT const & b, B beta, vector_view && c) { gemv(alpha,a,b,beta,c);} }}}// namespace diff --git a/triqs/arrays/blas_lapack/ger.hpp b/triqs/arrays/blas_lapack/ger.hpp index 88f53b4c..af825570 100644 --- a/triqs/arrays/blas_lapack/ger.hpp +++ b/triqs/arrays/blas_lapack/ger.hpp @@ -74,8 +74,8 @@ namespace triqs { namespace arrays { namespace blas { // to allow ger (alpha, x,y, M(..., ...)) i.e. a temporary view, which is not matched by previos templates // which require an lvalue - template< typename A, typename VTX, typename VTY, typename V, ull_t Opt, ull_t To, bool W> - void ger (A alpha, VTX const & x, VTY const & y, matrix_view && r) { ger(alpha,x,y,r);} + template< typename A, typename VTX, typename VTY, typename V, typename To, bool W> + void ger (A alpha, VTX const & x, VTY const & y, matrix_view && r) { ger(alpha,x,y,r);} }}}// namespace #endif diff --git a/triqs/arrays/cache.hpp b/triqs/arrays/cache.hpp index cd0b94f3..89795bc5 100644 --- a/triqs/arrays/cache.hpp +++ b/triqs/arrays/cache.hpp @@ -64,7 +64,7 @@ namespace arrays { exposed_view_type view2() const { return view1(bool_constant::value>()); } bool need_copy_dynamic(DataType const& x, std::false_type) const { return false; } - bool need_copy_dynamic(DataType const& x, std::true_type) const { return (x.indexmap().memory_indices_layout() != ml); } + bool need_copy_dynamic(DataType const& x, std::true_type) const { return (x.indexmap().get_memory_layout() != ml); } void back_update(std::true_type) {} void back_update(std::false_type) { @@ -98,12 +98,12 @@ namespace arrays { template using const_cache = cache_impl; template - const_cache make_const_cache(D const& x, memory_layout ml = memory_layout{'C'}) { + const_cache make_const_cache(D const& x, memory_layout ml = memory_layout{}) { return {x, ml}; } template - cache make_cache(D const& x, memory_layout ml = memory_layout{'C'}) { + cache make_cache(D const& x, memory_layout ml = memory_layout{}) { return {x, ml}; } } diff --git a/triqs/arrays/group_indices.hpp b/triqs/arrays/group_indices.hpp new file mode 100644 index 00000000..ba1434de --- /dev/null +++ b/triqs/arrays/group_indices.hpp @@ -0,0 +1,77 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014 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 . + * + ******************************************************************************/ +#pragma once +#include +#include + +namespace triqs { +namespace arrays { + + /** + * Regroup indices for a C array + * Usage : group_indices_view(A, {0,1}, {2,3}) + * Precondition : + * - every indices is listed in the {...} exactly once. + * - the indices in one group are consecutive in memory. + * Routine is not extra-fast, since it does check the preconditions + */ + template + array_view group_indices_view(A const& a, std::initializer_list... args) { + + auto indices_groups = std::vector>{std::vector(args)...}; + + static const int Rf = sizeof...(U); + if (Rf != indices_groups.size()) + TRIQS_RUNTIME_ERROR << "Return rank is " << Rf << " but " << indices_groups.size() << " indices are given"; + + /// First check indices are contiguous in memory + auto mem_pos = a.indexmap().get_memory_layout().get_memory_positions(); + for (auto const& gr : indices_groups) { + auto gr2 = std::vector(gr.size()); + for (int i = 0; i < int(gr.size()); ++i) gr2[i] = mem_pos[gr[i]]; + if (*std::max_element(begin(gr2), end(gr2)) - *std::min_element(begin(gr2), end(gr2)) + 1 != gr2.size()) + TRIQS_RUNTIME_ERROR << "Indices not contiguous in memory"; + } + + /// Now compute the new lengths and strides. + mini_vector l; + mini_vector s; + auto la = a.indexmap().lengths(); + l = 1; + s = a.num_elements(); // init : superior to all strides, so a good start for the min ... + for (int i = 0; i < Rf; ++i) { + for (auto ind : indices_groups[i]) { + if (la[ind] == -1) TRIQS_RUNTIME_ERROR << "Index "<< ind << "given multiple times"; + l[i] *= la [ind]; + la[ind] = -1; + s[i] = std::min(s[i], a.indexmap().strides()[ind]); + } + } + + /// Check all indices have been used + for (int i = 0; i < la.size(); ++i) + if (la[i] != -1) TRIQS_RUNTIME_ERROR << "Grouping incomplete : index " << i << " is to taken into account"; + + /// Return the new view (indexmap, storage) + return {{l, s, 0}, a.storage()}; + } +} +} diff --git a/triqs/arrays/impl/assignment.hpp b/triqs/arrays/impl/assignment.hpp index dfaafd7f..0d752bd7 100644 --- a/triqs/arrays/impl/assignment.hpp +++ b/triqs/arrays/impl/assignment.hpp @@ -114,9 +114,19 @@ namespace triqs { namespace arrays { TRIQS_REJECT_MATRIX_COMPOUND_MUL_DIV_NON_SCALAR; typedef typename LHS::value_type value_type; LHS & lhs; const RHS & rhs; - impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_) {} //, p(*(lhs_.data_start())) {} + impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_) {} template void operator()(Args const & ... args) const { _ops_::invoke(lhs(args...),rhs(args...));} - void invoke() { foreach(lhs,*this); } + FORCEINLINE void invoke() { foreach(lhs,*this); } + }; + + // help compiler : in the most common case, less things to inline.. + template + struct impl::value && (!is_scalar_for::value) && (!is_isp::value)&& (!is_special::value)) > { + TRIQS_REJECT_ASSIGN_TO_CONST; + typedef typename LHS::value_type value_type; + LHS & lhs; const RHS & rhs; + impl(LHS & lhs_, const RHS & rhs_): lhs(lhs_), rhs(rhs_) {} + FORCEINLINE void invoke() { assign_foreach(lhs, rhs); } }; // ----------------- assignment for scalar RHS, except some matrix case -------------------------------------------------- diff --git a/triqs/arrays/impl/flags.hpp b/triqs/arrays/impl/flags.hpp deleted file mode 100644 index 4cf098b7..00000000 --- a/triqs/arrays/impl/flags.hpp +++ /dev/null @@ -1,87 +0,0 @@ -/******************************************************************************* - * - * TRIQS: a Toolbox for Research in Interacting Quantum Systems - * - * 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 - * 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 . - * - ******************************************************************************/ -#ifndef TRIQS_ARRAYS_FLAGS_H -#define TRIQS_ARRAYS_FLAGS_H -#include "../indexmaps/permutation.hpp" -namespace triqs { namespace arrays { - - typedef unsigned long long ull_t; - - namespace Tag {struct no_init {}; struct default_init {};} - - // Flags is a 64 bit unsigned int. - // 0 is the default option. - // Meaning of the bits : - // 0 -> Const - // 1,2 -> Predefined order : - // 0 : None, 1 : C, 2 : Fortran - // 3 -> Init the array - // 5 -> Boundcheck - - constexpr ull_t TraversalOrderC = 1ull << 1; - constexpr ull_t TraversalOrderFortran = 1ull << 2; - constexpr ull_t DefaultInit = 1ull << 3; - constexpr ull_t BoundCheck = 1ull << 5; - -#define BOUND_CHECK triqs::arrays::BoundCheck -#define TRAVERSAL_ORDER_C triqs::arrays::TraversalOrderC -#define TRAVERSAL_ORDER_FORTRAN triqs::arrays::TraversalOrderFortran -#define DEFAULT_INIT triqs::arrays::DefaultInit - - // NB : flags MUST be insensitive to slicing ... - // i.e. when I slice, the flags does not change. - - namespace flags { - - constexpr ull_t get(ull_t f, ull_t a) { return ((f & (1ull<> a);} - constexpr ull_t get2(ull_t f, ull_t a) { return ((f & ((1ull<> a);} - -#ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK - constexpr bool bound_check (ull_t f) { return true;} -#else - constexpr bool bound_check (ull_t f) { return get (f, 5)!=0;} -#endif - - constexpr bool traversal_order_c (ull_t f) { return get (f,1ull)!=0ull;} - constexpr bool traversal_order_fortran (ull_t f) { return get (f,2ull)!=0ull;} - - template< ull_t F> struct bound_check_trait { static constexpr bool value = bound_check(F); }; - - constexpr ull_t init_mode (ull_t f) { return get (f,3);} - - template struct init_tag1; - template<> struct init_tag1<0> { typedef Tag::no_init type;}; - template<> struct init_tag1<1> { typedef Tag::default_init type;}; - - // for the init_tag, we pass the *whole* option flag. - template struct init_tag : init_tag1 < init_mode(F)> {}; - - template struct assert_make_sense { - 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_assert ( (!( is_c && is_f)), "You asked C and Fortran traversal order at the same time..."); - static_assert ( (!( (is_c || is_f) && To )), "You asked C or Fortran traversal order and gave a traversal order ..."); - }; - - } -}}//namespace triqs::arrays -#endif diff --git a/triqs/arrays/impl/indexmap_storage_pair.hpp b/triqs/arrays/impl/indexmap_storage_pair.hpp index 36753c75..618c1f74 100644 --- a/triqs/arrays/impl/indexmap_storage_pair.hpp +++ b/triqs/arrays/impl/indexmap_storage_pair.hpp @@ -18,12 +18,11 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAY_IMPL_INDEX_STORAGE_PAIR_H -#define TRIQS_ARRAY_IMPL_INDEX_STORAGE_PAIR_H +#pragma once #include "./common.hpp" -#include "./flags.hpp" #include "../storages/shared_block.hpp" #include "./assignment.hpp" +#include "./print.hpp" #include "../indexmaps/cuboid/foreach.hpp" #include "triqs/utility/exceptions.hpp" #include "triqs/utility/typeid_name.hpp" @@ -36,7 +35,7 @@ namespace triqs { namespace arrays { - template auto get_shape (A const & x) DECL_AND_RETURN(x.domain().lengths()); + template AUTO_DECL get_shape (A const & x) RETURN(x.domain().lengths()); template size_t first_dim (A const & x) { return x.domain().lengths()[0];} template size_t second_dim (A const & x) { return x.domain().lengths()[1];} @@ -45,23 +44,25 @@ namespace triqs { namespace arrays { template size_t fifth_dim (A const & x) { return x.domain().lengths()[4];} template size_t sixth_dim (A const & x) { return x.domain().lengths()[5];} template size_t seventh_dim (A const & x) { return x.domain().lengths()[6];} + template size_t eighth_dim (A const & x) { return x.domain().lengths()[7];} + template size_t ninth_dim (A const & x) { return x.domain().lengths()[8];} template class iterator_adapter; - template struct ISPViewType; + template struct ISPViewType; - template + template class indexmap_storage_pair : Tag::indexmap_storage_pair, TRIQS_CONCEPT_TAG_NAME(MutableCuboidArray) { public : - typedef typename StorageType::value_type value_type; + using value_type=typename StorageType::value_type ; static_assert(std::is_constructible::value, "array/array_view and const operate only on values"); static_assert(!std::is_const::value, "no const type"); - typedef StorageType storage_type; - typedef IndexMapType indexmap_type; - static constexpr unsigned int rank = IndexMapType::domain_type::rank; - static constexpr ull_t opt_flags = OptionFlags; - static constexpr ull_t traversal_order = TraversalOrder; + using storage_type=StorageType ; + using indexmap_type=IndexMapType ; + using traversal_order_t = typename _get_traversal_order_t::type; + using traversal_order_arg = TraversalOrder; + static constexpr int rank = IndexMapType::domain_type::rank; static constexpr bool is_const = IsConst; protected: @@ -72,30 +73,25 @@ namespace triqs { namespace arrays { // ------------------------------- constructors -------------------------------------------- indexmap_storage_pair() {} - indexmap_storage_pair (indexmap_type const & IM, storage_type const & ST): indexmap_(IM),storage_(ST) {deleg(IM,ST);} - indexmap_storage_pair (indexmap_type && IM, storage_type && ST) : indexmap_(std::move(IM)),storage_(std::move(ST)){deleg(IM,ST);} - indexmap_storage_pair (indexmap_type const & IM, storage_type && ST) : indexmap_(IM),storage_(std::move(ST)) {deleg(IM,ST);} - - /// The storage is allocated from the size of IM. - indexmap_storage_pair (const indexmap_type & IM): indexmap_(IM),storage_(){ - this->storage_ = StorageType(this->indexmap_.domain().number_of_elements(), typename flags::init_tag::type() ); - } - private : - void deleg (const indexmap_type & IM, const storage_type & ST) { - //std::cerr << " construct ISP && ST "<< storage_type::is_weak<< std::endl; + indexmap_storage_pair(indexmap_type IM, storage_type ST) : indexmap_(std::move(IM)), storage_(std::move(ST)) { #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 - } + } - public: - /// Shallow copy - indexmap_storage_pair(const indexmap_storage_pair & X):indexmap_(X.indexmap()),storage_(X.storage_){} - indexmap_storage_pair(indexmap_storage_pair && X):indexmap_(std::move(X.indexmap())),storage_(std::move(X.storage_)){} - protected: + /// The storage is allocated from the size of IM. + indexmap_storage_pair(const indexmap_type &IM) : indexmap_(IM), storage_() { + storage_ = StorageType(indexmap_.domain().number_of_elements()); + } + public: + // Shallow copy + indexmap_storage_pair(indexmap_storage_pair const &X) = default; + indexmap_storage_pair(indexmap_storage_pair &&X) = default; + + protected: #ifdef TRIQS_WITH_PYTHON_SUPPORT indexmap_storage_pair (PyObject * X, bool enforce_copy, const char * name ) { try { @@ -109,11 +105,9 @@ namespace triqs { namespace arrays { // linker search for IndexMapType::domain_type::rank in triqs.so // and cannot resolve it ??? //<<"\n rank = "<< IndexMapType::domain_type::rank//this->rank - <<"\n OptionFlags = "<< OptionFlags <<"\nfrom the python object \n"<< numpy_interface::object_to_string(X) <<"\nThe error was :\n "<::type type; + typedef typename ISPViewType < V2, IM2::domain_type::rank, typename IM2::traversal_order_in_template, ViewTag, + ForceBorrowed || StorageType::is_weak, is_const > ::type type; }; template // non const version @@ -243,13 +238,13 @@ namespace triqs { namespace arrays { } /// Equivalent to make_view - typename ISPViewType::type + typename ISPViewType::type operator()() const & { return *this; } - typename ISPViewType::type + typename ISPViewType::type operator()() & { return *this; } - typename ISPViewType::type + typename ISPViewType::type operator()() && { return *this; } // Interaction with the CLEF library : calling with any clef expression as argument build a new clef expression @@ -297,8 +292,8 @@ namespace triqs { namespace arrays { // template friend void triqs_clef_auto_assign (indexmap_storage_pair & x, Fnt f) { assign_foreach(x,f);} // ------------------------------- Iterators -------------------------------------------- - typedef iterator_adapter const_iterator; - typedef iterator_adapter iterator; + using const_iterator=iterator_adapter ; + using iterator=iterator_adapter ; const_iterator begin() const {return const_iterator(indexmap(),storage(),false);} const_iterator end() const {return const_iterator(indexmap(),storage(),true);} const_iterator cbegin() const {return const_iterator(indexmap(),storage(),false);} @@ -311,11 +306,11 @@ namespace triqs { namespace arrays { // ------------------------------- resize -------------------------------------------- // void resize (domain_type const & d) { - this->indexmap_ = IndexMapType(d,this->indexmap_.memory_indices_layout()); + indexmap_ = IndexMapType(d,indexmap_.get_memory_layout()); // build a new one with the lengths of IND BUT THE SAME layout ! // optimisation. Construct a storage only if the new index is not compatible (size mismatch). - if (this->storage_.size() != this->indexmap_.domain().number_of_elements()) - this->storage_ = StorageType(this->indexmap_.domain().number_of_elements(), typename flags::init_tag::type() ); + if (storage_.size() != indexmap_.domain().number_of_elements()) + storage_ = StorageType(indexmap_.domain().number_of_elements()); } template @@ -332,11 +327,10 @@ namespace triqs { namespace arrays { // pretty print of the array friend std::ostream & operator << (std::ostream & out, const indexmap_storage_pair & A) { if (A.storage().size()==0) out<<"empty "; - else indexmaps::pretty_print(out,A); + else pretty_print(out,A); return out; } };// end class }}//namespace triqs::arrays -#endif diff --git a/triqs/arrays/impl/print.hpp b/triqs/arrays/impl/print.hpp new file mode 100644 index 00000000..96b31122 --- /dev/null +++ b/triqs/arrays/impl/print.hpp @@ -0,0 +1,77 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include +namespace triqs { +namespace arrays { + + /// ------------ Pretty Printing : specializing the default behaviour for d=1,2 ------------------------- + namespace PrettyPrint_details { + template struct print_impl { + std::ostream &out; + A const &a; + print_impl(std::ostream &out_, A const &a_) : out(out_), a(a_) {} + template void operator()(A0 const &a0, Args const &... args) const { + out << a(a0, args...) << " "; + } + void print() const { + out << "["; + // indexmaps::cuboid::foreach(a.domain(), + // std::cref(*this)); // foreach(a, std::cref(*this)); + indexmaps::cuboid::foreach_impl(_traversal_fortran{}, a.domain(), FORTRAN_LAYOUT, *this); + // foreach(a, *this); + out << "]"; + } + }; + + template struct print_impl<1, A> { + std::ostream &out; + A const &a; + print_impl(std::ostream &out_, A const &a_) : out(out_), a(a_) {} + void print() const { + auto d = a.indexmap().domain(); + out << "["; + for (size_t i = 0; i < d.lengths()[0]; ++i) out << (i > 0 ? "," : "") << a(i); + out << "]"; + } + }; + + template struct print_impl<2, A> { + std::ostream &out; + A const &a; + print_impl(std::ostream &out_, A const &a_) : out(out_), a(a_) {} + void print() const { + auto d = a.indexmap().domain(); + out << "\n["; + for (size_t i = 0; i < d.lengths()[0]; ++i) { + out << (i == 0 ? "[" : " ["); + for (size_t j = 0; j < d.lengths()[1]; ++j) out << (j > 0 ? "," : "") << a(i, j); + out << "]" << (i == d.lengths()[0] - 1 ? "" : "\n"); + } + out << "]"; + } + }; + } + template void pretty_print(std::ostream &out, A const &a) { + PrettyPrint_details::print_impl(out, a).print(); + } +} +} diff --git a/triqs/arrays/indexmaps/cuboid/domain.hpp b/triqs/arrays/indexmaps/cuboid/domain.hpp index 90727b6d..1b9d5d35 100644 --- a/triqs/arrays/indexmaps/cuboid/domain.hpp +++ b/triqs/arrays/indexmaps/cuboid/domain.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011-2013 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,180 +18,132 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_CUBOID_DOMAIN_H -#define TRIQS_ARRAYS_INDEXMAP_CUBOID_DOMAIN_H +#pragma once #include "../common.hpp" #include "../range.hpp" #include "./mem_layout.hpp" +#include "../permutation.hpp" #include #include "../../impl/exceptions.hpp" #include #include -#include -#include -#include -#include -namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { - using namespace triqs::arrays::permutations; +namespace triqs { +namespace arrays { + namespace indexmaps { + namespace cuboid { + using namespace triqs::arrays::permutations; - /// Standard hyper_rectangular domain for arrays - template - class domain_t { - typedef mini_vector n_uple; - n_uple lengths_; - friend class boost::serialization::access; - template void serialize(Archive & ar, const unsigned int version) { ar & TRIQS_MAKE_NVP("dimensions",lengths_);} - public : - static constexpr unsigned int rank = Rank; - typedef n_uple index_value_type; + /// Standard hyper_rectangular domain for arrays + template class domain_t { + using n_uple = mini_vector; + n_uple lengths_; + friend class boost::serialization::access; + template void serialize(Archive &ar, const unsigned int version) { + ar &TRIQS_MAKE_NVP("dimensions", lengths_); + } - domain_t () =default; - domain_t (const domain_t & C) = default; - domain_t (domain_t && C) = default; - domain_t & operator =( domain_t const &) = default; - domain_t & operator =( domain_t && x)= default; + public: + static constexpr unsigned int rank = Rank; + using index_value_type = n_uple; - domain_t (n_uple lengths):lengths_(std::move(lengths)) {} - domain_t (mini_vector const & lengths):lengths_(lengths) {} - domain_t (std::vector const & l):lengths_() { - if (!(l.size()==Rank)) TRIQS_RUNTIME_ERROR << "cuboid domain_t construction : vector size incorrect : got "< domain_t(size_t i0, T... t): lengths_(i0, t...){} + domain_t() = default; + domain_t(const domain_t &C) = default; + domain_t(domain_t &&C) = default; + domain_t &operator=(domain_t const &) = default; + domain_t &operator=(domain_t &&x) = default; - size_t number_of_elements() const { return lengths_.product_of_elements();} - bool operator==(domain_t const & X) const { return this->lengths_ == X.lengths_;} - bool operator!=(domain_t const & X) const { return !(*this==X);} - n_uple const & lengths() const & { return lengths_;} - n_uple lengths() && { return lengths_;} - /** Generates the value of the indices of a cuboid_domain. */ - static constexpr ull_t iteration_order_default = permutations::identity(Rank); - template - class gal_generator { - typedef index_value_type indices_type; - const domain_t * dom; + domain_t(n_uple lengths) : lengths_(std::move(lengths)) {} + domain_t(mini_vector const &lengths) : lengths_(lengths) {} + domain_t(std::vector const &l) : lengths_() { + if (!(l.size() == Rank)) + TRIQS_RUNTIME_ERROR << "cuboid domain_t construction : vector size incorrect : got " << l.size() << " while expected " + << Rank; + lengths_ = n_uple(l); + } + template domain_t(size_t i0, T... t) : lengths_(i0, t...) {} + + size_t number_of_elements() const { return lengths_.product_of_elements(); } + bool operator==(domain_t const &X) const { return this->lengths_ == X.lengths_; } + bool operator!=(domain_t const &X) const { return !(*this == X); } + n_uple const &lengths() const &{ return lengths_; } + n_uple lengths() && { return lengths_; } + /** Generates the value of the indices of a cuboid_domain. */ + static constexpr ull_t iteration_order_default = permutations::identity(Rank); + template class gal_generator { + using indices_type = index_value_type; + const domain_t *dom; indices_type indices_tuple; bool atend; + public: - gal_generator (const domain_t & P, bool atEnd=false): dom(&P), atend(atEnd) {} - bool operator==(const gal_generator & IT2) const { assert((IT2.dom == dom)); return ((IT2.atend==atend) );} - bool operator!=(const gal_generator & IT2) const { return (!operator==(IT2));} - indices_type const & operator *() const { return indices_tuple;} - explicit operator bool () const { return !atend;} - gal_generator & operator++(){ assert(!atend); atend = advance_impl(std::integral_constant()); return *this;} - private: - template bool advance_impl(std::integral_constant) { - constexpr int p = permutations::apply(IterationOrder, r); - if (indices_tuple[p] < dom->lengths()[p]-1) { ++(indices_tuple[p]); return false;} - indices_tuple[p] = 0; - return advance_impl(std::integral_constant()); + gal_generator(const domain_t &P, bool atEnd = false) : dom(&P), atend(atEnd) {} + bool operator==(const gal_generator &IT2) const { + assert((IT2.dom == dom)); + return ((IT2.atend == atend)); } - bool advance_impl(std::integral_constant) { return true;} + bool operator!=(const gal_generator &IT2) const { return (!operator==(IT2)); } + indices_type const &operator*() const { return indices_tuple; } + explicit operator bool() const { return !atend; } + gal_generator &operator++() { + assert(!atend); + atend = advance_impl(std::integral_constant()); + return *this; + } + + private: + template bool advance_impl(std::integral_constant) { + constexpr int p = permutations::apply(IterationOrder, r); + if (indices_tuple[p] < dom->lengths()[p] - 1) { + ++(indices_tuple[p]); + return false; + } + indices_tuple[p] = 0; + return advance_impl(std::integral_constant()); + } + bool advance_impl(std::integral_constant) { return true; } }; - typedef gal_generator<> generator; + using generator = gal_generator<>; - generator begin() const { return generator(*this,false);} - generator end() const { return generator(*this,true);} - /* End of generator */ + generator begin() const { return generator(*this, false); } + generator end() const { return generator(*this, true); } + /* End of generator */ - // Check that key in in the domain - template void assert_key_in_domain(KeyType const & key) const { - std::stringstream fs; - bool res = key_check_impl(std::integral_constant(), key,this->lengths_,fs); - if (!res) TRIQS_ARRAYS_KEY_ERROR << " key out of domain \n" < - bool key_check_impl (std::integral_constant, KeyType const & key, n_uple const & L, std::stringstream & fs ) const { - //bool cond = ( ( size_t(key[r]) < L[r])); - bool cond = ( ( size_t(std::get(key)) < L[r])); - //if (!cond) fs << "key ["<(), key,L,fs) && cond; + // Check that key in in the domain + template void assert_key_in_domain(KeyType const &key) const { + std::stringstream fs; + bool res = key_check_impl(std::integral_constant(), key, this->lengths_, fs); + if (!res) TRIQS_ARRAYS_KEY_ERROR << " key out of domain \n" << fs.str(); } - template bool key_check_impl (std::integral_constant, KeyType const &, n_uple const &, std::stringstream &) const { return true;} - - // Check that key in in the domain : variadic form. No need for speed optimisation here, it is just for debug - template void assert_key_in_domain_v (Args const & ... args) const { assert_key_in_domain( std::make_tuple(args...));} - - friend std::ostream & operator<<(std::ostream & out, domain_t const & x){return out<<"Cuboid of rank "< void foreach(domain_t const & dom, FntType F) {\ - BOOST_PP_REPEAT(RR,AUX0,BOOST_PP_DEC(RR))\ - mini_vector t;\ - const mini_vector l(dom.lengths());\ - BOOST_PP_REPEAT(RR,AUX1,nil){ F(BOOST_PP_REPEAT(RR,AUX3,nil)); }} - BOOST_PP_REPEAT_FROM_TO(1,ARRAY_NRANK_MAX , IMPL, nil); -#undef IMPL -#undef AUX0 -#undef AUX1 -#undef AUX3 -} -/// ------------ Pretty Printing : specializing the default behaviour for d=1,2 ------------------------- -namespace PrettyPrint_details { - template - struct print_impl { - std::ostream & out; A const & a; - print_impl( std::ostream & out_, A const & a_) : out(out_), a(a_){} - template void operator()(A0 const & a0, Args const & ... args) const { out << a(a0,args...)<< " ";} - void print() const { out<<"["; - indexmaps::cuboid::foreach (a.domain(), std::cref(*this)); //foreach(a, std::cref(*this)); - out<<"]"; } - }; - - template - struct print_impl <1,A> { - std::ostream & out; A const & a; - print_impl( std::ostream & out_, A const & a_) : out(out_), a(a_){} - void print() const { - auto d = a.indexmap().domain(); - out<<"["; for (size_t i=0; i< d.lengths()[0]; ++i) out<<(i>0 ? ",": "")< - struct print_impl <2,A> { - std::ostream & out; A const & a; - print_impl( std::ostream & out_, A const & a_) : out(out_), a(a_){} - void print() const { - auto d = a.indexmap().domain(); - out<<"\n["; - for (size_t i=0; i< d.lengths()[0]; ++i) { - out<<(i==0 ? "[" : " ["); - for (size_t j=0; j< d.lengths()[1]; ++j) out<<(j>0 ? ",": "")< + bool key_check_impl(std::integral_constant, KeyType const &key, n_uple const &L, std::stringstream &fs) const { + bool cond = ((size_t(std::get(key)) < L[r])); + if (!cond) fs << "key [" << r << "] = " << std::get(key) << " is not within [0," << L[r] << "[\n"; + return key_check_impl(std::integral_constant(), key, L, fs) && cond; + } + template + bool key_check_impl(std::integral_constant, KeyType const &, n_uple const &, std::stringstream &) const { + return true; } - out<<"]"; - } - }; + // Check that key in in the domain : variadic form. No need for speed optimisation here, it is just for debug + template void assert_key_in_domain_v(Args const &... args) const { + assert_key_in_domain(std::make_tuple(args...)); + } + + friend std::ostream &operator<<(std::ostream &out, domain_t const &x) { + return out << "Cuboid of rank " << x.rank << " and dimensions " << x.lengths(); + } + }; + } + } + + template indexmaps::cuboid::domain_t make_cuboid_domain(U... u) { + return {size_t(u)...}; + } + + using matrix_shape_t = indexmaps::cuboid::domain_t<2>; + using vector_shape_t = indexmaps::cuboid::domain_t<1>; } -template -void pretty_print (std::ostream & out, A const & a ) { PrettyPrint_details::print_impl(out,a).print();} - } - -template -indexmaps::cuboid::domain_t make_cuboid_domain(U ... u) { return {size_t(u)...};} -//cuboid_array_domain make_cuboid_domain(U ... u) { return {u...};} - -typedef indexmaps::cuboid::domain_t<2> matrix_shape_t; -typedef indexmaps::cuboid::domain_t<1> vector_shape_t; - - -}} -#endif diff --git a/triqs/arrays/indexmaps/cuboid/foreach.hpp b/triqs/arrays/indexmaps/cuboid/foreach.hpp index a3d50f7d..3c94317a 100644 --- a/triqs/arrays/indexmaps/cuboid/foreach.hpp +++ b/triqs/arrays/indexmaps/cuboid/foreach.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,30 +18,125 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_CUBOID_FOREACH_H -#define TRIQS_ARRAYS_INDEXMAP_CUBOID_FOREACH_H +#pragma once #include "./map.hpp" +#include "./mem_layout.hpp" +#include +#include +#include +#include -namespace triqs { namespace arrays { +namespace triqs { +namespace arrays { - template struct get_traversal_order : indexmaps::mem_layout::fortran_order_tr{}; - //template struct get_traversal_order : indexmaps::mem_layout::c_order_tr{}; - template struct get_traversal_order : std::integral_constant{}; +#define FORCEINLINE __inline__ __attribute__((always_inline)) + // ------------------------- foreach ----------------------------------------------------- + namespace indexmaps { + namespace cuboid { + + // using foreach_int_type=size_t ; + using foreach_int_type = std::ptrdiff_t; +// better to be signed here : 1) on some machine/compiler, it is a lot faster ! +// When used with clef auto assign, e.g. A(i_,j_) = i -2*j, one needs signed arithmetics + +#define AUX_FOR(X) for (ind[X] = 0; ind[X] < l[X]; ++ind[X]) +#define AUX_C(z, P, unused) for (foreach_int_type i##P = 0; i##P < l[P]; ++i##P) +#define AUX_F(z, P, RANK) AUX_FOR(RANK - P) +#define AUX_Dynamical(z, P, unused) AUX_FOR(ml[P]) +#define AUX_Custom1(z, P, NNN) constexpr int p##P = permutations::apply(traversal_order_perm, P); +#define AUX_Custom2(z, P, unused) AUX_FOR(p##P) +#define AUX3(z, p, unused) BOOST_PP_COMMA_IF(p) ind[p] +#define AUX3C(z, p, unused) BOOST_PP_COMMA_IF(p) i##p +#define IMPL(z, RR, unused) \ + template \ + FORCEINLINE void foreach_impl(_traversal_c, domain_t const& dom, memory_layout const& ml, FntType F) { \ + const mini_vector l(dom.lengths()); \ + BOOST_PP_REPEAT(RR, AUX_C, nil) { F(BOOST_PP_REPEAT(RR, AUX3C, nil)); } \ + } \ + template \ + FORCEINLINE void foreach_impl(_traversal_fortran, domain_t const& dom, memory_layout ml, FntType F) { \ + foreach_int_type ind[RR]; \ + const mini_vector l(dom.lengths()); \ + BOOST_PP_REPEAT(RR, AUX_F, BOOST_PP_DEC(RR)) { F(BOOST_PP_REPEAT(RR, AUX3, nil)); } \ + } \ + template \ + FORCEINLINE void foreach_impl(_traversal_dynamical, domain_t const& dom, memory_layout ml, FntType F) { \ + foreach_int_type ind[RR]; \ + const mini_vector l(dom.lengths()); \ + BOOST_PP_REPEAT(RR, AUX_Dynamical, nil) { F(BOOST_PP_REPEAT(RR, AUX3, nil)); } \ + } \ + template \ + FORCEINLINE void foreach_impl(_traversal_custom, domain_t const& dom, memory_layout ml, FntType F) { \ + constexpr ull_t traversal_order_perm = permutations::permutation(Is...); \ + BOOST_PP_REPEAT(RR, AUX_Custom1, BOOST_PP_DEC(RR)); \ + foreach_int_type ind[RR]; \ + const mini_vector l(dom.lengths()); \ + BOOST_PP_REPEAT(RR, AUX_Custom2, nil) { F(BOOST_PP_REPEAT(RR, AUX3, nil)); } \ + } + BOOST_PP_REPEAT_FROM_TO(1, ARRAY_NRANK_MAX, IMPL, nil); +#undef IMPL +#undef AUX_C +#undef AUX_F +#undef AUX_FOR +#undef AUX_Dynamical +#undef AUX_Custom +#undef AUX_Custom1 +#undef AUX3 +#undef AUX3C + } + } + + /// Get the traversal order + template struct _get_traversal_order { + using traversal_order_t = _traversal_c; + static memory_layout invoke(A const& a) { + return memory_layout{}; + } + }; + template struct _get_traversal_order { + using traversal_order_t = typename A::traversal_order_t; + static memory_layout invoke(A const& a) { return a.indexmap().get_memory_layout(); } + }; + + /// --------------- FOREACH ------------------------ template - typename std::enable_if::value >::type - foreach (T const & x, Function const & F) { indexmaps::cuboid::foreach::value> (x.domain(), F); } - - template - struct assign_foreach_adapter { - T& x; Function const & f; - assign_foreach_adapter( T& x_, Function const & ff): x(x_), f(ff){} - template void operator()(Args const & ... args) const { x(args...) = f(args...);} - }; - - template - typename std::enable_if::value >::type - assign_foreach (T & x, Function const & F) { indexmaps::cuboid::foreach::value> (x.domain(),assign_foreach_adapter(x,F)); } - -}}//namespace + FORCEINLINE std14::enable_if_t::value> foreach(T const& x, Function const& F) { + using S = _get_traversal_order; +#ifndef TRIQS_ARRAYS_FOREACH_C_OR_DYNAMICAL + indexmaps::cuboid::foreach_impl(typename S::traversal_order_t{}, x.domain(), S::invoke(x), F); +#else + if (ml.is_c()) { + indexmaps::cuboid::foreach_impl(_traversal_c{}, x.domain(), S::invoke(x), F); + } else { + indexmaps::cuboid::foreach_impl(_traversal_dynamical{}, x.domain(), S::invoke(x), F); + } #endif + } + +/// --------------- ASSIGN FOREACH ------------------------ +#ifndef TRIQS_C11 + template + std14::enable_if_t::value> assign_foreach(T& x, Function const& f) { + using S = _get_traversal_order; + indexmaps::cuboid::foreach_impl(typename S::traversal_order_t{}, x.domain(), S::invoke(x), + [&x, &f](auto const&... args) { x(args...) = f(args...); }); + } + +#else + template struct assign_foreach_adapter { + T& x; + Function const& f; + assign_foreach_adapter(T& x_, Function const& ff) : x(x_), f(ff) {} + template void operator()(Args const&... args) const { x(args...) = f(args...); } + }; + + template + std14::enable_if_t::value> assign_foreach(T& x, Function const& F) { + using S = _get_traversal_order; + indexmaps::cuboid::foreach_impl(typename S::traversal_order_t{}, x.domain(), S::invoke(x), + assign_foreach_adapter(x, F)); + } +#endif +} +} // namespace diff --git a/triqs/arrays/indexmaps/cuboid/group_indices.hpp b/triqs/arrays/indexmaps/cuboid/group_indices.hpp deleted file mode 100644 index b800fef7..00000000 --- a/triqs/arrays/indexmaps/cuboid/group_indices.hpp +++ /dev/null @@ -1,107 +0,0 @@ -/******************************************************************************* - * - * TRIQS: a Toolbox for Research in Interacting Quantum Systems - * - * Copyright (C) 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 - * 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 . - * - ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_CUBOID_GROUP_INDICES_H -#define TRIQS_ARRAYS_INDEXMAP_CUBOID_GROUP_INDICES_H -#include -namespace triqs { namespace arrays { - - // a variadic push_back into a vector. - template void vector_push_back_v(V & v){} - template void vector_push_back_v(V & v, T0 && t0, T && ... t){ v.push_back(std::forward(t0)); vector_push_back_v(v,t...);} - - // a list of compile time int... - template struct m_index{}; - template struct m_index {static constexpr int head = I0; }; - template std::vector m_index_to_vector( m_index) { std::vector v; vector_push_back_v(v,Is...); return v;} - - // a trait to get the min, max of a m_index - template struct get_min_max; - template struct get_min_max> { - typedef get_min_max> r; - static constexpr int min = (I0 < r::min ? I0 : r::min); static constexpr int max = (I0 > r::max ? I0 : r::max); - }; - template struct get_min_max> { static constexpr int min = I0; static constexpr int max = I0; }; - - // given a m_index of indices, a metafunction to map to their position into memory - template struct index_group_to_mem_pos_list; - template struct index_group_to_mem_pos_list> { - typedef m_index < indexmaps::mem_layout::index_to_memory_rank(ML,Is)...> type; - static_assert( get_min_max::max - get_min_max::min + 1 == sizeof...(Is), "Indices not contiguous in memory"); - }; - -#ifndef TRIQS_WORKAROUND_INTEL_COMPILER_BUGS - // count the number of Is strictly smaller than C - constexpr int count_pos(int C){ return 0;} - template constexpr int count_pos(int C, int i, T... t){ return (i struct count_pos_list { typedef m_index< count_pos(Is, Is...)...> type; }; -#else - template struct count_pos; - template struct count_pos> : std::integral_constant>::value + (I0 {}; - template struct count_pos> : std::integral_constant{}; - template struct count_pos_list { typedef m_index tmp; typedef m_index< count_pos::value...> type; }; -#endif - - // a simple foreach - template void for_each (Callee & C, T0 const & t0, T const & ... t ) { C(t0); for_each(C, t...);} - template void for_each (Callee & C){} - - // make a permutation out of a m_index - template struct to_permu; - template struct to_permu> { static constexpr ull_t value = permutations::permutation(Is...);}; - - // the main class - template struct group_indices_impl { - static constexpr int new_dim = sizeof...(MIndex); - typedef typename count_pos_list< index_group_to_mem_pos_list::type::head ...>::type new_memory_pos; - - static constexpr ull_t traversal_layout = to_permu::value; - typedef array_view type; - - static type invoke(A const & a) { - if (a.indexmap().memory_indices_layout_ull() != a.indexmap().traversal_order) - TRIQS_RUNTIME_ERROR << "Grouping indices is only possible for arrays when the memory_layout is the same as the traversal order \n" - << "But here your memory_layout is "<< a.indexmap().memory_indices_layout() << " while the traversal order is "<< a.indexmap().traversal_order_indices_layout(); - std::vector< std::vector > Indices; - vector_push_back_v(Indices, m_index_to_vector(MIndex())...); - mini_vector l; - mini_vector s; - for (size_t u=0; u - typename group_indices_impl::type group_indices(A const & a, MIndex... ) { return group_indices_impl::invoke(a);} -}}//namespace triqs::arrays -#endif diff --git a/triqs/arrays/indexmaps/cuboid/map.hpp b/triqs/arrays/indexmaps/cuboid/map.hpp index 2fa2ef05..14626f91 100644 --- a/triqs/arrays/indexmaps/cuboid/map.hpp +++ b/triqs/arrays/indexmaps/cuboid/map.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,50 +18,49 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_CUBOID_MAP_H -#define TRIQS_ARRAYS_INDEXMAP_CUBOID_MAP_H +#pragma once #include "./domain.hpp" #include "./mem_layout.hpp" -#include "../../impl/flags.hpp" #include #include -namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { +namespace triqs { namespace arrays { - template< bool BC, class D, class K> struct _chk; - template< class D, class K> struct _chk { static void invoke (D const & d, K const & k) {d.assert_key_in_domain(k);} }; - template< class D, class K> struct _chk { static void invoke (D const & d, K const & k) {} }; + struct _traversal_c {}; + struct _traversal_fortran {}; + struct _traversal_dynamical {}; + template struct _traversal_custom {}; - template< bool BC, class D, class ... K> struct _chk_v; - template< class D, class ... K> struct _chk_v { static void invoke (D const & d, K const & ... k) {d.assert_key_in_domain_v(k...);} }; - template< class D, class ... K> struct _chk_v { static void invoke (D const & d, K const & ... k) {} }; + template struct _get_traversal_order_t { + using type = T; + }; + template <> struct _get_traversal_order_t { + using type = _traversal_c; + }; - template ull_t memory_layout_from_strides(mini_vector const & strides) { - int c[Rank]; for (size_t i=0; i strides[j]) c[i]++; else c[j]++; - // we computed the map : index -> memory_rank, which is the inverse map, cf mem_layout - return permutations::inverse(permutations::permutation_from_array(c, Rank)); + constexpr ull_t _get_traversal_order_permutation(int R, _traversal_c) { return permutations::identity(R); } + constexpr ull_t _get_traversal_order_permutation(int R, _traversal_fortran) { return permutations::ridentity(R); } + template constexpr ull_t _get_traversal_order_permutation(int R,_traversal_custom) { + static_assert(sizeof...(Is) == R, " Rank mismatch"); + return permutations::permutation(Is...); } + namespace indexmaps { namespace cuboid { + /** Standard hyper_rectangular arrays, implementing the IndexMap concept. */ - template - class map { + template class map { public : - static constexpr bool CheckBounds = flags::bound_check_trait::value; - static constexpr ull_t traversal_order_in_template = TraversalOrder; - static constexpr ull_t traversal_order = indexmaps::mem_layout::get_traversal_order::value; - typedef void has_traversal_order_tag; - static const unsigned int rank = Rank; - typedef mini_vector lengths_type; - typedef mini_vector strides_type; - typedef domain_t domain_type; - domain_type const & domain() const { return mydomain;} + static const int rank = Rank; + using lengths_type=mini_vector ; + using strides_type=mini_vector ; + using domain_type = domain_t; + using traversal_order_in_template = TraversalOrder; + using has_traversal_order_tag = void; + domain_type const& domain() const { return mydomain; } // semi-regular type - map (): start_shift_(0), memory_order_(memory_layout(traversal_order)) {} + map (): start_shift_(0), memory_order_() {} map (map const & C) = default; map (map && C) = default; map & operator = (map const & m) = default; @@ -69,8 +68,8 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { // basic construction map(memory_layout const & ml):mydomain(), start_shift_(0), memory_order_(ml) {} - map(domain_type const & C): mydomain(C), start_shift_(0), memory_order_(traversal_order) {compute_stride_compact();} - map(domain_type const & C, memory_layout ml): mydomain(C), start_shift_(0), memory_order_(ml) {compute_stride_compact();} + map(domain_type const & C): mydomain(C), start_shift_(0), memory_order_() {compute_stride_compact();} + map(domain_type const & C, memory_layout ml): mydomain(C), start_shift_(0), memory_order_(std::move(ml)) {compute_stride_compact();} /// Construction from the length, the stride, start_shift map(lengths_type Lengths, strides_type strides, std::ptrdiff_t start_shift ): @@ -81,30 +80,45 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { map(lengths_type Lengths, strides_type strides, std::ptrdiff_t start_shift, memory_layout const & ml ): mydomain(std::move(Lengths)), strides_(std::move(strides)), start_shift_(start_shift), memory_order_ (ml) {} - /// Construction from another map with the same order (used in grouping indices) - template map (map const & C): - mydomain(C.domain()), strides_(C.strides()), start_shift_(C.start_shift()), memory_order_ (C.memory_indices_layout()) {} + /// Construction from another map with the same order + template + map(map const& C) + : mydomain(C.domain()), strides_(C.strides()), start_shift_(C.start_shift()), memory_order_(C.get_memory_layout()) {} + + template map& operator=(map const& m) { + *this = map{m}; + } // transposition - template friend map transpose(map const& m, INT... is) { - return map{{m.domain().lengths()[is]...}, {m.strides_[is]...}, m.start_shift_, transpose(m.memory_order_, is...)}; + friend map transpose(map const& m, mini_vector const & perm) { + lengths_type l; + strides_type s; + for (int u = 0; u < Rank; ++u) { + l[perm[u]] = m.domain().lengths()[u]; + s[perm[u]] = m.strides_[u]; + } + return map{l, s, m.start_shift_, transpose(m.memory_order_, perm)}; } /// Returns the shift in position of the element key. - template - size_t operator[] (KeyType const & key ) const { - _chk::invoke (this->domain(),key); - return start_shift_ + dot_product(key,this->strides()); + template size_t operator[](KeyType const& key) const { +#ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK + this->domain().assert_key_in_domain(key); +#endif + return start_shift_ + dot_product(key, this->strides()); } friend std::ostream & operator << (std::ostream & out, const map & P) { - return out <<" ordering = {"< size_t operator()(Args const & ... args) const { - _chk_v::invoke (this->domain(),args...); +#ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK + this->domain().assert_key_in_domain_v(args...); +#endif return start_shift_ + _call_impl<0>(args...); } private : @@ -115,19 +129,17 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { /// bool is_contiguous() const { - const size_t last_index = mem_layout::memory_rank_to_index(memory_order_.value, rank-1); - return (strides()[last_index] * this->lengths()[last_index] == mydomain.number_of_elements()); + int slowest_index = memory_order_[0]; + return (strides()[slowest_index] * this->lengths()[slowest_index] == mydomain.number_of_elements()); } size_t start_shift() const { return start_shift_;} lengths_type const & lengths() const { return mydomain.lengths();} strides_type const & strides() const { return this->strides_;} - memory_layout const & memory_indices_layout() const { return memory_order_;} - memory_layout traversal_order_indices_layout() const { return memory_layout(traversal_order);} - ull_t memory_indices_layout_ull() const { return memory_order_.value;} - bool memory_layout_is_c() const { return memory_indices_layout().value == mem_layout::c_order(Rank);} - bool memory_layout_is_fortran() const { return memory_indices_layout().value == mem_layout::fortran_order(Rank);} + memory_layout const & get_memory_layout() const { return memory_order_;} + bool memory_layout_is_c() const { return get_memory_layout().is_c();} + bool memory_layout_is_fortran() const { return get_memory_layout().is_fortran();} private : domain_type mydomain; @@ -143,19 +155,30 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { ar & TRIQS_MAKE_NVP("start_shift",start_shift_); } // for construction + // TODO : use tupletools void compute_stride_compact() { size_t str = 1; - csc_impl(memory_order_.value, str, std::integral_constant()); + csc_impl(memory_order_, str, std::integral_constant()); assert(this->domain().number_of_elements()==str); } - template - void csc_impl(ull_t order, size_t & str, std::integral_constant) { - size_t u = mem_layout::memory_rank_to_index(order, rank-v); - this->strides_[u] = str; - str *= this->lengths() [u]; - csc_impl(order, str, std::integral_constant()); - } - void csc_impl(ull_t order,size_t & str, std::integral_constant) {} + // call for indices fastest (rank -1) to slowest (0) + template void csc_impl(memory_layout const& ml, size_t& str, std::integral_constant) { + // size_t u = mem_layout::memory_rank_to_index(order, rank-v); + int u = ml[v - 1]; + this->strides_[u] = str; + str *= this->lengths()[u]; + csc_impl(ml, str, std::integral_constant()); + } + void csc_impl(memory_layout const&, size_t&, std::integral_constant) {} + + // iterator helper impl. + static constexpr int __iter_get_p(int p, map const * im, _traversal_c) { return p;} + static constexpr int __iter_get_p(int p, map const* im, _traversal_fortran) { return Rank - p - 1; } + static int __iter_get_p(int p, map const* im, _traversal_dynamical) { return im->get_memory_layout()[p]; } + template static constexpr int __iter_get_p(int p, map const* im, _traversal_custom) { + return permutations::apply(_get_traversal_order_permutation(Rank, _traversal_custom{}), + p); + } public: @@ -166,9 +189,9 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { */ class iterator : public boost::iterator_facade< iterator, const std::ptrdiff_t, boost::forward_traversal_tag > { public: - typedef map indexmap_type; - typedef typename domain_type::index_value_type indices_type; - typedef const std::ptrdiff_t return_type; + using indexmap_type=map ; + using indices_type=typename domain_type::index_value_type ; + using return_type=const std::ptrdiff_t ; iterator (): im(NULL), pos(0),atend(true) {} iterator (const map & P, bool atEnd=false, ull_t iteration_order=0): im(&P), pos(im->start_shift()),atend(atEnd) {} @@ -178,7 +201,7 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { friend class boost::iterator_core_access; void increment(){ inc_ind_impl (std::integral_constant()); } template inline void inc_ind_impl(std::integral_constant) { - constexpr size_t p = mem_layout::memory_rank_to_index(traversal_order, rank-v); + int p = __iter_get_p(v - 1, im, typename _get_traversal_order_t::type{}); #ifdef TRIQS_ARRAYS_ENFORCE_BOUNDCHECK if (atend) TRIQS_RUNTIME_ERROR << "Iterator in cuboid cannot be pushed after end !"; #endif @@ -197,19 +220,18 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { }; }; //------------- end class --------------------- - - }//namespace cuboid -template -bool compatible_for_assignment (const cuboid::map & X1, const cuboid::map & X2) { return X1.lengths() == X2.lengths();} +template +bool compatible_for_assignment(const cuboid::map& X1, const cuboid::map& X2) { + return X1.lengths() == X2.lengths(); +} -template - bool raw_copy_possible (const cuboid::map & X1,const cuboid::map & X2) { - return ( (X1.memory_indices_layout() == X2.memory_indices_layout()) +template +bool raw_copy_possible(const cuboid::map& X1, const cuboid::map& X2) { + return ( (X1.get_memory_layout() == X2.get_memory_layout()) && X1.is_contiguous() && X2.is_contiguous() && (X1.domain().number_of_elements()==X2.domain().number_of_elements())); } }}}//namespace triqs::arrays::indexmaps -#endif diff --git a/triqs/arrays/indexmaps/cuboid/mem_layout.hpp b/triqs/arrays/indexmaps/cuboid/mem_layout.hpp index f8c83e96..d13ae92a 100644 --- a/triqs/arrays/indexmaps/cuboid/mem_layout.hpp +++ b/triqs/arrays/indexmaps/cuboid/mem_layout.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,93 +18,110 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_MEMORY_LAYOUT_H -#define TRIQS_ARRAYS_INDEXMAP_MEMORY_LAYOUT_H -#include "../permutation.hpp" -#include "../../impl/flags.hpp" +#pragma once +namespace triqs { +namespace arrays { -namespace triqs { namespace arrays { + struct __mem_layout_fortran_order {}; + struct __mem_layout_c_order {}; +#define C_LAYOUT (triqs::arrays::__mem_layout_c_order()) +#define FORTRAN_LAYOUT (triqs::arrays::__mem_layout_fortran_order()) - namespace indexmaps { namespace mem_layout { - /* The storage order is given by a permutation P stored in a ull_t (unsigned long long) as in permutations::.. - * P[0] : the fastest index, - * P[RANK-1] : the slowest index - * Example : - * 012 : Fortran, the first index is the fastest - * 210: C the last index is the fastest - * 120 : storage (i,j,k) is : index j is fastest, then k, then i - * - * index_to_memory_rank : i ---> r : to the index (0,1, ... etc), associates the rank in memory - * e.g. r=0 : fastest index, r = RANK-1 : the slowest - * memory_rank_to_index : the inverse mapping : r---> i : - * 0-> the fastest index, etc.. - * - * All these computations can be done *at compile time* (constexpr) - */ + enum class memory_order_type { + C, + Fortran, + Custom + }; - constexpr int memory_rank_to_index(ull_t p, int r) { return permutations::apply(p, r);} - constexpr int index_to_memory_rank(ull_t p, int r) { return permutations::apply(permutations::inverse(p), r);} + /* The storage order is given by a permutation P + * P[0] : the slowest index, + * P[Rank-1] : the fastest index + * Example : + * 210 : Fortran, the first index is the fastest + * 012 : C the last index is the fastest + * 120 : storage (i,j,k) is : index j is slowest, then k, then i + */ + template class memory_layout { - constexpr bool is_fortran (ull_t p){ return p == permutations::identity(permutations::size(p));} - constexpr bool is_c (ull_t p){ return p == permutations::ridentity(permutations::size(p));} + memory_order_type order = memory_order_type::C; + mini_vector p = utility::no_init_tag{}; - constexpr ull_t fortran_order (int n){ return permutations::identity(n);} - constexpr ull_t c_order (int n){ return permutations::ridentity(n);} - - template struct fortran_order_tr { static constexpr ull_t value = permutations::identity(n);}; - template struct c_order_tr { static constexpr ull_t value = permutations::ridentity(n);}; + public: + static const int rank = Rank; - // From the flag in the template definition to the real traversal_order - // 0 -> C order - // 1 -> Fortran Order - // Any other number interpreted as a permutation ? - - constexpr ull_t _get_traversal_order (int rank, ull_t fl, ull_t to) { return (flags::traversal_order_c(fl) ? c_order(rank) : - (flags::traversal_order_fortran(fl) ? fortran_order(rank) : (to==0 ? c_order(rank) : to )));} - - template< int rank, ull_t fl, ull_t to> struct get_traversal_order { static constexpr ull_t value = _get_traversal_order (rank,fl,to); }; - }} - - struct memory_layout_fortran {}; - struct memory_layout_c {}; - -#define FORTRAN_LAYOUT (triqs::arrays::memory_layout_fortran()) -#define C_LAYOUT (triqs::arrays::memory_layout_fortran()) - - // stores the layout == order of the indices in memory - // wrapped into a little type to make constructor unambigous. - template - struct memory_layout { - ull_t value; - explicit memory_layout (ull_t v) : value(v) {assert((permutations::size(v)==Rank));} - explicit memory_layout (char ml='C') { - assert( (ml=='C') || (ml == 'F')); - value = (ml=='F' ? indexmaps::mem_layout::fortran_order(Rank) : indexmaps::mem_layout::c_order(Rank)); - } - memory_layout (memory_layout_fortran) { value = indexmaps::mem_layout::fortran_order(Rank); } - memory_layout (memory_layout_c) { value = indexmaps::mem_layout::c_order(Rank); } - template - explicit memory_layout(int i0, int i1, INT ... in) : value (permutations::permutation(i0,i1,in...)){ - static_assert( sizeof...(in)==Rank-2, "Error"); - } - memory_layout (const memory_layout & C) = default; - memory_layout (memory_layout && C) = default; - memory_layout & operator =( memory_layout const &) = default; - memory_layout & operator =( memory_layout && x) = default; - - bool operator ==( memory_layout const & ml) const { return value == ml.value;} - bool operator !=( memory_layout const & ml) const { return value != ml.value;} - - friend std::ostream &operator<<(std::ostream &out, memory_layout const &s) { - permutations::print(out, s.value); - return out; - - } - }; - - template memory_layout transpose(memory_layout ml, INT... is) { - static_assert(sizeof...(INT)==R, "!"); - return memory_layout{permutations::compose(ml.value, permutations::inverse(permutations::permutation(is...)))}; + explicit memory_layout() { + for (int u = 0; u < Rank; ++u) p[u] = u; } -}}//namespace triqs::arrays -#endif + + memory_layout(__mem_layout_fortran_order) { + order = memory_order_type::Fortran; + for (int u = 0; u < Rank; ++u) p[u] = Rank - u - 1; + } + + memory_layout(__mem_layout_c_order) : memory_layout() {} + + // internal use only : no check , we KNOW the order is not C or Fortran + template explicit memory_layout(mini_vector v, bool) : p(std::move(v)) { + order = memory_order_type::Custom; + } + + // For the user : we find out whether it is C or Fortran + template explicit memory_layout(mini_vector v) : p(std::move(v)) { + bool c = true, f = true; + for (int u = 0; u < Rank; ++u) { + c = c && (v[u] == u); + f = f && (v[u] == Rank - u - 1); + } + if (c) return; + order = (f ? memory_order_type::Fortran : memory_order_type::Custom); + } + + template + explicit memory_layout(int i0, INT... in) + : memory_layout(mini_vector{i0, in...}) {} + + bool operator==(memory_layout const &ml) const { return p == ml.p; } + bool operator!=(memory_layout const &ml) const { return !operator==(ml); } + + int operator[](int u) const { return p[u]; } + + bool is_c() const { return order == memory_order_type::C; } + bool is_fortran() const { return order == memory_order_type::Fortran; } + + mini_vector get_memory_positions() const { + auto r = mini_vector(utility::no_init_tag{}); + for (int u = 0; u < Rank; ++u) r[p[u]] = u; + return r; + } + + mini_vector get_indices_in_order() const { return p; } + + friend std::ostream &operator<<(std::ostream &out, memory_layout const &s) { return out << s.p; } + }; + + // ------------------------- functions ------------------------------------- + + /// Make a memory_layout from the indices + template memory_layout make_memory_layout(T... x) { return memory_layout(x...); } + + /// Make a memory_layout from the strides + template memory_layout memory_layout_from_strides(mini_vector const &strides) { + mini_vector c; // rely on init to 0 by default of mini_vector + for (int i = 0; i < Rank; ++i) + for (int j = i + 1; j < Rank; ++j) + if (strides[i] > strides[j]) + c[i]++; + else + c[j]++; + // we computed the map : index -> memory_rank, which is the inverse map, cf mem_layout + return memory_layout{c}; + } + + /// Apply a permutation to the indices + template memory_layout transpose(memory_layout const &ml, mini_vector const &perm) { + mini_vector res; + for (int u = 0; u < R; ++u) res[u] = perm[ml[u]]; + return memory_layout{res}; + } +} +} // namespace triqs::arrays diff --git a/triqs/arrays/indexmaps/cuboid/slice.hpp b/triqs/arrays/indexmaps/cuboid/slice.hpp index 30f37507..677a7102 100644 --- a/triqs/arrays/indexmaps/cuboid/slice.hpp +++ b/triqs/arrays/indexmaps/cuboid/slice.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,98 +18,118 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_CUBOID_SLICE_H -#define TRIQS_ARRAYS_INDEXMAP_CUBOID_SLICE_H +#pragma once #include -#include "./slice_traversal_order.hpp" -namespace triqs { namespace arrays { namespace indexmaps { - namespace cuboid_details { -#define LISILOSO l_type const * li, s_type const * si, l_type * lo, s_type * so - typedef size_t l_type; - typedef std::ptrdiff_t s_type; +namespace triqs { +namespace arrays { + namespace indexmaps { + namespace cuboid_details { +#define LISILOSO l_type const* li, s_type const* si, l_type* lo, s_type* so, int* imap + using l_type = size_t; + using s_type = std::ptrdiff_t; - template inline void _check_BC ( int N, int ind, size_t B, int ind_min =0) { } - template <> inline void _check_BC (int N, int ind, size_t B, int ind_min) { - if (!((ind >= ind_min) && (ind < int(B)))) TRIQS_ARRAYS_KEY_ERROR << " index "<= ind_min) && (ind < int(B)))) + TRIQS_ARRAYS_KEY_ERROR << " index " << N << " is out of domain: \n " << ind << " is not within [0," << B << "[\n"; +#endif + } - struct slice_calc { // group in a struct to avoid declaration order issue with the cross call of templates... + struct slice_calc { // group in a struct to avoid declaration order issue with the cross call of templates... - template - static void one_step(LISILOSO, std::ptrdiff_t& offset, size_t R){ - _check_BC (N, R, li[N]); - offset += R*si[N]; + template static void one_step(LISILOSO, std::ptrdiff_t& offset, size_t R) { + _check_BC(N, R, li[N]); + offset += R * si[N]; + imap[N] = -1; } - template - static void one_step(LISILOSO, std::ptrdiff_t& offset, range R){ - _check_BC (N, R.first(),li[N]); - lo[P] = ((R.last()==-1 ? li[N] : R.last()) - R.first() + R.step()-1 )/R.step(); // python behaviour - _check_BC (N, R.first() + (lo[P]!=0 ? (lo[P]-1) : 0)*R.step() ,li[N], -1); - so[P] = si[N] * R.step(); + template static void one_step(LISILOSO, std::ptrdiff_t& offset, range R) { + _check_BC(N, R.first(), li[N]); + lo[P] = ((R.last() == -1 ? li[N] : R.last()) - R.first() + R.step() - 1) / R.step(); // python behaviour + _check_BC(N, R.first() + (lo[P] != 0 ? (lo[P] - 1) : 0) * R.step(), li[N], -1); + so[P] = si[N] * R.step(); offset += R.first() * si[N]; + imap[N] = P; } - template - static typename std::enable_if<((EllipsisLength==1) || (!std::is_base_of::type::value)), void>::type - invoke(LISILOSO, s_type & offset, Arg0 const & arg0, Args const & ... args ) { - constexpr bool dP = (std::is_base_of::type::value ? 1 : 0); // Arg0 is range or ellipsis - one_step(li,si,lo,so,offset, arg0); - invoke(li,si,lo,so, offset, args...); + // ellipsis is a total range, we can simplify the computation in that case... + template static void one_step(LISILOSO, std::ptrdiff_t& offset, ellipsis) { + _check_BC(N, 0, li[N]); + lo[P] = li[N]; + _check_BC(N, (lo[P] != 0 ? (lo[P] - 1) : 0), li[N], -1); + so[P] = si[N]; + imap[N] = P; } - template - static typename std::enable_if<((EllipsisLength==0) && std::is_base_of::type::value), void>::type - invoke(LISILOSO, s_type & offset, Arg0 const & arg0, Args const & ... args ) { - invoke(li,si,lo,so, offset, args...); + template + static std14::enable_if_t<((EllipsisLength == 1) || (!std::is_base_of::type::value)), void> + invoke(LISILOSO, s_type& offset, Arg0 const& arg0, Args const&... args) { + constexpr bool dP = (std::is_base_of::type::value ? 1 : 0); // Arg0 is range or ellipsis + one_step(li, si, lo, so, imap, offset, arg0); + invoke(li, si, lo, so, imap, offset, args...); } - template - static typename std::enable_if<((EllipsisLength>1) && std::is_base_of::type::value), void>::type - invoke(LISILOSO, s_type & offset, Arg0 const & arg0, Args const & ... args ) { - //constexpr int dP = 1;//(std::is_base_of::type::value ? 1 : 0); // Arg0 is range or ellipsis - one_step(li,si,lo,so,offset, arg0); - invoke(li,si,lo,so, offset, arg0, args...); + template + static std14::enable_if_t<((EllipsisLength == 0) && std::is_base_of::type::value), void> + invoke(LISILOSO, s_type& offset, Arg0 const& arg0, Args const&... args) { + invoke(li, si, lo, so, imap, offset, args...); } - template static void invoke(LISILOSO, s_type & offset ) {} - }; + template + static std14::enable_if_t<((EllipsisLength > 1) && std::is_base_of::type::value), void> + invoke(LISILOSO, s_type& offset, Arg0 const& arg0, Args const&... args) { + one_step(li, si, lo, so, imap, offset, arg0); + invoke(li, si, lo, so, imap, offset, arg0, args...); + } + + template static void invoke(LISILOSO, s_type& offset) {} + }; #undef LISILOSO - }//namespace cuboid_details + } // namespace cuboid_details - // special case of no argument : - template struct slicer < cuboid::map > { typedef cuboid::map < R, Opt,To > r_type; }; - - // general case - template struct slicer < cuboid::map, Args...> { - - static const unsigned int len = sizeof...(Args); - static constexpr bool has_ellipsis = (count_type_occurrence::value>0); - static_assert((count_type_occurrence::value < 2), "Only one ellipsis is permitted"); - static_assert((len>=R || has_ellipsis), "Too few arguments in slice"); - static_assert((len<=R+(has_ellipsis ?1:0)), "Too many arguments in slice"); // + one to allow an empty ellipsis - - typedef cuboid::map im_t; - static constexpr int Rf = R - count_type_occurrence_not::value; - static constexpr ull_t To_i = im_t::traversal_order; - - // compute a new traversal order, only if there is no ellipsis - // the computation with ellipsis is not yet implemented... - static constexpr ull_t Tof = (has_ellipsis ? 0 : cuboid::slicing_TO_order::sliced_memory_order::value); - typedef cuboid::map < Rf , Opt, Tof> r_type; - -#ifndef TRIQS_WORKAROUND_INTEL_COMPILER_BUGS - static_assert(( (has_ellipsis) || (permutations::size(Tof) == Rf)), "Mismatch between Rank and the TraversalOrder"); -#endif - - static r_type invoke (im_t const & X, Args ... args) { - typename r_type::lengths_type newlengths; - typename r_type::strides_type newstrides; - std::ptrdiff_t newstart= X.start_shift(); - constexpr int EllipsisLength = R -len+1; - cuboid_details::slice_calc::invoke<0,0,flags::bound_check_trait::value,EllipsisLength>(&X.lengths()[0],&X.strides()[0],&newlengths[0],&newstrides[0],newstart, args...); - return r_type(std::move(newlengths),std::move(newstrides),newstart);// use move construction ? + // special case of no argument : + template struct slicer> { + using r_type = cuboid::map; }; - }; -}}}//namespaces triqs::arrays::indexmaps -#endif + + // general case + template struct slicer, Args...> { + + static const unsigned int len = sizeof...(Args); + static constexpr bool has_ellipsis = (count_type_occurrence::value > 0); + static_assert((count_type_occurrence::value < 2), "Only one ellipsis is permitted"); + static_assert((len >= R || has_ellipsis), "Too few arguments in slice"); + static_assert((len <= R + (has_ellipsis ? 1 : 0)), "Too many arguments in slice"); // + one to allow an empty ellipsis + + using im_t = cuboid::map; + static constexpr int Rf = R - count_type_occurrence_not::value; + // TO DO : compute the new traversal order using the same routine, constexpr, if needed ? + // Af the moment, using simply the C traversal order. + using Tof = void; //_traversal_c; + using r_type = cuboid::map; + + static r_type invoke(im_t const& X, Args... args) { + typename r_type::lengths_type newlengths; + typename r_type::strides_type newstrides; + std::ptrdiff_t newstart = X.start_shift(); + constexpr int EllipsisLength = R - len + 1; + mini_vector imap(utility::no_init_tag{}); + cuboid_details::slice_calc::invoke<0, 0, EllipsisLength>(&X.lengths()[0], &X.strides()[0], &newlengths[0], &newstrides[0], + imap.ptr(), newstart, args...); + + if (X.get_memory_layout().is_c()) return r_type(std::move(newlengths), std::move(newstrides), newstart, memory_layout{}); + if (X.get_memory_layout().is_fortran()) return r_type(std::move(newlengths), std::move(newstrides), newstart, memory_layout{FORTRAN_LAYOUT}); + + // Compute the new memory index order + mini_vector p(utility::no_init_tag{}); + for (int i = 0, j = 0; j < R; ++j) { + auto k = imap[X.get_memory_layout()[j]]; + if (k != -1) p[i++] = k; + } + + return r_type(std::move(newlengths), std::move(newstrides), newstart, memory_layout{p, bool()}); + }; + }; + } +} +} // namespaces triqs::arrays::indexmaps diff --git a/triqs/arrays/indexmaps/cuboid/slice_traversal_order.hpp b/triqs/arrays/indexmaps/cuboid/slice_traversal_order.hpp deleted file mode 100644 index b6b8835b..00000000 --- a/triqs/arrays/indexmaps/cuboid/slice_traversal_order.hpp +++ /dev/null @@ -1,48 +0,0 @@ -/******************************************************************************* - * - * TRIQS: a Toolbox for Research in Interacting Quantum Systems - * - * Copyright (C) 2012 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 . - * - ******************************************************************************/ -#ifndef TRIQS_ARRAYS_INDEXMAP_STORAGE_ORDER_H -#define TRIQS_ARRAYS_INDEXMAP_STORAGE_ORDER_H -#include "./mem_layout.hpp" -namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid { namespace slicing_TO_order { - - template struct _impl{ - static constexpr int n_range_ellipsis=0; - static constexpr ull_t mask(ull_t mo, int P, int c) { return 0;} - static constexpr ull_t smo (ull_t mo, ull_t ma, int P, int c) { return 0;} - }; - template struct _impl< A0,SliceArgs...> { - typedef _impl next_t; - static constexpr bool is_range = std::is_base_of::value; // range and ellipsis (derived from range) - static constexpr int n_range_ellipsis= next_t::n_range_ellipsis + is_range; - static constexpr ull_t mask(ull_t mo, int P=0, int c=0) { return ((is_range ? P+1 : 0) << (4*c)) + next_t::mask(mo,P + is_range,c+1);} - static constexpr ull_t r2 (ull_t mo, ull_t ma, int c) { return (ma >> (4*(permutations::apply(mo,c)))) & 0xF; } - static constexpr ull_t aux (ull_t mo,ull_t ma,ull_t r,int P,int c) { return (r!=0? ( (r-1) << (4*P)) : 0) + next_t::smo(mo,ma,P+(r!=0), c+1); } - static constexpr ull_t smo (ull_t mo, ull_t ma, int P=0, int c=0) { return aux(mo,ma,r2(mo,ma,c),P,c); } - }; - template - constexpr ull_t sliced_memory_order1(ull_t mo) { return _impl::n_range_ellipsis + 0x10*_impl::smo(mo, _impl::mask(mo));} - template - constexpr ull_t sliced_memory_order2(ull_t mo) { return (mem_layout::is_c (mo) ? c_order(_impl::n_range_ellipsis) : - ( mem_layout::is_fortran(mo) ? fortran_order(_impl::n_range_ellipsis) : sliced_memory_order1(mo) )); - } - template struct sliced_memory_order { static constexpr ull_t value = sliced_memory_order1(mo); }; -}}}}}//namespace triqs::arrays -#endif diff --git a/triqs/arrays/linalg/det_and_inverse.hpp b/triqs/arrays/linalg/det_and_inverse.hpp index 6f544b23..2f16cb7f 100644 --- a/triqs/arrays/linalg/det_and_inverse.hpp +++ b/triqs/arrays/linalg/det_and_inverse.hpp @@ -175,7 +175,7 @@ namespace arrays { ENABLE_IF(is_matrix_or_view::type>) triqs_arrays_assign_delegation(MT &lhs, inverse_lazy const &rhs) { static_assert(is_matrix_or_view::value, "Can only assign an inverse matrix to a matrix or a matrix_view"); - bool M_eq_inverse_M = ((lhs.indexmap().memory_indices_layout() == rhs.input().indexmap().memory_indices_layout()) && + bool M_eq_inverse_M = ((lhs.indexmap().get_memory_layout() == rhs.input().indexmap().get_memory_layout()) && (lhs.data_start() == rhs.input().data_start()) && (has_contiguous_data(lhs))); if (!M_eq_inverse_M) { lhs = rhs(); diff --git a/triqs/arrays/matrix.hpp b/triqs/arrays/matrix.hpp index 37d9a4c1..08732c74 100644 --- a/triqs/arrays/matrix.hpp +++ b/triqs/arrays/matrix.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,236 +18,257 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_MATRIX_H -#define TRIQS_ARRAYS_MATRIX_H +#pragma once #include "indexmaps/cuboid/map.hpp" #include "indexmaps/cuboid/slice.hpp" #include "impl/indexmap_storage_pair.hpp" #include "impl/assignment.hpp" #include "vector.hpp" -namespace triqs { namespace arrays { +namespace triqs { +namespace arrays { - template class matrix_view; - template class matrix; - - // ---------------------- matrix -------------------------------- - // -#define _IMPL_MATRIX_COMMON \ - bool is_square() const { return this->shape()[0] == this->shape()[1];}\ - \ - view_type transpose() const {\ - typename indexmap_type::lengths_type l; l[0] = this->indexmap().lengths()[1];l[1] = this->indexmap().lengths()[0];\ - typename indexmap_type::strides_type s; s[0] = this->indexmap().strides()[1];s[1] = this->indexmap().strides()[0];\ - return view_type( indexmap_type(l,s, this->indexmap().start_shift()), this->storage());\ - }\ - bool memory_layout_is_c() const { return this->indexmap().strides()[0] >= this->indexmap().strides()[1]; } \ + template class matrix_view; + template class matrix; + +// ---------------------- matrix -------------------------------- +// +#define _IMPL_MATRIX_COMMON \ + bool is_square() const { return this->shape()[0] == this->shape()[1]; } \ + \ + view_type transpose() const { \ + typename indexmap_type::lengths_type l; \ + l[0] = this->indexmap().lengths()[1]; \ + l[1] = this->indexmap().lengths()[0]; \ + typename indexmap_type::strides_type s; \ + s[0] = this->indexmap().strides()[1]; \ + s[1] = this->indexmap().strides()[0]; \ + return view_type(indexmap_type(l, s, this->indexmap().start_shift()), this->storage()); \ + } \ + bool memory_layout_is_c() const { return this->indexmap().strides()[0] >= this->indexmap().strides()[1]; } \ 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, Opt, TraversalOrder, IsConst, true, Tag::matrix_view > +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, TraversalOrder, IsConst, true, \ + Tag::matrix_view> - template - class matrix_view : Tag::matrix_view, TRIQS_CONCEPT_TAG_NAME(MutableMatrix), public IMPL_TYPE { - public : - typedef matrix regular_type; - typedef matrix_view view_type; - typedef matrix_view const_view_type; - typedef matrix_view weak_view_type; + template + class matrix_view : Tag::matrix_view, TRIQS_CONCEPT_TAG_NAME(MutableMatrix), public IMPL_TYPE { + public: + using regular_type = matrix; + using view_type = matrix_view; + using const_view_type = matrix_view; + using weak_view_type = matrix_view; + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using storage_type = typename IMPL_TYPE::storage_type; - typedef typename IMPL_TYPE::indexmap_type indexmap_type; - typedef typename IMPL_TYPE::storage_type storage_type; + /// Build from an IndexMap and a storage + template matrix_view(typename IMPL_TYPE::indexmap_type const& Ind, S const& Mem) : IMPL_TYPE(Ind, Mem) {} - /// Build from an IndexMap and a storage - template matrix_view (typename IMPL_TYPE::indexmap_type const & Ind,S const & Mem): IMPL_TYPE(Ind, Mem) {} - - /// Build from anything that has an indexmap and a storage compatible with this class - template matrix_view(const ISP & X): IMPL_TYPE(X.indexmap(),X.storage()) {} + /// Build from anything that has an indexmap and a storage compatible with this class + template matrix_view(const ISP& X) : IMPL_TYPE(X.indexmap(), X.storage()) {} #ifdef TRIQS_WITH_PYTHON_SUPPORT - /// Build from a numpy.array : throws if X is not a numpy.array - explicit matrix_view (PyObject * X): IMPL_TYPE(X, false, "matrix_view "){} + /// Build from a numpy.array : throws if X is not a numpy.array + explicit matrix_view(PyObject* X) : IMPL_TYPE(X, false, "matrix_view ") {} #endif - /// Copy construction - matrix_view( matrix_view const & X): IMPL_TYPE(X.indexmap(),X.storage()) {} + /// Copy construction + matrix_view(matrix_view const& X) : IMPL_TYPE(X.indexmap(), X.storage()) {} - matrix_view () = delete; + matrix_view() = delete; - // Move - matrix_view(matrix_view && X) { this->swap_me(X); } + // Move + matrix_view(matrix_view&& X) { this->swap_me(X); } - /// Swap - friend void swap( matrix_view & A, matrix_view & B) { A.swap_me(B);} + /// Swap + friend void swap(matrix_view& A, matrix_view& B) { A.swap_me(B); } - /// Rebind the view - void rebind(matrix_view const& X) { - this->indexmap_ = X.indexmap_; - this->storage_ = X.storage_; - } + /// Rebind the view + void rebind(matrix_view const& X) { + this->indexmap_ = X.indexmap_; + this->storage_ = X.storage_; + } - // rebind the other view, iif this is const, and the other is not. - template ENABLE_IFC(C) rebind(matrix_view const& X) { - this->indexmap_ = X.indexmap_; - this->storage_ = X.storage_; - } + // rebind the other view, iif this is const, and the other is not. + template ENABLE_IFC(C) rebind(matrix_view const& X) { + this->indexmap_ = X.indexmap_; + this->storage_ = X.storage_; + } - /** Assignement. The size of the array MUST match exactly. */ - template matrix_view & operator=(const RHS & X) {triqs_arrays_assign_delegation(*this,X); return *this; } + /** Assignement. The size of the array MUST match exactly. */ + template matrix_view& operator=(const RHS& X) { + triqs_arrays_assign_delegation(*this, X); + return *this; + } - matrix_view & operator=(matrix_view const & X) {triqs_arrays_assign_delegation(*this,X); return *this; }//cf array_view class comment + matrix_view& operator=(matrix_view const& X) { + triqs_arrays_assign_delegation(*this, X); + return *this; + } // cf array_view class comment - // Move assignment not defined : will use the copy = since view must copy data + // Move assignment not defined : will use the copy = since view must copy data - TRIQS_DEFINE_COMPOUND_OPERATORS(matrix_view); - _IMPL_MATRIX_COMMON; - }; + TRIQS_DEFINE_COMPOUND_OPERATORS(matrix_view); + _IMPL_MATRIX_COMMON; + }; //--------------------------------------------------------------------- // 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, bool Borrowed, bool IsConst> - struct ISPViewType< V, R, OptionFlags, To, Tag::matrix_view, Borrowed,IsConst> { - typedef typename std::conditional, - matrix_view>::type type; - }; + template + struct ISPViewType : std::conditional, + matrix_view> {}; #undef IMPL_TYPE - template - using matrix_const_view = matrix_view; + template + using matrix_const_view = matrix_view; - // ---------------------- matrix -------------------------------- -#define IMPL_TYPE indexmap_storage_pair < indexmaps::cuboid::map<2,Opt,TraversalOrder>, \ - storages::shared_block, Opt, TraversalOrder, false, false, Tag::matrix_view > +// ---------------------- matrix -------------------------------- +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, TraversalOrder, false, \ + false, Tag::matrix_view> - template - class matrix: Tag::matrix, TRIQS_CONCEPT_TAG_NAME(MutableMatrix), 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 matrix regular_type; - typedef matrix_view view_type; - typedef matrix_view const_view_type; - typedef matrix_view weak_view_type; + template + class matrix : Tag::matrix, TRIQS_CONCEPT_TAG_NAME(MutableMatrix), public IMPL_TYPE { + public: + using value_type = typename IMPL_TYPE::value_type; + using storage_type = typename IMPL_TYPE::storage_type; + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using regular_type = matrix; + using view_type = matrix_view; + using const_view_type = matrix_view; + using weak_view_type = matrix_view; - /// Empty matrix. - matrix(memory_layout<2> ml = memory_layout<2>(IMPL_TYPE::indexmap_type::traversal_order) ): IMPL_TYPE(indexmap_type(ml)) {} + /// Empty matrix. + matrix(memory_layout<2> ml = memory_layout<2>{}) : IMPL_TYPE(indexmap_type(ml)) {} - /// Move - explicit matrix(matrix && X) { this->swap_me(X); } + /// Move + explicit matrix(matrix&& X) { this->swap_me(X); } - /// - matrix(size_t dim1, size_t dim2, memory_layout<2> ml = memory_layout<2>(IMPL_TYPE::indexmap_type::traversal_order) ) : - IMPL_TYPE(indexmap_type(mini_vector(dim1,dim2),ml)) {} + /// + matrix(size_t dim1, size_t dim2, memory_layout<2> ml = memory_layout<2>{}) + : IMPL_TYPE(indexmap_type(mini_vector(dim1, dim2), ml)) {} - /// - matrix(mini_vector const & sha, memory_layout<2> ml = memory_layout<2>(IMPL_TYPE::indexmap_type::traversal_order)) : - IMPL_TYPE(indexmap_type(sha,ml)) {} + /// + matrix(mini_vector const& sha, memory_layout<2> ml = memory_layout<2>{}) : IMPL_TYPE(indexmap_type(sha, ml)) {} - /** Makes a true (deep) copy of the data. */ - matrix(const matrix & X): IMPL_TYPE(X.indexmap(),X.storage().clone()) {} + /** Makes a true (deep) copy of the data. */ + matrix(const matrix& X) : IMPL_TYPE(X.indexmap(), X.storage().clone()) {} - /// Build a new matrix from X.domain() and fill it with by evaluating X. X can be : - template - matrix(const T & X, TYPE_ENABLE_IF(memory_layout<2>, ImmutableCuboidArray) ml = memory_layout<2>(IMPL_TYPE::indexmap_type::traversal_order)): - IMPL_TYPE(indexmap_type(X.domain(),ml)) { triqs_arrays_assign_delegation(*this,X); } + /// Build a new matrix from X.domain() and fill it with by evaluating X. X can be : + template + matrix(const T& X, TYPE_ENABLE_IF(memory_layout<2>, ImmutableCuboidArray) ml = memory_layout<2>{}) + : IMPL_TYPE(indexmap_type(X.domain(), ml)) { + triqs_arrays_assign_delegation(*this, X); + } #ifdef TRIQS_WITH_PYTHON_SUPPORT - ///Build from a numpy.array X (or any object from which numpy can make a numpy.array). Makes a copy. - explicit matrix (PyObject * X): IMPL_TYPE(X, true, "matrix "){} + /// Build from a numpy.array X (or any object from which numpy can make a numpy.array). Makes a copy. + explicit matrix(PyObject* X) : IMPL_TYPE(X, true, "matrix ") {} #endif - // build from a init_list - template - matrix (std::initializer_list> const & l): - IMPL_TYPE(memory_layout<2>(IMPL_TYPE::indexmap_type::traversal_order)) { - size_t i=0,j=0; int s=-1; - for (auto const & l1 : l) { if (s==-1) s= l1.size(); else if (s != l1.size()) TRIQS_RUNTIME_ERROR << "initializer list not rectangular !";} - IMPL_TYPE::resize(typename IMPL_TYPE::domain_type (mini_vector(l.size(),s))); - for (auto const & l1 : l) { - for (auto const & x : l1) { (*this)(i,j++) = x;} - j=0; ++i; - } - } + // build from a init_list + template matrix(std::initializer_list> const& l) : IMPL_TYPE(memory_layout<2>()) { + size_t i = 0, j = 0; + int s = -1; + for (auto const& l1 : l) { + if (s == -1) + s = l1.size(); + else if (s != l1.size()) + TRIQS_RUNTIME_ERROR << "initializer list not rectangular !"; + } + IMPL_TYPE::resize(typename IMPL_TYPE::domain_type(mini_vector(l.size(), s))); + for (auto const& l1 : l) { + for (auto const& x : l1) { + (*this)(i, j++) = x; + } + j = 0; + ++i; + } + } - /** - * Resizes the matrix. NB : all references to the storage is invalidated. - * Does not initialize the matrix by default - */ - matrix & resize (size_t n1, size_t n2) { IMPL_TYPE::resize(typename IMPL_TYPE::domain_type (mini_vector(n1,n2))); return *this;} + /** + * Resizes the matrix. NB : all references to the storage is invalidated. + * Does not initialize the matrix by default + */ + matrix& resize(size_t n1, size_t n2) { + IMPL_TYPE::resize(typename IMPL_TYPE::domain_type(mini_vector(n1, n2))); + return *this; + } - /** - * Resizes the matrix. NB : all references to the storage is invalidated. - * Does not initialize the matrix by default - */ - matrix & resize (const indexmaps::cuboid::domain_t & l) { IMPL_TYPE::resize(l); return *this; } + /** + * Resizes the matrix. NB : all references to the storage is invalidated. + * Does not initialize the matrix by default + */ + matrix& resize(const indexmaps::cuboid::domain_t& l) { + IMPL_TYPE::resize(l); + return *this; + } - /// Assignement resizes the matrix. All references to the storage are therefore invalidated. - matrix & operator=(const matrix & X) { IMPL_TYPE::resize_and_clone_data(X); return *this; } + /// Assignement resizes the matrix. All references to the storage are therefore invalidated. + matrix& operator=(const matrix& X) { + IMPL_TYPE::resize_and_clone_data(X); + return *this; + } - /// Move assignment - matrix & operator=(matrix && X) { this->swap_me(X); return *this;} + /// Move assignment + matrix& operator=(matrix&& X) { + this->swap_me(X); + return *this; + } - /// Swap - friend void swap( matrix & A, matrix & B) { A.swap_me(B);} + /// Swap + friend void swap(matrix& A, matrix& B) { A.swap_me(B); } - /** - * Assignement resizes the matrix. All references to the storage are therefore invalidated. - * NB : to avoid that, do make_view(A) = X instead of A = X - */ - template - matrix & operator=(const RHS & X) { - static_assert( ImmutableCuboidArray::value, "Assignment : RHS not supported"); - IMPL_TYPE::resize(X.domain()); - triqs_arrays_assign_delegation(*this,X); - return *this; - } + /** + * Assignement resizes the matrix. All references to the storage are therefore invalidated. + * NB : to avoid that, do make_view(A) = X instead of A = X + */ + template matrix& operator=(const RHS& X) { + static_assert(ImmutableCuboidArray::value, "Assignment : RHS not supported"); + IMPL_TYPE::resize(X.domain()); + triqs_arrays_assign_delegation(*this, X); + return *this; + } - TRIQS_DEFINE_COMPOUND_OPERATORS(matrix); - _IMPL_MATRIX_COMMON; - };//matrix class + TRIQS_DEFINE_COMPOUND_OPERATORS(matrix); + _IMPL_MATRIX_COMMON; + }; // matrix class #undef _IMPL_MATRIX_COMMON #undef IMPL_TYPE - template - matrix make_unit_matrix(int dim) { + template matrix make_unit_matrix(int dim) { matrix r(dim, dim); r() = 1; return r; - } + } - template - matrix_view - make_matrix_view(ArrayType const & a) { - static_assert(ArrayType::rank ==2, "make_matrix_view only works for array of rank 2"); - return a; - } + template + matrix_view make_matrix_view(ArrayType const& a) { + static_assert(ArrayType::rank == 2, "make_matrix_view only works for array of rank 2"); + return a; + } - template - matrix - make_matrix(ArrayType const & a) { - static_assert(ArrayType::domain_type::rank ==2, "make_matrix only works for array of rank 2"); - return a; - } + template matrix make_matrix(ArrayType const& a) { + static_assert(ArrayType::domain_type::rank == 2, "make_matrix only works for array of rank 2"); + return a; + } - template - TYPE_ENABLE_IF(typename M::value_type, ImmutableMatrix) trace(M const &m) { - auto r = typename M::value_type{}; - if (first_dim(m) != second_dim(m)) - TRIQS_RUNTIME_ERROR << " Trace of a non square matrix"; - auto d = first_dim(m); - for (int i = 0; i < d; ++i) - r += m(i, i); - return r; - } -}}//namespace triqs::arrays + template TYPE_ENABLE_IF(typename M::value_type, ImmutableMatrix) trace(M const& m) { + auto r = typename M::value_type{}; + if (first_dim(m) != second_dim(m)) TRIQS_RUNTIME_ERROR << " Trace of a non square matrix"; + auto d = first_dim(m); + for (int i = 0; i < d; ++i) r += m(i, i); + return r; + } +} +} // namespace triqs::arrays // The std::swap is WRONG for a view because of the copy/move semantics of view. // Use swap instead (the correct one, found by ADL). -namespace std { - template - void swap( triqs::arrays::matrix_view & a , triqs::arrays::matrix_view & b)= delete; +namespace std { +template +void swap(triqs::arrays::matrix_view& a, triqs::arrays::matrix_view& b) = delete; } #include "./expression_template/matrix_algebra.hpp" -#endif diff --git a/triqs/arrays/matrix_tensor_proxy.hpp b/triqs/arrays/matrix_tensor_proxy.hpp index 0641e5e3..1fe9854b 100644 --- a/triqs/arrays/matrix_tensor_proxy.hpp +++ b/triqs/arrays/matrix_tensor_proxy.hpp @@ -44,6 +44,7 @@ namespace arrays { long n; typedef typename A_t::value_type value_type; + using traversal_order_t = typename A_t::traversal_order_t; typedef indexmaps::slicer slicer_t; typedef typename slicer_t::r_type indexmap_type; typedef typename indexmap_type::domain_type domain_type; @@ -98,6 +99,7 @@ namespace arrays { long n; typedef typename A_t::value_type value_type; + using traversal_order_t = typename A_t::traversal_order_t; typedef indexmaps::slicer slicer_t; typedef typename slicer_t::r_type indexmap_type; typedef typename indexmap_type::domain_type domain_type; diff --git a/triqs/arrays/python/numpy_extractor.cpp b/triqs/arrays/python/numpy_extractor.cpp index a99b9f3b..eedc3aaa 100644 --- a/triqs/arrays/python/numpy_extractor.cpp +++ b/triqs/arrays/python/numpy_extractor.cpp @@ -20,7 +20,6 @@ ******************************************************************************/ #include "../impl/common.hpp" #include "../impl/traits.hpp" -#include "../impl/flags.hpp" #include "../../utility/mini_vector.hpp" #include "./numpy_extractor.hpp" #ifdef TRIQS_WITH_PYTHON_SUPPORT diff --git a/triqs/arrays/storages/shared_block.hpp b/triqs/arrays/storages/shared_block.hpp index 9d899ec1..4c3987cf 100644 --- a/triqs/arrays/storages/shared_block.hpp +++ b/triqs/arrays/storages/shared_block.hpp @@ -49,12 +49,12 @@ namespace triqs { namespace arrays { namespace storages { 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) { 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; usize(); for (size_t u=0; u(arr,weak)) { data_ = sptr->p; s= sptr->size(); } diff --git a/triqs/arrays/vector.hpp b/triqs/arrays/vector.hpp index 13a010c2..c0e996f3 100644 --- a/triqs/arrays/vector.hpp +++ b/triqs/arrays/vector.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -18,197 +18,214 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAY_VECTOR_H -#define TRIQS_ARRAY_VECTOR_H +#pragma once #include "indexmaps/cuboid/map.hpp" #include "indexmaps/cuboid/slice.hpp" #include "impl/indexmap_storage_pair.hpp" #include "impl/assignment.hpp" -#include "impl/flags.hpp" #include -namespace triqs { namespace arrays { +namespace triqs { +namespace arrays { - template class vector_view; - template class vector; + template class vector_view; + template class vector; - // ---------------------- vector_view -------------------------------- +// ---------------------- vector_view -------------------------------- -#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<1,Opt,0> , storages::shared_block, Opt, 0, IsConst, true, Tag::vector_view > +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, void, IsConst, true, \ + Tag::vector_view> /** */ - template - class vector_view : Tag::vector_view, TRIQS_CONCEPT_TAG_NAME(MutableVector), public IMPL_TYPE { - public : - typedef vector regular_type; - typedef vector_view view_type; - typedef vector_view const_view_type; - typedef vector_view weak_view_type; + template + class vector_view : Tag::vector_view, TRIQS_CONCEPT_TAG_NAME(MutableVector), public IMPL_TYPE { + public: + using regular_type = vector; + using view_type = vector_view; + using const_view_type = vector_view; + using weak_view_type = vector_view; + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using storage_type = typename IMPL_TYPE::storage_type; - typedef typename IMPL_TYPE::indexmap_type indexmap_type; - typedef typename IMPL_TYPE::storage_type storage_type; + /// Build from an IndexMap and a storage + template vector_view(indexmaps::cuboid::map<1> const& Ind, S const& Mem) : IMPL_TYPE(Ind, Mem) {} - /// Build from an IndexMap and a storage - template vector_view (indexmaps::cuboid::map<1, Opt2,To2> const & Ind,S const & Mem): IMPL_TYPE(Ind, Mem) {} - - /// Build from anything that has an indexmap and a storage compatible with this class - template - vector_view(const ISP & X): IMPL_TYPE(X.indexmap(),X.storage()) {} + /// Build from anything that has an indexmap and a storage compatible with this class + template vector_view(const ISP& X) : IMPL_TYPE(X.indexmap(), X.storage()) {} #ifdef TRIQS_WITH_PYTHON_SUPPORT - /// Build from a numpy.array : throws if X is not a numpy.array - explicit vector_view (PyObject * X): IMPL_TYPE(X, false, "vector_view "){} + /// Build from a numpy.array : throws if X is not a numpy.array + explicit vector_view(PyObject* X) : IMPL_TYPE(X, false, "vector_view ") {} #endif - /// Copy construction - vector_view(vector_view const & X): IMPL_TYPE(X.indexmap(),X.storage()) {} + /// Copy construction + vector_view(vector_view const& X) : IMPL_TYPE(X.indexmap(), X.storage()) {} - vector_view () = delete; + vector_view() = delete; - // Move - vector_view(vector_view && X) { this->swap_me(X);} + // Move + vector_view(vector_view&& X) { this->swap_me(X); } - /// Swap - friend void swap( vector_view & A, vector_view & B) { A.swap_me(B);} + /// Swap + friend void swap(vector_view& A, vector_view& B) { A.swap_me(B); } - /// Rebind the view - void rebind(vector_view const& X) { - this->indexmap_ = X.indexmap_; - this->storage_ = X.storage_; - } + /// Rebind the view + void rebind(vector_view const& X) { + this->indexmap_ = X.indexmap_; + this->storage_ = X.storage_; + } - // rebind the other view, iif this is const, and the other is not. - template ENABLE_IFC(C) rebind(vector_view const& X) { - this->indexmap_ = X.indexmap_; - this->storage_ = X.storage_; - } + // rebind the other view, iif this is const, and the other is not. + template ENABLE_IFC(C) rebind(vector_view const& X) { + this->indexmap_ = X.indexmap_; + this->storage_ = X.storage_; + } - /** Assignment. The size of the array MUST match exactly. */ - template vector_view & operator=(const RHS & X) { triqs_arrays_assign_delegation(*this,X); return *this; } + /** Assignment. The size of the array MUST match exactly. */ + template vector_view& operator=(const RHS& X) { + triqs_arrays_assign_delegation(*this, X); + return *this; + } - vector_view & operator=(vector_view const & X) {triqs_arrays_assign_delegation(*this,X); return *this; }//cf array_view class comment + vector_view& operator=(vector_view const& X) { + triqs_arrays_assign_delegation(*this, X); + return *this; + } // cf array_view class comment - // Move assignment not defined : will use the copy = since view must copy data + // Move assignment not defined : will use the copy = since view must copy data - size_t size() const { return this->shape()[0];} + size_t size() const { return this->shape()[0]; } + std::ptrdiff_t stride() const { return this->indexmap().strides()[0]; } - std::ptrdiff_t stride() const { return this->indexmap().strides()[0];} + TRIQS_DEFINE_COMPOUND_OPERATORS(vector_view); - TRIQS_DEFINE_COMPOUND_OPERATORS(vector_view); - - // to make interface similar to std::vector : forward [] to () - template auto operator[](Arg && arg) const DECL_AND_RETURN((*this)(std::forward(arg))); - template auto operator[](Arg && arg) DECL_AND_RETURN((*this)(std::forward(arg))); - }; + // to make interface similar to std::vector : forward [] to () + template auto operator[](Arg&& arg) const DECL_AND_RETURN((*this)(std::forward(arg))); + template auto operator[](Arg&& arg)DECL_AND_RETURN((*this)(std::forward(arg))); + }; #undef IMPL_TYPE - template < class V, int R, ull_t OptionFlags, ull_t To, bool Borrowed, bool IsConst> - struct ISPViewType< V, R,OptionFlags,To, Tag::vector_view, Borrowed, IsConst> { typedef vector_view type; }; + template + struct ISPViewType { + using type = vector_view; + }; - template - using vector_const_view = vector_view; + template using vector_const_view = vector_view; - // ---------------------- vector-------------------------------- -#define IMPL_TYPE indexmap_storage_pair< indexmaps::cuboid::map<1,Opt,0> , storages::shared_block, Opt, 0, false, false,Tag::vector_view > +// ---------------------- vector-------------------------------- +#define IMPL_TYPE \ + indexmap_storage_pair, storages::shared_block, void, false, false, Tag::vector_view> - template - class vector: Tag::vector, TRIQS_CONCEPT_TAG_NAME(MutableVector), 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 vector regular_type; - typedef vector_view view_type; - typedef vector_view const_view_type; - typedef vector_view weak_view_type; + template class vector : Tag::vector, TRIQS_CONCEPT_TAG_NAME(MutableVector), public IMPL_TYPE { + public: + using value_type = typename IMPL_TYPE::value_type; + using storage_type = typename IMPL_TYPE::storage_type; + using indexmap_type = typename IMPL_TYPE::indexmap_type; + using regular_type = vector; + using view_type = vector_view; + using const_view_type = vector_view; + using weak_view_type = vector_view; - /// Empty vector. - vector(){} + /// Empty vector. + vector() {} - // Move - explicit vector(vector && X) { this->swap_me(X); } + // Move + explicit vector(vector&& X) { this->swap_me(X); } - /// - vector(size_t dim):IMPL_TYPE(indexmap_type(mini_vector(dim))) {} + /// + vector(size_t dim) : IMPL_TYPE(indexmap_type(mini_vector(dim))) {} - /// to mimic std vector - template - vector(size_t dim, Arg && arg):IMPL_TYPE(indexmap_type(mini_vector(dim))) { (*this)() = std::forward(arg);} + /// to mimic std vector + template vector(size_t dim, Arg&& arg) : IMPL_TYPE(indexmap_type(mini_vector(dim))) { + (*this)() = std::forward(arg); + } - /** Makes a true (deep) copy of the data. */ - vector(const vector & X): IMPL_TYPE(X.indexmap(),X.storage().clone()) {} + /** Makes a true (deep) copy of the data. */ + vector(const vector& X) : IMPL_TYPE(X.indexmap(), X.storage().clone()) {} - /** - * Build a new vector from X.domain() and fill it with by evaluating X. X can be : - * - another type of array, array_view, matrix,.... (any pair) - * - a expression : e.g. array > A( B+ 2*C); - * - ml : useless directly, since there only one ml, but used in generic code it maintains the same constructor as array, matrix - */ - template - //vector(const T & X, typename boost::enable_if< ImmutableCuboidArray >::type *dummy =0): - vector(const T & X, TYPE_ENABLE_IF(memory_layout<1>, ImmutableCuboidArray) ml = memory_layout<1>(IMPL_TYPE::indexmap_type::traversal_order)): - IMPL_TYPE(indexmap_type(X.domain(),ml)) { triqs_arrays_assign_delegation(*this,X); } + /** + * Build a new vector from X.domain() and fill it with by evaluating X. X can be : + * - another type of array, array_view, matrix,.... (any pair) + * - a expression : e.g. array > A( B+ 2*C); + * - ml : useless directly, since there only one ml, but used in generic code it maintains the same constructor as array, + * matrix + */ + template + // vector(const T & X, std14::enable_if_t< ImmutableCuboidArray ::value> *dummy =0): + vector(const T& X, TYPE_ENABLE_IF(memory_layout<1>, ImmutableCuboidArray) ml = memory_layout<1>{}) + : IMPL_TYPE(indexmap_type(X.domain(), ml)) { + triqs_arrays_assign_delegation(*this, X); + } #ifdef TRIQS_WITH_PYTHON_SUPPORT - ///Build from a numpy.array X (or any object from which numpy can make a numpy.array). Makes a copy. - explicit vector (PyObject * X): IMPL_TYPE(X, true, "vector "){} + /// Build from a numpy.array X (or any object from which numpy can make a numpy.array). Makes a copy. + explicit vector(PyObject* X) : IMPL_TYPE(X, true, "vector ") {} #endif - // build from a init_list - template - vector(std::initializer_list const & l): - IMPL_TYPE(indexmap_type(mini_vector(l.size()),memory_layout<1>(IMPL_TYPE::indexmap_type::traversal_order))) { - size_t i=0; - for (auto const & x : l) (*this)(i++) = x; - } + // build from a init_list + template + vector(std::initializer_list const& l) + : IMPL_TYPE(indexmap_type(mini_vector(l.size()), memory_layout<1>{})) { + size_t i = 0; + for (auto const& x : l) (*this)(i++) = x; + } - /** - * Resizes the vector. NB : all references to the storage is invalidated. - * Does not initialize the vector by default: to resize and init, do resize(IND).init() - */ - vector & resize (size_t L) { IMPL_TYPE::resize(typename IMPL_TYPE::domain_type(mini_vector(L))); return *this; } + /** + * Resizes the vector. NB : all references to the storage is invalidated. + * Does not initialize the vector by default: to resize and init, do resize(IND).init() + */ + vector& resize(size_t L) { + IMPL_TYPE::resize(typename IMPL_TYPE::domain_type(mini_vector(L))); + return *this; + } - /** - * Resizes the vector. NB : all references to the storage is invalidated. - * Does not initialize the vector by default: to resize and init, do resize(IND).init() - */ - vector & resize (const indexmaps::cuboid::domain_t & l) { IMPL_TYPE::resize(l); return *this; } + /** + * Resizes the vector. NB : all references to the storage is invalidated. + * Does not initialize the vector by default: to resize and init, do resize(IND).init() + */ + vector& resize(const indexmaps::cuboid::domain_t& l) { + IMPL_TYPE::resize(l); + return *this; + } - /// Assignement resizes the vector. All references to the storage are therefore invalidated. - vector & operator=(const vector & X) { IMPL_TYPE::resize_and_clone_data(X); return *this; } + /// Assignement resizes the vector. All references to the storage are therefore invalidated. + vector& operator=(const vector& X) { + IMPL_TYPE::resize_and_clone_data(X); + return *this; + } - /** - * Assignement resizes the vector. All references to the storage are therefore invalidated. - * NB : to avoid that, do make_view(A) = X instead of A = X - */ - template - vector & operator=(const RHS & X) { - static_assert( ImmutableCuboidArray::value, "Assignment : RHS not supported"); - IMPL_TYPE::resize(X.domain()); - triqs_arrays_assign_delegation(*this,X); - return *this; - } + /** + * Assignement resizes the vector. All references to the storage are therefore invalidated. + * NB : to avoid that, do make_view(A) = X instead of A = X + */ + template vector& operator=(const RHS& X) { + static_assert(ImmutableCuboidArray::value, "Assignment : RHS not supported"); + IMPL_TYPE::resize(X.domain()); + triqs_arrays_assign_delegation(*this, X); + return *this; + } - /// Move assignment - vector & operator=(vector && X) { this->swap_me(X); return *this;} + /// Move assignment + vector& operator=(vector&& X) { + this->swap_me(X); + return *this; + } - /// Swap - friend void swap( vector & A, vector & B) { A.swap_me(B);} + friend void swap(vector& A, vector& B) { A.swap_me(B); } - /// - size_t size() const { return this->shape()[0];} + size_t size() const { return this->shape()[0]; } + std::ptrdiff_t stride() const { return this->indexmap().strides()[0]; } - /// - std::ptrdiff_t stride() const { return this->indexmap().strides()[0];} + TRIQS_DEFINE_COMPOUND_OPERATORS(vector); - TRIQS_DEFINE_COMPOUND_OPERATORS(vector); + // to make interface similar to std::vector : forward [] to () + template auto operator[](Arg&& arg) const DECL_AND_RETURN((*this)(std::forward(arg))); + template auto operator[](Arg&& arg)DECL_AND_RETURN((*this)(std::forward(arg))); - // to make interface similar to std::vector : forward [] to () - template auto operator[](Arg && arg) const DECL_AND_RETURN((*this)(std::forward(arg))); - template auto operator[](Arg && arg) DECL_AND_RETURN((*this)(std::forward(arg))); - - };//vector class -}}//namespace triqs::arrays + }; // vector class +} +} // namespace triqs::arrays #undef IMPL_TYPE @@ -216,112 +233,121 @@ namespace triqs { namespace arrays { #include "./blas_lapack/copy.hpp" #include "./blas_lapack/swap.hpp" #include "./blas_lapack/axpy.hpp" -namespace triqs { namespace arrays { +namespace triqs { +namespace arrays { // norm2 squared - template typename std::enable_if::value, typename V::value_type>::type norm2_sqr(V const& a) { + template std14::enable_if_t::value, typename V::value_type> norm2_sqr(V const& a) { int dim = a.size(); auto r = typename V::value_type{}; for (int i = 0; i < dim; ++i) r += a(i) * a(i); return r; } - + // norm2 - template typename std::enable_if::value, typename V::value_type>::type norm2(V const& a) { + template std14::enable_if_t::value, typename V::value_type> norm2(V const& a) { using std::sqrt; return sqrt(norm2(a)); } // lexicographical comparison operators - template - typename std::enable_if< ImmutableVector::value && ImmutableVector::value , bool>::type - operator < (V1 const & a, V2 const & b) { return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());} + template + std14::enable_if_t::value&& ImmutableVector::value, bool> operator<(V1 const& a, V2 const& b) { + return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end()); + } - template - typename std::enable_if< ImmutableVector::value && ImmutableVector::value , bool>::type - operator > (V1 const & a, V2 const & b) { return (b + std14::enable_if_t::value&& ImmutableVector::value, bool> operator>(V1 const& a, V2 const& b) { + return (b < a); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_assign_delegation (vector & lhs, RHS const & rhs) { blas::copy(rhs,lhs); } + template + std14::enable_if_t::value> triqs_arrays_assign_delegation(vector& lhs, RHS const& rhs) { + blas::copy(rhs, lhs); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_compound_assign_delegation (vector & lhs, RHS const & rhs, char_<'A'>) { - T a = 1.0; blas::axpy(a,rhs,lhs); - } + template + std14::enable_if_t::value> triqs_arrays_compound_assign_delegation(vector& lhs, RHS const& rhs, + char_<'A'>) { + T a = 1.0; + blas::axpy(a, rhs, lhs); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_compound_assign_delegation (vector & lhs, RHS const & rhs, char_<'S'>) { - T a = -1.0; blas::axpy(a,rhs,lhs); - } + template + std14::enable_if_t::value> triqs_arrays_compound_assign_delegation(vector& lhs, RHS const& rhs, + char_<'S'>) { + T a = -1.0; + blas::axpy(a, rhs, lhs); + } - template - typename boost::enable_if< is_scalar_for > >::type - triqs_arrays_compound_assign_delegation (vector & lhs, RHS const & rhs, char_<'M'>) { - T a = rhs; blas::scal(a,lhs); - } + template + std14::enable_if_t>::value> triqs_arrays_compound_assign_delegation(vector& lhs, RHS const& rhs, + char_<'M'>) { + T a = rhs; + blas::scal(a, lhs); + } - template - typename boost::enable_if< is_scalar_for > >::type - triqs_arrays_compound_assign_delegation (vector & lhs, RHS const & rhs, char_<'D'>) { - T a = 1/rhs; blas::scal(a,lhs); - } + template + std14::enable_if_t>::value> triqs_arrays_compound_assign_delegation(vector& lhs, RHS const& rhs, + char_<'D'>) { + T a = 1 / rhs; + blas::scal(a, lhs); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_assign_delegation (vector_view & lhs, RHS const & rhs) { blas::copy(rhs,lhs); } + template + std14::enable_if_t::value> triqs_arrays_assign_delegation(vector_view& lhs, RHS const& rhs) { + blas::copy(rhs, lhs); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_compound_assign_delegation (vector_view & lhs, RHS const & rhs, char_<'A'>) { - T a = 1.0; blas::axpy(a,rhs,lhs); - } + template + std14::enable_if_t::value> triqs_arrays_compound_assign_delegation(vector_view& lhs, RHS const& rhs, + char_<'A'>) { + T a = 1.0; + blas::axpy(a, rhs, lhs); + } - template - typename boost::enable_if< is_vector_or_view >::type - triqs_arrays_compound_assign_delegation (vector_view & lhs, RHS const & rhs, char_<'S'>) { - T a = -1.0; blas::axpy(a,rhs,lhs); - } + template + std14::enable_if_t::value> triqs_arrays_compound_assign_delegation(vector_view& lhs, RHS const& rhs, + char_<'S'>) { + T a = -1.0; + blas::axpy(a, rhs, lhs); + } - template - typename boost::enable_if< is_scalar_for > >::type - triqs_arrays_compound_assign_delegation (vector_view & lhs, RHS const & rhs, char_<'M'>) { - T a = rhs; blas::scal(a,lhs); - } + template + std14::enable_if_t>::value> + triqs_arrays_compound_assign_delegation(vector_view& lhs, RHS const& rhs, char_<'M'>) { + T a = rhs; + blas::scal(a, lhs); + } - template - typename boost::enable_if< is_scalar_for > >::type - triqs_arrays_compound_assign_delegation (vector_view & lhs, RHS const & rhs, char_<'D'>) { - T a = 1/rhs; blas::scal(a,lhs); - } + template + std14::enable_if_t>::value> + triqs_arrays_compound_assign_delegation(vector_view& lhs, RHS const& rhs, char_<'D'>) { + T a = 1 / rhs; + blas::scal(a, lhs); + } - template - void triqs_arrays_assign_delegation(vector_view & av, std::vector const& vec) { - std::size_t size = vec.size(); - for(std::size_t n = 0; n < size; ++n) av(n) = vec[n]; - } + template + void triqs_arrays_assign_delegation(vector_view& av, std::vector const& vec) { + std::size_t size = vec.size(); + for (std::size_t n = 0; n < size; ++n) av(n) = vec[n]; + } // swapping 2 vector - template - void deep_swap(vector_view x, vector_view y) { - blas::swap(x,y); - } + template void deep_swap(vector_view x, vector_view y) { + blas::swap(x, y); + } - template - void deep_swap(vector & x, vector & y) { blas::swap(x,y);} - -}} + template void deep_swap(vector& x, vector& y) { blas::swap(x, y); } +} +} // The std::swap is WRONG for a view because of the copy/move semantics of view. // Use swap instead (the correct one, found by ADL). namespace std { -template -void swap(triqs::arrays::vector_view& a, triqs::arrays::vector_view& b) = delete; +template +void swap(triqs::arrays::vector_view& a, triqs::arrays::vector_view& b) = delete; } #include "./expression_template/vector_algebra.hpp" -#endif - diff --git a/triqs/gfs/data_proxies.hpp b/triqs/gfs/data_proxies.hpp index ebb32db9..e0314dde 100644 --- a/triqs/gfs/data_proxies.hpp +++ b/triqs/gfs/data_proxies.hpp @@ -25,7 +25,7 @@ #include #include "../arrays/matrix_tensor_proxy.hpp" -#define TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS +//#define TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS namespace triqs { namespace gfs { @@ -51,12 +51,12 @@ namespace triqs { namespace gfs { #ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS template auto operator()(S& data, long i) const DECL_AND_RETURN(data(i, arrays::ellipsis())); #else - auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); - auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); - auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); - auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); - auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); - auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(typename B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); + auto operator()(typename B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(typename B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); + auto operator()(typename B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(typename B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(typename B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); #endif }; @@ -68,12 +68,12 @@ namespace triqs { namespace gfs { template auto operator()(S & data, long i) const RETURN(make_matrix_view(data(i, arrays::ellipsis()))); #else /// The data access - auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); - auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); - auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); - auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); - auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); - auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(typename B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); + auto operator()(typename B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(typename B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); + auto operator()(typename B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(typename B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(typename B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); #endif }; diff --git a/triqs/utility/mini_vector.hpp b/triqs/utility/mini_vector.hpp index 816925e0..38cc11ef 100644 --- a/triqs/utility/mini_vector.hpp +++ b/triqs/utility/mini_vector.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by O. Parcollet + * Copyright (C) 2011-2014 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 @@ -22,13 +22,14 @@ #define TRIQS_ARRAYS_MINI_VECTOR_H #include "./first_include.hpp" #include +#include "./macros.hpp" +#include "./c14.hpp" #include "./compiler_details.hpp" #include "./exceptions.hpp" #include #include #include #include -#include #define TRIQS_MINI_VECTOR_NRANK_MAX 10 @@ -36,6 +37,8 @@ //#define TRIQS_MAKE_NVP(NAME,X) boost::serialization::make_nvp(NAME,X) namespace triqs { namespace utility { + struct no_init_tag{}; + template class mini_vector { T _data[(Rank!=0 ? Rank : 1)]; @@ -45,17 +48,107 @@ namespace triqs { namespace utility { public : typedef T value_type; - - mini_vector(){init();} -#define AUX(z,p,unused) _data[p] = x_##p; -#define IMPL(z, NN, unused) \ - mini_vector (BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(NN), T x_)){ \ - static_assert(Rank-1==NN,"mini_vector : incorrect number of variables in constructor");\ - BOOST_PP_REPEAT(BOOST_PP_INC(NN),AUX,nil) } - BOOST_PP_REPEAT(TRIQS_MINI_VECTOR_NRANK_MAX , IMPL, nil); -#undef IMPL -#undef AUX + mini_vector(){init();} + + mini_vector(no_init_tag){} + + mini_vector(T x_0) { + static_assert(Rank == 1, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + } + mini_vector(T x_0, T x_1) { + static_assert(Rank == 2, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + } + mini_vector(T x_0, T x_1, T x_2) { + static_assert(Rank == 3, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + } + mini_vector(T x_0, T x_1, T x_2, T x_3) { + static_assert(Rank == 4, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4) { + static_assert(Rank == 5, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4, T x_5) { + static_assert(Rank == 6, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + _data[5] = x_5; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4, T x_5, T x_6) { + static_assert(Rank == 7, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + _data[5] = x_5; + _data[6] = x_6; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4, T x_5, T x_6, T x_7) { + static_assert(Rank == 8, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + _data[5] = x_5; + _data[6] = x_6; + _data[7] = x_7; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4, T x_5, T x_6, T x_7, T x_8) { + static_assert(Rank == 9, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + _data[5] = x_5; + _data[6] = x_6; + _data[7] = x_7; + _data[8] = x_8; + } + mini_vector(T x_0, T x_1, T x_2, T x_3, T x_4, T x_5, T x_6, T x_7, T x_8, T x_9) { + static_assert(Rank == 10, "mini_vector : incorrect number of variables in constructor"); + _data[0] = x_0; + _data[1] = x_1; + _data[2] = x_2; + _data[3] = x_3; + _data[4] = x_4; + _data[5] = x_5; + _data[6] = x_6; + _data[7] = x_7; + _data[8] = x_8; + _data[9] = x_9; + }; + + /* + #define AUX(z,p,unused) _data[p] = x_##p; + #define IMPL(z, NN, unused) \ + mini_vector (BOOST_PP_ENUM_PARAMS(BOOST_PP_INC(NN), T x_)){ \ + static_assert(Rank-1==NN,"mini_vector : incorrect number of variables in constructor");\ + BOOST_PP_REPEAT(BOOST_PP_INC(NN),AUX,nil) } + BOOST_PP_REPEAT(TRIQS_MINI_VECTOR_NRANK_MAX , IMPL, nil); + #undef IMPL + #undef AUX + */ mini_vector(const mini_vector & x){ *this = x; } mini_vector(mini_vector && x){ *this = std::move(x); } @@ -77,7 +170,12 @@ namespace triqs { namespace utility { template mini_vector & operator=(const mini_vector & x){ for (int i=0;i class mini_vector { T _data[1]; public: + mini_vector(){} + mini_vector(no_init_tag){} + T & operator[](size_t i) { return _data[i];} const T & operator[](size_t i) const { return _data[i];} - }; + T * restrict ptr() { return _data;} + const T * restrict ptr() const { return _data;} + friend std::ostream & operator << ( std::ostream & out, mini_vector const & v ) {return out<<"[]";} + }; // ------------- Comparison -------------------------------------- @@ -161,13 +265,13 @@ namespace triqs { namespace utility { return res; } -#ifndef TRIQS_C11 - // ------------- transform a tuple into a minivector -------------------------------------- - - template mini_vector::value> tuple_to_mini_vector(TU const &t) { - return tuple::apply_construct_parenthesis::value>>(t); - } -#endif +//#ifndef TRIQS_C11 +// // ------------- transform a tuple into a minivector -------------------------------------- +// +// template mini_vector::value> tuple_to_mini_vector(TU const &t) { +// return tuple::apply_construct_parenthesis::value>>(t); +// } +//#endif }}//namespace triqs::arrays @@ -175,6 +279,8 @@ namespace std { // overload std::get to work with it template AUTO_DECL get(triqs::utility::mini_vector &v) RETURN(v[i]); template AUTO_DECL get(triqs::utility::mini_vector const &v) RETURN(v[i]); +template struct tuple_size>: std::integral_constant{}; + } #endif diff --git a/triqs/utility/mpi1.hpp b/triqs/utility/mpi1.hpp index f1803286..7283e6e6 100644 --- a/triqs/utility/mpi1.hpp +++ b/triqs/utility/mpi1.hpp @@ -154,11 +154,14 @@ namespace triqs { namespace mpi { template struct mpi_impl::ok && arrays::is_amv_value_or_view_class::value)> : boost_mpi_impl{}; // overload for views rvalues (created on the fly) - template - void reduce_in_place( communicator _c, arrays::array_view && a, int root =0) { reduce_in_place(_c,a,root);} + template void reduce_in_place(communicator _c, arrays::array_view&& a, int root = 0) { + reduce_in_place(_c, a, root); + } - template - void reduce( communicator _c, A const & a, arrays::array_view && b, int root =0) { reduce(_c,a,b,root);} + template + void reduce(communicator _c, A const& a, arrays::array_view&& b, int root = 0) { + reduce(_c, a, b, root); + } // to be implemented : scatter, gather for arrays diff --git a/triqs/utility/tuple_tools.hpp b/triqs/utility/tuple_tools.hpp index 999b0e0a..c2c05004 100644 --- a/triqs/utility/tuple_tools.hpp +++ b/triqs/utility/tuple_tools.hpp @@ -23,6 +23,7 @@ #include "./c14.hpp" #include #include +#include "./mini_vector.hpp" // adding to the std lib the reversed lazy tuple... // overloading & specializing only the functions needed here.