diff --git a/doc/reference/c++/arrays/assignment.rst b/doc/reference/c++/arrays/assignment.rst deleted file mode 100644 index 58851195..00000000 --- a/doc/reference/c++/arrays/assignment.rst +++ /dev/null @@ -1,57 +0,0 @@ -.. highlight:: c - -Assignment -========================= - -The `value classes` and the `view classes` have a quite general assignment operator. -We will illustrate it on the `array` class, it is the same for `matrix` and `vector`. - -* **Syntax** : the syntax is the same in both cases:: - - template array & operator=(const RHS & X); - template array_view & operator=(const RHS & X); - -* **What can be RHS ?** - - RHS can be ... anything that models the :ref:`ImmutableCuboidArray` concept ! - - e.g. : array, array_view, matrix, matrix_view, - but also formal expression (See , e.g. A+B), or any custom object of your choice. - -* **Behaviour** - - ================= ======================================================================= ====================================================================================== - Topic value class : array, matrix, vector view: array_view, matrix_view, vector_view - ================= ======================================================================= ====================================================================================== - RHS ... - RHS models the concept :ref:`ImmutableCuboidArray`. - RHS models :ref:`ImmutableCuboidArray` - [or less ? : RHS can be evaluated in the domain_type::value_type, no domain needed.]. - Effect - array is resized to have a domain X.domain() - View's domain must match X's domain or behaviour is undefined. - - array is filled with the evaluation of X. - array is filled with the evaluation of X. - ================= ======================================================================= ====================================================================================== - - * By evaluation of X, we mean : - - - a copy if X is an array. - - computing the value if X is an expression. - -Compound operators (+=, -=, * =, /=) -------------------------------------------------- - -* **Syntax** - - The syntax is natural :: - - template array & operator += (const RHS & X); - template array & operator -= (const RHS & X); - template array & operator *= (const Scalar & S); - template array & operator /= (const Scalar & S); - -* **What can be RHS ?** - - Same as for = operator. - -* **Behaviour** - - Similar to for = operator, with obvious changes. - - diff --git a/doc/reference/c++/arrays/concepts.rst b/doc/reference/c++/arrays/concepts.rst index a964efe2..56eb3a8c 100644 --- a/doc/reference/c++/arrays/concepts.rst +++ b/doc/reference/c++/arrays/concepts.rst @@ -18,17 +18,14 @@ i.e. whether one can or not modify a(i,j). * An `Immutable` array is just a pure function on the domain of definition. a(i,j) returns a int, or a int const &, that can not be modified (hence immutable). -* A `Mutable` array is on the other hand, a piece of memory, addressed by the indices, - which can be modified. a(i,j) can return a int &. - -The formal definition is given below - -.. note:: +* A `Mutable` array is an Immutable array, which can also be modified. Non const object return a reference, + e.g. a(i,j) can return a int &. Typically this is a piece of memory, with a integer coordinate system on it. - The tag (made by derivation) is use to quickly recognize that an object - models a given concept, due to the current lack of concept support in C++. - - It is used e.g. to detect to which algebra (Array or Matrix/Vector) the object belongs... +The main point here is that `Immutable` array is a much more general notion : +a formal expression made of array (E.g. A + 2*B) models this concept, but not the `Mutable` one. +Most algorithms only use the `Immutable` notion, when then are pure (mathematical) function +that returns something depending on the value of an object, without side effects. + .. _ImmutableCuboidArray: @@ -37,9 +34,9 @@ ImmutableCuboidArray * **Purpose** : - The most abstract definition of something that behaves like an immutable array. + The most abstract definition of something that behaves like an immutable array on a cuboid domain. - * it has a domain (hence a rank). + * it has a cuboid domain (hence a rank). * it can be evaluated on any value of the indices in the domain * NB : It does not need to be stored in memory. A formal expression, e.g. model this concept. @@ -220,7 +217,7 @@ Example : auto d = immutable_diagonal_matrix_view{a}; std::cout << "domain = " << d.domain()<< std::endl; std::cout << "d = "<< d << std::endl; - std::cout << "2*d = "<< matrix(2*d) << std::endl; + std::cout << "2*d = "<< make_matrix(2*d) << std::endl; std::cout << "d*d = "<< matrix(d*d) << std::endl; } diff --git a/doc/reference/c++/arrays/containers/call.rst b/doc/reference/c++/arrays/containers/call.rst index e43bb290..ae776596 100644 --- a/doc/reference/c++/arrays/containers/call.rst +++ b/doc/reference/c++/arrays/containers/call.rst @@ -92,6 +92,10 @@ Example } +.. toctree:: + :hidden: + + range_ell .. highlight:: c diff --git a/doc/reference/c++/arrays/containers/regular.rst b/doc/reference/c++/arrays/containers/regular.rst index ce19b875..072bdbba 100644 --- a/doc/reference/c++/arrays/containers/regular.rst +++ b/doc/reference/c++/arrays/containers/regular.rst @@ -52,7 +52,7 @@ NB: Rank is only present for array, since matrix have rank 2 and vector rank 1. :hidden: template_parameters - + Member types -------------------------------------- @@ -108,7 +108,7 @@ Member functions reg_constructors reg_assign - comp_ops + compound_ops call resize STL @@ -132,4 +132,4 @@ Non-member functions :hidden: stream - + swap diff --git a/doc/reference/c++/arrays/containers/views.rst b/doc/reference/c++/arrays/containers/views.rst index cd52a9e3..4d4a70c6 100644 --- a/doc/reference/c++/arrays/containers/views.rst +++ b/doc/reference/c++/arrays/containers/views.rst @@ -38,6 +38,12 @@ where triqs::ull_t is the type defined by : Views offers a strong guarantee of the existence of the data, like a std::shared_ptr. Cf :ref:`arr_view_memory` for details. + .. toctree:: + :hidden: + + memory + + Template parameters ---------------------------- diff --git a/doc/reference/c++/arrays/implementation_notes/centralconcept.rst b/doc/reference/c++/arrays/implementation_notes/centralconcept.rst deleted file mode 100644 index aacc28a6..00000000 --- a/doc/reference/c++/arrays/implementation_notes/centralconcept.rst +++ /dev/null @@ -1,39 +0,0 @@ -The central concept : HasImmutableArrayInterface -======================================================= - -The central concept used by the library (both in the colloquial and in the technical acceptance of the word) -is the notion of **Immutable Array Interface**. - -This is just a formal definition of what has an object has to implement ("which concept it has to model") -to look like an array. - -Basically, the answer is simple, an array is a (discrete) map between a domain of indices into an algebra -(typically R, ou C). Therefore it has to defined : - -* a domain of indices, which can be enumerated, e.g. cartesian products of integer ranges. -* an evaluation method [] so that a[ I ] has a `value_type` type for any index I in the domain. - -The precise concept if defined at :ref:`HasImmutableArrayInterface` concept. - -Examples : - -* a matrix has a domain of range(0,n1) x range(0,n2). -* a general array has a domain of dimension rank and values in R, C, other arrays... -* a triangular matrix has a more complex domain... -* any algebraic expression of arrays, matrices are like this (and those one are certainly NOT mutable objects !). - -* What is then a (partial) view/slice of an array ? - - Simply another map on the same data. - -* Do you want to see a 4 indices arrays as a matrix of double indices ? - - Simply change the map... - -* Why is this concept useful ? - - Because if you program another object of your own, which models it, - it is guaranteed to interact properly with arrays.... - - - diff --git a/doc/reference/c++/arrays/implementation_notes/contents.rst b/doc/reference/c++/arrays/implementation_notes/contents.rst index 2e87d181..5d8b7d5a 100644 --- a/doc/reference/c++/arrays/implementation_notes/contents.rst +++ b/doc/reference/c++/arrays/implementation_notes/contents.rst @@ -13,6 +13,7 @@ in a way suitable for studying/reading/maintening... :numbered: concepts + strategy storage indexmaps isp @@ -24,4 +25,6 @@ in a way suitable for studying/reading/maintening... hdf5 reinterpretation python + cuboid_formula + slicing diff --git a/doc/reference/c++/arrays/implementation_notes/expresssion_template.rst b/doc/reference/c++/arrays/implementation_notes/expression_template.rst similarity index 100% rename from doc/reference/c++/arrays/implementation_notes/expresssion_template.rst rename to doc/reference/c++/arrays/implementation_notes/expression_template.rst diff --git a/doc/reference/c++/arrays/slicing.rst b/doc/reference/c++/arrays/slicing.rst deleted file mode 100644 index 34ce1835..00000000 --- a/doc/reference/c++/arrays/slicing.rst +++ /dev/null @@ -1,107 +0,0 @@ -.. highlight:: c - -Partial views -================================== - -Various kind of partial views and slices can be made on arrays and matrices. - -* A `partial view` is defined as a view of a restricted portion of the array while - a `slice` is strictly speaking a partial view of a lower dimension of the original array, - e.g. a column of a matrix. - -* Partial views uses the ( ) operator, as the evaluation of the array:: - - array A(10,10); // defines an array - - A(1, range(0,2) ) // 1d slice - A(1, range()) // 1d slice taking all the second dim - - A(range(0,10,2), range(0,10,2)) // a 2d slice viewing every each elements with even coordinates. - - array_view SL = A(0,range(0,3)); // naming the view. No data copied here ! - array_view SL ( A(0,range(0,3))); // same thing ! - - auto SL = A(0,range(0,3)); // even simpler with C++11. - // CAREFUL : this is a weak view !!!! -> to be explained. - -* **Return type** : - - * Partial views of array or array_view return an array_view. - * Partial views of vector or vector_view return an vector_view. - * 2d partial views of matrix or matrix_view return matrix_view. - * BUT : (1d) slices of matrix or matrix_view return vector_view. - * 0d slices of anything are converted to the `value_type` of the array. - -Memory Weak view -^^^^^^^^^^^^^^^^^^^^^ - - The `range` type -^^^^^^^^^^^^^^^^^^^^^ - - `range` mimics the python `range`. It can be constructed with : - - * no argument : it then takes the whole set of indices in the dimension (like `:` in python) :: - - A(range(), 0) // take the first column of A - - * two arguments to specify a range :: - - A(range (0,3), 0) // means A(0,0), A(1,0), A(2,0) - - .. warning:: - the second element is excluded : range(0,3) is 0,1,2, like in Python. - - * three arguments : a range with a step :: - - A(range(0,4,2), 0) // means A(0,0), A(2,0) - - -The `ellipsis` type -^^^^^^^^^^^^^^^^^^^^^^ - -* Ellipsis can be provided in place of `range`, as in python. The type `ellipsis` is similar to range - except that it is implicitely repeated to as much as necessary. - -* Example: - - .. compileblock :: - - #include - using triqs::arrays::array; using triqs::arrays::ellipsis; - int main(){ - array B(2,3,4,5) ; - B(0,ellipsis(),3) ; // same as B(0, range(),range(), 3 ) - B(0,ellipsis(),2,3); // same as B(0, range(), 2, 3 ) - B(ellipsis(),2,3) ; // same as B( range(),range(), 2, 3 ) - } - - -* NB : there can be at most one ellipsis per expression (otherwise it would be meaningless). - -* Example of usage : - - Ellipsis are useful to write generic algorithms. For example, imagine that you want to sum - arrays on their first index : - - .. compileblock :: - - #include - using triqs::arrays::array; using triqs::arrays::ellipsis; - - // a generic function that sum array, array_view or in fact anything - // with the right concept on its first dimension - template - array sum0 (ArrayType const & A) { - array res = A(0,ellipsis()); - for (size_t u =1; u< first_dim(A); ++u) res += A(u,ellipsis()); - return res; - } - - // test - int main(){ - array A(5,2); A()=2; - array B(5,2,3); B() = 1; - std::cout<< sum0(A) << sum0(B) <` for these classes and functions. - - - diff --git a/doc/reference/c++/arrays_old/IO.rst b/doc/reference/c++/arrays_old/IO.rst deleted file mode 100644 index e07d3d36..00000000 --- a/doc/reference/c++/arrays_old/IO.rst +++ /dev/null @@ -1,44 +0,0 @@ -.. highlight:: c - -MPI and serialization -########################################## - -Serialization -============================ - -The `value classes` and the views are compatible with Boost Serialization library. - -MPI -============================ - -The `value classes` (array, matrix, vector) can be bcasted, reduced, -with the boost.mpi library. - - Example:: - - #include < - #include - #include - using namespace triqs::arrays; - namespace mpi=boost::mpi; - - int main() { - mpi::environment env(argc, argv); - mpi::communicator world; - - array A (2,2), B(2,2),C(2,2); - for (int i =0; i<2; ++i) - for (int j=0; j<2; ++j) - { A(i,j) = (1+world.rank())*(10*i+ j);} - - if (world.rank() ==0) cout<<" A = "< >(),0); - - int s= world.size(); - if (world.rank() ==0) cout<<" C = "<( (s*(s+1)/2) * A) < :: - - boost::python::object C_to_Py::convert::invoke (T const & x); // convert T into a python object - -* The opposite conversion consists in constructing a C++ object from the boost::python::object:: - - Py_to_C::convert::possible(boost::python::object x); // returns true iif the conversion of x into a T is possible - T Py_to_C::convert::invoke(boost::python::object x); // do it ! - -* T can be : - - * array, array_view, matrix, matrix_view, vector, vector_view - - * [this is not array lib but triqs] any std::vector, std::pair, boost::unordered_map of those, recursively. - ---> cf the doc of this converters - - -Copy vs View ------------------------ - - The python conversion of the array and array_view follows the policy of the C++ classes : - - * `value classes` (array, matrix, vector) **always** make copies : - - * when returning to python : they present a fresh copy of the data. - * when being contructed from python data, they make their own copy. - - * `view classes` **never** make copies : - - * when returning to python : they return a numpy, which is a view of their data. - By doing this, they transfer ownership of the data to the python interpreter. - Cf memorymanagement. - - * when being contructed from python data, they make view. - If this is not possible (e.g. the python object is not a numpy, but a list, the type are not exactly the same) - they throw an exception (`triqs::runtime_error`). - -Interfacing, wrapping ---------------------------- - -* boost.python - - These converters can be registered to build regular boost.python converter, using, e.g. :: - - triqs::python_tools::register_converter< triqs::arrays::array > (); - -* swig : not implemented [ to do : add the typemap] - - -Examples ------------------------ - - -Example 1 : wrapping a simple function with Boost Python -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -* The problem we want to solve : - -- write a simple function in C++ taking an array, and compute something from it. -- wrap this function and use it from python. - -* The function :: - - array_view f( array_view A) { - A(0,0) = 34; // do something - } - -* The wrapping code :: - - #include - #include - BOOST_PYTHON_MODULE(_testarray) { - using namespace boost::python; - triqs::python_tools::register_converter< triqs::arrays::array_view > (); - def("f", f); - } - -* Use it in a python code. - - .. code-block:: python - - import numpy,_testarray - a=numpy.array([[1.0,2],[3,4]]) - _testarray.f(a) - - -Example 2 : wrapping a simple function with Boost Python -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -* A variant: we can be more explicit in the type conversion - (in this case, we don't need to register the converter, since we do the job ourselves). - - -* The function :: - - boost::python::object f( boost::python::object M ) { - array_view A(M.ptr()); // make a view of M. Throws if impossible - A(0,0) = 100; // do something - return M; // return it - } - - -* The wrapping code :: - - #include - #include - BOOST_PYTHON_MODULE(_testarray) { - using namespace boost::python; - def("f", f); - } - -* Use it in a python code. - - .. code-block:: python - - import numpy,_testarray - a=numpy.array([[1.0,2],[3,4]]) - print _testarray.f(a) - - - - diff --git a/doc/reference/c++/arrays_old/Iterators.rst b/doc/reference/c++/arrays_old/Iterators.rst deleted file mode 100644 index 8ad8b176..00000000 --- a/doc/reference/c++/arrays_old/Iterators.rst +++ /dev/null @@ -1,2 +0,0 @@ -.. highlight:: c - diff --git a/doc/reference/c++/arrays_old/STL.rst b/doc/reference/c++/arrays_old/STL.rst deleted file mode 100644 index 86c38c78..00000000 --- a/doc/reference/c++/arrays_old/STL.rst +++ /dev/null @@ -1,92 +0,0 @@ -.. highlight:: c - -Compatibility with STL and iterators, Interaction with STL -################################################################## - -The arrays are compatible with STL containers and algorithms. - -Iterators -================ - - Standard iterators are provided that model the boost Mutable_ForwardIterator and ForwardIterator concepts - (and hence are STL compliant). - - The iterator implements also two additionnal methods : - - * it can be casted to a bool : it is true iif the iteration is not at end. - * it.indices() : returns const & to the indices at the point of the iteration. - - Examples:: - - array A (2,3); - for (auto it = A.begin(); it; ++it) *it =it.indices()[0] + 10 *it.indices()[1] ; - //for (array::iterator it = A.begin(); it; ++it) *it =it.indices()[0] + 10 *it.indices()[1] ; - -Indices generators -================================ - -The domain have an index generator that enumerates the indices in the domain. - -Examples:: - - array A (2,3); - - for (auto it = A.indexmap().domain().begin(); it; ++it) cout<<" "<<*it<::indexmap_type::domain_type, Permutations::permutation<0,1> > it_type; - for (it_type it(A.indexmap().domain()); it; ++it) cout<<" "<<*it<()); it; ++it) cout<<" "<<*it< A (2,3); - std::vector > VV; - VV.push_back(A); - -* ... or a map :: - - std::map > MAP; - MAP["1"] = A; - -* We can put a std::vector in an array ... :: - - std::vector V (10); - array B(V.size()), C(V.size()); - std::copy(V.begin(),V.end(),B.begin()); - -* ... or in reverse :: - - std::copy(B.begin(),B.end(),V.begin()); - -* ... or use other algorithms of std:: - - bool te(int x) { return (x<25);} - //.... - cout<<" max(B) "<< *std::max_element(B.begin(),B.end())< or - files. - -For example :: - - array A (2,2), B(2,2),C; - C= A + 2*B; - array D( A+ 2*B); - array F( 0.5 * A); // Type promotion is automatic - -The technique is called `expression templates`. It allows the elimination of temporaries, so that the clear -and readable code:: - - Z= A + 2*B + C/2; - -is in fact rewritten by the compiler into :: - - for (i,j,...) indices of A, B : - C(i,j) = A(i,j) + 2* B(i,j) + C(i,j)/2 - -instead of making a chain of temporaries (C/2, 2*B, 2*B + C/2...) that "ordinary" object-oriented programming would produce. -As a result, the produced code is as fast as if you were writing the loop yourself, -but with several advantages : - -* It is more **compact** and **readable** : you don't have to write the loop, and he indices range are computed automatically. -* It is much better for **optimization** : - - * What you want is to tell the compiler/library to compute this expression, not *how* to do it optimally on a given machine. - * For example, since the memory layout is decided at compile time, the library can traverse the data - in an optimal way, allowing machine-dependent optimization. - * The library can perform easy optimisations, e.g. for vector it will use blas if possible. - -Arrays vs matrices ----------------------- - -Because their multiplication is not the same, arrays and matrices don't form the same algebra. -Mixing them in expression would therefore be meaningless and it is therefore not allowed :: - - array A; - matrix M; - M + A; // --> ERROR - -However, you can always make a matrix_view from a array of rank 2 :: - - A + make_matrix_view(M); //--> OK. - -.. note:: - - Making view is very cheap, it only copies the index systems. - -Expressions are lazy ---------------------------- - -This means that constructing an expression is separated from evaluating it :: - - auto e = A + 2*B; // expression, purely formal, no computation is done - cout<< e < D(e); // now really makes the computation and store the result in D. - -The expression type is complicated (the expression in stored in the C++ type), so we used here -the C++0x `auto` to make simple things simple... - -FAQ ----------- - -Where can expressions be used ? -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Expressions models :ref:`HasImmutableArrayInterface` concept, so they can used *anywhere* -an object modeling this concept is accepted, e.g. : - -* array, matrix contruction -* operator =, +=, -=, ... - -They behave like an immutable array : they have a domain, they can be evaluated. -When `C` is assigned to the expression in the previous example, -the compiler just needs to compute the new domain for `C`, resize it and fill it by evaluation of the expression. - -What is the cost of this technique ? -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Thanks to the Boost Proto library, this can be done : - -* with an acceptable increase in compilation time (try it !). -* in about *2 pages of readable code* ! - diff --git a/doc/reference/c++/arrays_old/all.rst b/doc/reference/c++/arrays_old/all.rst deleted file mode 100644 index 3b6a4376..00000000 --- a/doc/reference/c++/arrays_old/all.rst +++ /dev/null @@ -1,361 +0,0 @@ -.. highlight:: c - -array and array_view -============================ - -array and array_view are the class for standard d-dimensional cuboid array -and the corresponding view. - -Template parameters ----------------------------- - -* The class has four template parameters (same for array_view). - - .. code-block:: c - - array - - - ============================ ========================== ======================================= - Template parameter Access in the class Meaning - ============================ ========================== ======================================= - ValueType value_type The type of the element of the array - Rank rank The rank of the array *[int]* - IndexOrderTag indexmap_type The ordering in memory : can be Tag::C, Tag::Fortran or a permutation - StorageTag storage_type The storage : Tag::shared_block or Tag::numpy - ============================ ========================== ======================================= - - -* Various IndexOrderTag are possible : - - ================= ==================================================================================================== - IndexOrderTag Meaning - ================= ==================================================================================================== - Tag::C C-style order *[default]* - Tag::Fortran Fortran-style order - P - P is a permutation - - Determined by a permutation P at compile time. Explain here the permutation, the convention. - ================= ==================================================================================================== - - -* Two possible storages : - - ================== ============================================================================ - StorageTag Meaning - ================== ============================================================================ - Tag::shared_block a (shared_ptr on a) C++ block *[default]* - Tag::numpy stored in a numpy array, in which case the array is also a numpy array - and read numpy, be returned as numpy, sliced into a numpy, etc... - ================== ============================================================================ - - -.. _array_constructors: - -Constructors ------------------ - -Intentionally, array and array_view have only a few constructors : - -========================================== =========================================================================================== -Constructors of array Comments -========================================== =========================================================================================== -array() - empty array of size 0 -array(size_t, ...., size_t) - from the dimensions -array(cuboid_domain const &) - a new array with the corresponding cuboid -array(const array &) - copy construction -array(const T & X) - Type T models the :ref:`HasImmutableArrayInterface` concept. - - X must have the appropriate domain (checked at compile time). - - Enabled iif has_immutable_array_interface::value == true. - - Constructs a new array of domain X.domain() and fills it with evaluation of X. -========================================== =========================================================================================== - -====================================================================== ======================================================================================= -Constructors of array_view Comments -====================================================================== ======================================================================================= -array_view(indexmap_type const & I, S_type const &) from a couple of indexmap I and storage of type S_type -array_view(const T & X) T is any type such that X.indexmap() and X.storage(). Models ISP ... -====================================================================== ======================================================================================= - -array_view are typically constructed by slicing (Cf below). - -* Examples :: - - array A(10,2); - array Af (2,2); - - //Higher dim, custom order : - 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); - - -Access to data, domain, simple evaluation, ... --------------------------------------------------------- - -array, array_view, matrix, matrix_view, vector, vector_view model HasImmutableArrayInterface. - - -Assignment & Copy --------------------- - -Every classes comes in two flavors: C and C_view (with C = array, matrix, vector, etc...). -These two flavors differ in the way they handle their data in construction, copy construction, assignement. -Basically, C owns its data, while C_view if only a view. - -array, matrix, vector -^^^^^^^^^^^^^^^^^^^^^^^^ - -They own their data. In many aspects, they are similar to like std::vector. - -* The data are contiguous in memory. -* Constructors and copy constructors all create a new memory block. If needed, they - make a *true* copy of the data. -* The assignment operator may create a new Storage if size do not match. -* As a result, /pointers to the data/ and reference to the storage are invalid after assignment. -* They can be resized, again invalidating all references. - -array_view, matrix_view, vector_view -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -These classes do not own their data, but only present a view of them. - -* The data may not be contiguous in memory (e.g. if the view is the result of a slice). -* Constructors only make another view of the data. -* They *never* copy data, so they are quite quick. - In particular, copy constructor makes shallow copy (i.e. return another view). -* The assignement operator just copy data into the view. Behaviour is undefined if the - size of the view is too small (define the macro ARRAY_CHECK for dynamical debugging checks). -* Pointers to data taken from the views are still valid after assignement. -* Views can be not be resized. - -.. warning:: **Memory management** - - Views carry a reference to the memory block they view, - which guarantees that memory will not be - dellocated before the destruction of the view. - Indeed, the Storage types implement incorporated a reference counting mechanism, - either using boost::shared_ptr for the C++ arrays, or using the python references - for the numpy storage. - The memory block will be dellocated when its array and all array_view - pointing to it or to a portion of it will be destroyed, and only at that moment. - -Examples:: - - array *A = new array (Matrix(2,3)); // create an array A - array_view B(*A); // making a view - delete A; // A is gone... - cout< array & operator=(const RHS & X); - - * RHS models HasImmutableArrayInterface. - * array is first resized to have a domain X.domain(), and then filled - with the evaluation of X (e.g. a copy if X is an array, computing the value if X is an expression). - - -array_view, matrix_view, vector_view -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -//.. cpp:function:: - template array_view & operator=(const RHS & X); - - * RHS models HasImmutableArrayInterface [ or less ? : RHS can be evaluated in the domain_type::value_type, no domain needed.]. - * Dimension of the view must match or behaviour is undefined. - -Iterators and interaction with STL containers and algorithms ----------------------------------------------------------------- - -STL compliant iterators, hence STL algorithms work... - -Examples:: - - array A (2,3); - - // first print the index generator - for (array::indexmap_type::domain_type::generator it = A.indexmap().domain().begin(); !it.at_end(); ++it) - cout<<" "<<*it<::iterator it = A.begin(); !it.at_end(); ++it) - { *it =it.indices().get<0>() + 10 *it.indices().get<1>() ; } - - int u=0; - for (array::iterator it = A.begin(); !it.at_end(); ++it,++u) { *it =u; } - - array A (2,3); - std::vector > VV; VV.push_back(A); - map > MAP; MAP["1"] = A; - - // Trying to put a vector in an array - std::vector V (10); - array B(V.size()), C(V.size()); - for (unsigned int i =0; i<10; ++i) V[i] = 10+i; - - std::copy(V.begin(),V.end(),B.begin()); - std::copy(B.begin(),B.end(),V.begin()); - cout<<" Number of elements <25 : "<< std::count_if(B.begin(), B.end(),te)< A (2,3); - array_view SL( A(range(0,2),0)); - array_view SL2( A(1,range(0,2))); - - -It is the standard way to produce a view. - -NB : - -* we use here the python convention: range(0,3) is 0:3, i.e. 0,1,2 NOT 0,1,2,3. -* Todo : in fact we should wrap the range to python::slice for interoperability with python. - -Serialization -------------------------------------------------- - -* Boost.serialization -* Boost.mpi - -Examples:: - - array A (2,2), B(2,2),C(2,2); - boost::mpi::reduce (world, A,C, std::plus >(),0); - - -* HDF5 (ALPS), eg. - -Examples:: - - array A (2,3),B,vc; - array Af,Bf,vf; - - alps::hdf5::oarchive ar1("data.h5"); - ar1 << alps::make_pvp("Tableau", A); - ar1 << alps::make_pvp("Tableau2", Af); - ar1 << alps::make_pvp("Tableau_view", A(range(),range(1,3))); - - alps::hdf5::iarchive ar2("data.h5"); - ar2 >> alps::make_pvp("Tableau", B); - ar2 >> alps::make_pvp("Tableau", Bf); - ar2 >> alps::make_pvp("Tableau_view", vc); - ar2 >> alps::make_pvp("TableauC",C); - - -blas/lapack interface -------------------------------------------------- - -* matrix, vector and their views are interfaced with blas/lapack, via boost::numerics::bindings. -* If needed (for a view), a temporary (and silent) copy is made to reorganize the -data before calling blas/lapack (impl: cache class). -Of course, performance is altered, but code is simple... - -Examples:: - - namespace blas = boost::numeric::bindings::blas; - namespace lapack = boost::numeric::bindings::lapack; - - triqs_arrays::vector > V(5),V2(5); - triqs_arrays::vector V3(2); - triqs_arrays::matrix M1(2,2), M2(2,2), M3(2,2); - - blas::axpy(2.0,V,V2); - blas::gemm(1.0,M1, M2, 1.0, M3); - blas::ger(1.0,V3,V3,M2); - - // invert - triqs_arrays::vector ipiv(2); - lapack::getrf(M1, ipiv); - lapack::getri(M1, ipiv); - -Transparent use of python arrays -------------------------------------------------- - -* If the storage is Tag::numpy, the memory block is allocated/viewed through the numpy interface. -* One can mix arrays with any storage in expression (they have the same concepts). -* boost python converters are enable for those arrays into numpy and their views [impl :broken for views]. - - -Expression -------------------------------------------------- - -Simple expressions are made using boost.proto. -Examples :: - - array A (2,2), B(2,2),C; - C= A + 2*B; - array D( A+ 2*B); - - // or even in C++0x : - auto e = A + 2*B; // expression, purely formal - array D(e); // really makes the computation - cout<< e < A (2,2), B(2,2),C(2,2); - C= A + 2*B; - C= std::plus >()(A,B); - C = A + Transpose(B); // Transpose(X) returns a lazy object that models HasImmutableArrayInterface. - C = A + Transpose(B + B); // X can also be an expression... - C = Transpose(B); // - array F( 0.5 * A); // Type promotion is automatic - - // non square - array R(2,3),Rt(3,2); - cout<<" R = "<< array(Transpose(R)) < M1(2,2), M2(2,2), M3; - for (int i =0; i<2; ++i) for (int j=0; j<2; ++j) { M1(i,j) = i+j; M2(i,j) = 1 + i -j ; } - - // The central instruction : note that matmul returns a lazy object - // that has ImmutableArray interface, and defines a specialized version assignment - // As a result this is equivalent to matmul_with_lapack(M1,M2,M3) : there is NO intermediate copy. - M3 = matmul(M1,M2); - - std::cout<<"M3 = "< > V(5),V2(5); - triqs_arrays::vector V3(2); - triqs_arrays::matrix M1(2,2), M2(2,2), M3(2,2); - - blas::axpy(2.0,V,V2); - blas::gemm(1.0,M1, M2, 1.0, M3); - blas::ger(1.0,V3,V3,M2); - - // invert - triqs_arrays::vector ipiv(2); - lapack::getrf(M1, ipiv); - lapack::getri(M1, ipiv); - diff --git a/doc/reference/c++/arrays_old/centralconcept.rst b/doc/reference/c++/arrays_old/centralconcept.rst deleted file mode 100644 index aacc28a6..00000000 --- a/doc/reference/c++/arrays_old/centralconcept.rst +++ /dev/null @@ -1,39 +0,0 @@ -The central concept : HasImmutableArrayInterface -======================================================= - -The central concept used by the library (both in the colloquial and in the technical acceptance of the word) -is the notion of **Immutable Array Interface**. - -This is just a formal definition of what has an object has to implement ("which concept it has to model") -to look like an array. - -Basically, the answer is simple, an array is a (discrete) map between a domain of indices into an algebra -(typically R, ou C). Therefore it has to defined : - -* a domain of indices, which can be enumerated, e.g. cartesian products of integer ranges. -* an evaluation method [] so that a[ I ] has a `value_type` type for any index I in the domain. - -The precise concept if defined at :ref:`HasImmutableArrayInterface` concept. - -Examples : - -* a matrix has a domain of range(0,n1) x range(0,n2). -* a general array has a domain of dimension rank and values in R, C, other arrays... -* a triangular matrix has a more complex domain... -* any algebraic expression of arrays, matrices are like this (and those one are certainly NOT mutable objects !). - -* What is then a (partial) view/slice of an array ? - - Simply another map on the same data. - -* Do you want to see a 4 indices arrays as a matrix of double indices ? - - Simply change the map... - -* Why is this concept useful ? - - Because if you program another object of your own, which models it, - it is guaranteed to interact properly with arrays.... - - - diff --git a/doc/reference/c++/arrays_old/compound_ops.rst b/doc/reference/c++/arrays_old/compound_ops.rst deleted file mode 100644 index a495d3e9..00000000 --- a/doc/reference/c++/arrays_old/compound_ops.rst +++ /dev/null @@ -1,21 +0,0 @@ -.. highlight:: c - -Compound operators (+=, -=, *=, /=) -====================================== - -The `value classes` and the `view classes` behaves similarly. -We will illustrate it on the `array` class, it is the same for `matrix` and `vector`. - - * **Syntax** - - The syntax is natural :: - - template array & operator += (const RHS & X); - template array & operator -= (const RHS & X); - template array & operator *= (const Scalar & S); - template array & operator /= (const Scalar & S); - - * **Behaviour** - - - Domain must match. - - X is evaluated and added term by term. diff --git a/doc/reference/c++/arrays_old/contents.rst b/doc/reference/c++/arrays_old/contents.rst deleted file mode 100644 index 8a5fd233..00000000 --- a/doc/reference/c++/arrays_old/contents.rst +++ /dev/null @@ -1,33 +0,0 @@ -Array library -**************** - -.. warning:: - This library is beta. - - This manual is not finished: work in progress. - - -.. highlight:: c - -.. toctree:: - :maxdepth: 1 - :numbered: - - introduction - getting_started - view_or_not_view - basic_classes - slicing - shape - debug - assignment - algebras - foreach - functional - STL - H5 - Interop_Python - IO - blas_lapack - FAQ - design/contents diff --git a/doc/reference/c++/arrays_old/debug.rst b/doc/reference/c++/arrays_old/debug.rst deleted file mode 100644 index 36891a98..00000000 --- a/doc/reference/c++/arrays_old/debug.rst +++ /dev/null @@ -1,48 +0,0 @@ - -.. highlight:: python - -.. _Debug: - -Bound checking and debug macros -=================================== - -To be fast, by default, no check are made on the indices while accessing elements or slicing. -However, checks can be activated in two ways : - -* Adding the `Tag::BoundCheck` option (in the options<...> template of the array, matrix, vector) - -* Defining the debug macro TRIQS_ARRAYS_ENFORCE_BOUNDCHECK, which switches the default option from `Tag::NoBoundCheck` to `Tag::BoundCheck` - for all arrays, matrices and vectors. - -In both cases, if the indices are not within the domain of defintion, an exception triqs::arrays::key_error -will be thrown. It's .what() returns the file and line where the exception occurs, with the stack of all in C++, -e.g. :: - - BOX>./bound_check_nopy - - catching key error in B - triqs::arrays key_error at ..../triqs_source/triqs/arrays/test/C++/./src/IndexMaps/cuboid/./cuboid_domain.hpp : 104 - - Trace is : - - ./bound_check_nopy void triqs::arrays::IndexMaps::cuboid_domain<2>::assert_key_in_domain > \ - (boost::tuples::tuple const&) const 0x181 [0x403e11] - ./bound_check_nopy main 0x916 [0x403016] - /lib/libc.so.6 __libc_start_main 0xfd [0x7f389e6abc4d] - ./bound_check_nopy [0x4023c9] - - key out of domain - key [1] = 3 is not within [0,3[ - -Further information on the line of the error in the stack can be retrieved with the addr2line utility `[linux only]` -which retrieve the source file and line from the address of the function (if you compile with -g). The script in TRIQS_SOURCE_DIR/dev_tools/analyse_stack.py -does it automatically for you (when compiled with -g) :: - - BOX> ./bound_check_nopy 2>&1 | analyse_stack.py - - .... - ./bound_check_nopy main 0x8a7 [0x403477] - ---> /home/parcolle/triqs_source/triqs/arrays/test/C++/bound_check.cpp:83 - - - diff --git a/doc/reference/c++/arrays_old/design/concepts.rst b/doc/reference/c++/arrays_old/design/concepts.rst deleted file mode 100644 index fc5733a5..00000000 --- a/doc/reference/c++/arrays_old/design/concepts.rst +++ /dev/null @@ -1,270 +0,0 @@ - -Concepts -============================================================= - -BoostSerializable -------------------------------------------------- - -* define serialize for boost serialize library - -Printable -------------------------------------------------- - - ====================================================== =========================================================== - Elements Comment - ====================================================== =========================================================== - std::ostream & operator << (std::ostream & out, ...) Printing - ====================================================== =========================================================== - -Domain -------------------------------------------------- - -* **Purpose** : The domain of definition of the array, i.e. the possible value of the indices. -* **Refines** : BoostSerializable, Printable -* **Definition** : - - =========================================== =========================================================== - Elements Comment - =========================================== =========================================================== - * static const unsigned int rank rank - * index_value_type type of the multi-index - * default constructor, copy construction - * operator = - * operator ==, != - * size_t number_of_elements() const number of elements - * generator type of the IndexGenerator that generates the indices - * begin() const/ end() const a generator at start/end - =========================================== =========================================================== - -* **Examples** : - - Typically, this is is a multi-index for an array, matrix, ...., e.g. - - * Cuboid : standard hyperrectangular arrays. This is little more than the tuple of the lengths. - * triangle, .... ? - -.. _HasImmutableArrayInterface: - -HasImmutableArrayInterface -------------------------------------------------- - -* **Purpose** : The most abstract definition of something that behaves like an immutable array. - * it has a domain (hence a rank). - * it can be evaluated on any value of the indices in the domain - * By combining the generator of the domain with the evaluation, it is therefore easy to - iterate on the values of such an object. - * NB : It does not need to be stored in memory. A formal expression, e.g. model this concept. - -* **Definition** ([A|B] denotes that the return type maybe A or B). - - ================================================================================================== ============================================================= - Elements Comment - ================================================================================================== ============================================================= - value_type Type of the element of the array - domain_type Type of the domain. - [domain_type const & | domain_type] domain() const Access to the domain. - [ value_type | value_type const &] operator[] (domain_type::index_value_type const &) const Evaluation. - ================================================================================================== ============================================================= - -* **Examples** : - * array, array_view, matrix, matrix_view, vector, vector_view (all implemented as Indexmap_Storage_Pair) - * array expressions (in which case the returns are not const & since they are computed, not stored). - -IndexGenerator -------------------------------------------------- -* **Purpose** : Generate the indices of a domain. - -* **Definition** : - - ============================================================== ================================================================================================== - Elements Comment - ============================================================== ================================================================================================== - domain_type Type of the domain whose indices are generated. - default contruction, copy construction - construction from (domain_type const &, bool atend=false) - IndexGenerator & operator=(const IndexGenerator & ) - operator ++ - operator = - operator ==, != - domain_type::index_value_type const & operator * () const Access to the value of the multi-index - bool at_end() const True iif the generator has reached the last value - ============================================================== ================================================================================================== - - -* **Examples** : - - * cuboid_index_generator - -IndexMap -------------------------------------------------- - -* **Purpose** : - - * Store the mapping of the index domain to a linear array. - * It is basically a function : indices --> 1d position, the bijection between the indices - of an element and its position in the memory block. -* **Refines** : BoostSerializable, Printable -* **Definition** : - - ======================================================================== ================================================================================================== - Elements Comment - ======================================================================== ================================================================================================== - * domain_type The type of the domain. - * domain_type const & domain() const The domain. - * default constructor, copy construction Cpy is a true copy. - * can be constructed from domain_type const & - * size_t operator[] (domain_type::index_value_type const & key ) const The mapping itself. - * iterator A type modeling IndexMapIterator, which is the optimal memory traversal. - NB : the order of indices is chosen for optimal traversal of memory, it - does not need to be "natural". - cuboid_map also provides a natural_iterator for that purpose. - ======================================================================== ================================================================================================== - -* The type also has to define two free functions and to specialize a template : - - ========================================================== ================================================================================================== - Elements Comment - ========================================================== ================================================================================================== - * bool compatible_for_assignment (M1, M2) Returns whether an array/matrix/vector with map M1 can be equated to a array/matrix/vector with - map M2 - * bool raw_copy_possible (M1, M2) Is the assignment of an array/matrix/vector with map M2 into an array/matrix/vector with map M1 - doable with raw copy - * struct indexmap_iterator_adapter< It, I > Metafunction : - - - I is the IndexMap class - - It any similar IndexMapIterator which returns (in ::type) the IndexMapIterator on I - with the same order traversal as It. - - Example : It is a IndexMapIterator on I1 stored in C order, I is in Fortran order, - the result will be an IndexMapIterator on I that presents the data of I in C order - This is used in copying array with different indexmaps. - ========================================================== ================================================================================================== - - -* **Examples** : - * cuboid_map : a map of the cuboid indices in a fixed order in memory. - -IndexMapIterator -------------------------------------------------- - -* **Purpose** : - * A basic iterator on an IndexMap which can be dereferenced into the shift of successive elements compared to the start of the memory block. - * These iterators are kept as simple as possible, so that it is easy to implement new indices maps and their iterators. - * NB : In particular, they are *not* necessary STL-compliant. The array_iterator class will - take such an iterator and a Storage and produce a true, STL compliant iterator on the array (iterator_adapter). - -* **Definition** : - - ========================================================== ================================================================================================== - Elements Comment - ========================================================== ================================================================================================== - indexmap_type The index_map on which the iterator is iterating - domain_type Type of the domain whose indices are generated. - default contruction, copy construction - construction from (domain_type const &, bool atend=false) - IndexMapIterator & operator=(const IndexMapIterator & ) - IndexMapIterator & operator ++ - operator ==, != - std::ptrdiff_t operator*() const Dereference as a shift from the beginning of the array - domain_type::index_value_type const & indices () const Access to the value of the multi-index at the iterator position - bool at_end() const True iif the generator has reached the last value (in practice quicker that it = XX.end()). - ========================================================== ================================================================================================== - -* **Example(s)** : - - * cuboid_map_iterator - -Storage -------------------------------------------------- - -* **Purpose** : - * The storage of the array in memory, e.g. plain C++ array, a numpy, etc... - * A Storage keeps the reference to the memory block where the array is stored. - * NB : This memory block can be typically shared between various arrays and views, - so the Storage is just a reference. The memory is deallocated only - when all storages referencing it has been destroyed. -* **Refines** : BoostSerializable -* **Definition** : - - ====================================================== ================================================================================================== - Elements Comment - ====================================================== ================================================================================================== - value_type Type of the element stored, e.g. int, const int, double, const double, ... - default construction Makes a storage of size 0 - copy construction a shallow copy (another reference to the same data). - the copy construction is possible from another storage of the same value_type - up to the const qualifier. - The construction of a storage with value_type=T from a storage with value_type const T - is forbidden at compile time. - void operator = (const STO &) A shallow copy of the reference to the data. - - clone() const Create a clone of the data. - const_clone() const Create a clone of the data with const value_type (e.g. int--> const int). - - void raw_copy_from(const STO & X) Copy all the data from X to *this. Behaviour undefined if sizes do not match. - - size_t size() const Number of elements in the storage - - value_type & operator[](size_t p) const Access to the data. Behaviour is undefined if empty()==true. - ====================================================== ================================================================================================== - - -* **Examples** : - * shared_block : basically a shared_ptr to a basic C++ array, dressed to model the Storage concept. - * numpy : a dressing of a python numpy object to model the Storage concept. - - -StorageOrder concept -------------------------------------------------- - -* **Purpose** : - - * Store the order of indices in memory. - * Can be fixed at compile time, or dynamically (not implemented). - -* **Refines** : BoostSerializable -* **Definition** : - - ====================================================== ================================================================================================== - Elements Comment - ====================================================== ================================================================================================== - size_t index_number(size_t i) - static unsigned int rank - default construction - copy construction - bool is_Fortran() const Is it Fortran-style ordering ? - bool is_C() const Is it C-style ordering ? - - ====================================================== ================================================================================================== - - -* The type also has to define the == operator : - - ========================================================== ================================================================================================== - Elements Comment - ========================================================== ================================================================================================== - Operator == Defined between any of the ordering. - ========================================================== ================================================================================================== - - -* **Examples** :: - - storage_order_C // canonical C ordering - - storage_order_fortran // canonical Fortran ordering - - storage_order_custom

// custom storage given by a permutation P - - * Details on the custom order : - - * P = [p_1,... p_r] : p_1 is the fastest index, p_r the slowest. - - * Example:: - - storage_order_C == storage_order_custom< Permutations::reverse_identity > - - storage_order_fortran == storage_order_custom< Permutations::identity > - - storage_order_custom > // the fastest index in memory is 0, then 1, then 2. - - diff --git a/doc/reference/c++/arrays_old/design/contents.rst b/doc/reference/c++/arrays_old/design/contents.rst deleted file mode 100644 index 80e6c14c..00000000 --- a/doc/reference/c++/arrays_old/design/contents.rst +++ /dev/null @@ -1,10 +0,0 @@ -Design & Implementation -========================= - -.. toctree:: - :maxdepth: 2 - - strategy - slicing - concepts - cuboid_formula diff --git a/doc/reference/c++/arrays_old/design/cuboid_formula.rst b/doc/reference/c++/arrays_old/design/cuboid_formula.rst deleted file mode 100644 index e2185c11..00000000 --- a/doc/reference/c++/arrays_old/design/cuboid_formula.rst +++ /dev/null @@ -1,60 +0,0 @@ - -.. _cuboid_formula: - -Cuboid formula -====================== - -* **Notations** : - - * R = rank - * index : (i0,...,i_{R-1}) - * Lengths : (L0, ... , L_{R-1}) - * Strides : (S0, ... , S_{R-1}) - * position = start_shift + \sum_{j=0}^{R-1} i_j S_j - -* **Strides for compact cuboid** : - - If :math:`\sigma(i)` is the i-th index in memory (0 the first, R-1 the last one), we have - - .. math:: - :nowrap: - - \begin{align*} - S_{\sigma(0)} &= 1\\ - S_{\sigma(1)} &= L_{\sigma(0)} S_{\sigma(0)} = L_{\sigma(0)}\\ - S_{\sigma(i)} &= L_{\sigma(i-1)} S_{\sigma(i-1)} = \prod_{j=i-1}^{0} L_{\sigma(j)} - \end{align*} - - -* **Slicing the cuboid** : - - - * Argument of the slice : (R0, R1, ..., R_{R-1}), with R_i = range : (start,end,step) - for the moment, if R_i = E_i, a number, consider it as a special range of length 1. - * new index (k0,...., k_{R-1}). - i_u = R_u.start + k_u * R_u.step - * Compute the new strides : Sn_j - - .. math:: - :nowrap: - - \def\offset{\text{offset}} - \begin{align*} - position &= \offset + \sum_{j=0}^{R-1} i_j S_j \\ - &= \offset + \sum_{j=0}^{R-1} (k_j* R_j.step + R_j.start) S_j \\ - &= (\offset + \sum_{j=0}^{R-1} (k_j* R_j.step + R_j.start) S_j ) + \sum_{j=0}^{R-1} k_j (S_j * R_j.start) - \\ - \offset_\text{new} &= \offset + \sum_{j=0}^{R-1} (k_j* R_j.step + R_j.start) S_j \\ - &= position( R_j.start) \\ - Sn_j &= S_j * R_j.start - \end{align*} - - * Now : for the cases R_i = E_i, just drop the index. - - - - - - - - diff --git a/doc/reference/c++/arrays_old/design/slicing.rst b/doc/reference/c++/arrays_old/design/slicing.rst deleted file mode 100644 index 6f21a8d4..00000000 --- a/doc/reference/c++/arrays_old/design/slicing.rst +++ /dev/null @@ -1,51 +0,0 @@ -.. _DesignSlicing: - -Slicing -============================================================= - -In this section, IM denotes some IndexMaps. - -* We refer here to any partial view of the arrays as "slicing", i.e. subarrays, true slices, etc... - A slicing is any (meta)-function that take an indexmap and return another one. - - * It is a meta-function that computes the type of the resulting indexmap - * It is a function that computes the resulting indexmap. - -* The array/matrix/vector classes and their views, when called with the () operator, will : - - * forward all the arguments and their types to IndexMaps::slice, to compute the new indexmap IM2. - * If IM2 is of dimension 0, return a value_type & or const & (it is a simple number, not an array). - * Otherwise : return a new view, made of IM2 and the same data as for the original object. - -* Possible slices are defined by the IndexMap type. - Slicing an class C with IndexMap I1 produces a class C2_view, with IndexMap I2, - i.e. a new sliced IndexMap on the same data. - -* **Examples** : - * array and array_view can be sliced : -`` - array A(10,10); : defines an array - A(1,2) : element access. - A ( 1, range(0,2) ) : 1d slice - A( range(0,10,2), range(0,10,2)) : a 2d slice viewing every each elements with even coordinates. -`` - * matrix, matrix_view when sliced, return vector_view or matrix_view. - -* One can be much more general : e.g. slicing the diagonal of a matrix, etc... - - - -* Implementation is entirely done in the IndexMaps classes, by specializing 2 templates - (basically by writing the function IM1-> IM2). - The rest is done by indexmap_storage_pair class, which will compute the correct view class - depending on the view class and IM2 (view_from_tag_I_S template). - -:: - //In namespace IndexMaps::result_of - template - struct slice< IM, ArgsTuple> { typedef IM2 type; }; - - //In namespace IndexMaps : - template - typename result_of::slice::type slice(IM const &, ArgsTuple args); - diff --git a/doc/reference/c++/arrays_old/design/strategy.rst b/doc/reference/c++/arrays_old/design/strategy.rst deleted file mode 100644 index d985e716..00000000 --- a/doc/reference/c++/arrays_old/design/strategy.rst +++ /dev/null @@ -1,62 +0,0 @@ -.. _Design: - -Strategy -============================================================= - -All the classes are a combination of a system of indices (called IndexMap I in the following) -and a physical storage S in the computer (a block of memory), denoted as an IndexMap_Storage_Pair (I,S). - -* I models the IndexMap concept [REF below]. - - * It is the bijection between the a set of indices and the position in the memory block. It can be though as a coordinate system on the (linear) memory block. - * Various types of indices are possible (only the first is implemented now). - - * cuboid (the standard hypercubic array, the only one currently implemented) - * triangular arrays - * band matrix - * multi-indices, with indices made of pair e.g. - -* S models the Storage concept [REF]. - - * It is a handle to the memory block containing the actual data. - * It can be e.g.: - - * a C++ shared pointer to a memory block. - * a reference to a numpy array. - -This design has several consequences : - -* **Interoperability** : classes are widely interoperable, e.g. one can add a array and a matrix (if dimensions are ok of course). - one can also add a python numpy and a C++ array without any further coding. - -* It is straighforward to construct a matrix_view from an array, since it is the same couple , - just interpreted differently. - -* It is easy to view a array as a matrix by gathering indices (properly ordered in memory) : - one just has to provide a new IndexMap I2 to see the same data. - [ --> very useful for vertex computation in many body...] - -* Slicing, or partial view is very natural : it is just a function on indexmaps : I--> I2, - independantly of any storage. - - -Quick guide through the implementation -============================================================= - -The implementation is divided into basically 4 parts : - -* Storages : implements two storages shared_block and numpy - -* IndexMaps : implements cuboid index map, its domain and iterators - -* impl/indexmap_storage_pair.hpp : the basic implementation class for all user class. - It is basically just a couple of an indexmap and a storage, with shallow copy. - It also forward the slices to the indexmap and construct the correct views. - -* upper level : - * user class : array, array_view, matrix, matrix_view, vector, vector_view - * expression.hpp : boost proto expression - * numpy_interface.hpp : helper to get numpy into array - * lapack/blas interface - * hdf5 support. - diff --git a/doc/reference/c++/arrays_old/expr_arith.rst b/doc/reference/c++/arrays_old/expr_arith.rst deleted file mode 100644 index 87e824e5..00000000 --- a/doc/reference/c++/arrays_old/expr_arith.rst +++ /dev/null @@ -1,104 +0,0 @@ -.. highlight:: c - -.. _arith_expression: - -Arithmetic Expressions -------------------------------------------------- - -* **Definition** : - - By `expression`, we mean here an object, typically representing a mathematical expression, - which models the :ref:`HasImmutableArrayInterface` concept. - -* **Use** : - - Expression can be used : - - as RHS (Right Hand Side) of various operators (=, +=, -=, ...) - - at any place where an `expression` is expected. - -* **How to make expression ?** - - Expressions can be build easily in various ways : - - - Using arithmetic operators, e.g. A + 2*B - Examples :: - - array A (2,2), B(2,2),C; - C= A + 2*B; - array D( A+ 2*B); - array F( 0.5 * A); // Type promotion is automatic - - // or even in C++0x : - auto e = A + 2*B; // expression, purely formal - array D(e); // really makes the computation - cout<< e < M1(2,2), M2(2,2), M3; - M3 = matmul(M1,M2); - - will do the following : - - matmul returns a lazy object modelling the :ref:`Expression` concept. - - a result this is compiled as some `matmul_with_lapack(M1,M2,M3)` : **there is NO intermediate copy**. - - - mat_vec_mul - - - Building custom expressions. The transpose example... - - Describe the concept, what to implement, etc.... - - - Transpose:: - - array A (2,2), B(2,2),C(2,2); - C= A + 2*B; - C = A + Transpose(B); // Transpose(X) returns a lazy object that models HasImmutableArrayInterface. - C = A + Transpose(B + B); // X can also be an expression... - C = Transpose(B); // - - // non square - array R(2,3),Rt(3,2); - cout<<" R = "<< array(Transpose(R)) < - using namespace triqs::arrays; - - h5::H5File file( "ess.h5", H5F_ACC_TRUNC ); - - // first create a 'big' array in the file (no temporary is made in memory, it uses directly the HDF5 API). - const size_t N = 100; - h5::array_proxy Pn( file, "Z", make_shape(N,2,2) ); - - // a slice... - array the_slice (2,2); the_slice() = 1; - - // write the large array slice by slice - for (int u = 0; u - using namespace triqs::arrays; - - h5::H5File file ("ess.h5", H5F_ACC_RDONLY ); - h5::array_proxy P(file, "A"); - - array a_slice ; - for (int u = 0; u` for this class. - diff --git a/doc/reference/c++/arrays_old/h5_rw.rst b/doc/reference/c++/arrays_old/h5_rw.rst deleted file mode 100644 index db02a8e1..00000000 --- a/doc/reference/c++/arrays_old/h5_rw.rst +++ /dev/null @@ -1,20 +0,0 @@ -.. highlight:: c - -Simple read/write operations of an array (or a view) -============================================================================ - -Given an array (or an array_view), the functions `h5::write` and `h5::read` write and read it to/from the file -or any subgroup thereof. For example : - -.. literalinclude:: examples_code/h5_rw.cpp - -Note that : - - * The data in the hdf5 file are stored in C order. - - * The data is reordered in Fortran (or custom) order if necessary when reading/writing. - - * It also works with matrix and vector (but the storage is the same in HDF5). - - * It also works with the corresponding views. TO BE ILLUSTRATED. - diff --git a/doc/reference/c++/arrays_old/h5_stack.rst b/doc/reference/c++/arrays_old/h5_stack.rst deleted file mode 100644 index af79cc96..00000000 --- a/doc/reference/c++/arrays_old/h5_stack.rst +++ /dev/null @@ -1,42 +0,0 @@ -.. highlight:: c - -h5::array_stack : stacking arrays or scalars -================================================================ - -h5::array_stack writes a sequences of arrays of the same shape (or of scalars) into an hdf5 array with one more dimension, unlimited in the stacking direction. - -It is typically used to store a Monte-Carlo data series for later analysis. - -* If the base of the stack is an array of rank R, the resulting hdf5 array will be of rank R+1. - -* If the base of the stack is a simple number (double, int, ...), R=0. - -* The syntax is simple : - - * The << operator piles up an array/scalar onto the stack. - * The ++ operator advances by one slice in the stack. - * The () operator returns a view on the current slice of the stack. - -* The stack is bufferized in memory (`bufsize` parameter), so that the file access does not happen too often. - -* NB : beware to complex numbers ---> REF TO COMPLEX - -Reference ------------- - -Here is the :arrays_doxy:`full C++ documentation` for this class. - - -Tutorial ------------ - -A simple example with a stack of double: - -.. literalinclude:: examples_code/h5_stack_ex_sca.cpp - -A simple example with a stack of array of rank 2 : - -.. literalinclude:: examples_code/h5_stack_ex.cpp - - - diff --git a/doc/reference/c++/arrays_old/introduction.rst b/doc/reference/c++/arrays_old/introduction.rst deleted file mode 100644 index 5bf1aa54..00000000 --- a/doc/reference/c++/arrays_old/introduction.rst +++ /dev/null @@ -1,67 +0,0 @@ -Motivation -=========================== - -.. highlight:: c - - -This library provides a multi-dimensionnal array library -for numerical computations with the following characteristics/goals : - -* **Simplicity of use**. - - Arrays must be as simple to use as in python (numpy) or fortran. - This library is designed to be used by physicists, not only by professionnal programmers, - We do *a lot* of array manipulations, and we want to maintain *readable* codes. - - Examples:: - - using namespace triqs::arrays; // All classes of this library are here. We will omit in the following - array A(4,2,3); // creating an array - A(0,1,2) // access an element - A(range(1,2), range(), 4) // a partial_view (or slice). - -* Various interoperable forms of arrays : - - * array : the general cuboid n-dimensional array, with **custom memory ordering** at compile time. - * matrix : a matrix class, with C or Fortran order at compile time. - * vector : a vector class. - -* **Genericity, abstraction and performance** : - - We want to have simple, readeable code, with the same (or better) speed than manually written low level code. - Most optimisations should be delegated to the library and the compiler. - - For example, we allow the use of **expression templates** for the algebra of arrays or matrices :: - - array A(4,2),B(4,2),C; - // fill A, B.... - C = 2* A + B/2 - - Expressions templates allows the elimination of temporaries and hence compact and readable code - with the speed of manually written codes. Thanks to boost::proto library, this can nowadays be achieved - while maintaining a good compilation time, and with a short and readable code. - -* **Compilation time** is (at far as tested) kept under control. - -* **Complete interoperability with python numpy arrays**. - - This library is used a lot with mixed C++/python codes. - It handles quick conversion between the C++ and python world, e.g. : - - * work on a view of a numpy, - * create a array in C++, and return it as a numpy. - * mix the various kind of arrays transparently in C++ expressions. - -* **HDF5** : simple interface to hdf5 library for an easy storing/retrieving into/from HDF5 files. - -* **STL compliance**. - - The arrays are compliant with STL (algorithms, std::vector >, iterators, map, etc...). - -* **Simple interface to Blas/Lapack** for matrices, vectors. - -* **MPI** : compatibility with boost::mpi interface. - - - - diff --git a/doc/reference/c++/arrays_old/matrix.rst b/doc/reference/c++/arrays_old/matrix.rst deleted file mode 100644 index 8df0214e..00000000 --- a/doc/reference/c++/arrays_old/matrix.rst +++ /dev/null @@ -1,101 +0,0 @@ -.. highlight:: c - -matrix and matrix_view -============================ - -matrix and matrix_view : - -* are the class for matrix and the corresponding view. -* allow **custom memory ordering** (C or Fortran) at compile time. -* model the :ref:`HasImmutableArrayInterface` concept. - -Template parameters ----------------------------- - -* The classes have two template parameters. - - .. code-block:: c - - template class matrix_view; - template class matrix; - - ============================ ================================== ========================== ==================================================================== - Template parameter Accepted type Access in the class Meaning - ============================ ================================== ========================== ==================================================================== - ValueType normally a scalar, but any default value_type The type of the element of the array - constructible type (?). - Opt Option::options< ...> opt_type Compile time options, see below. - ============================ ================================== ========================== ==================================================================== - -Options template parameters are described :ref:`here `. -The memory order must obviously be only C or Fortran (memory_order is fine, but useless here...). - -.. _matrix_constructors: - -Constructors ------------------ - -Intentionally, matrix and matrix_view have only a few constructors : - -========================================== =========================================================================================== -Constructors of matrix Comments -========================================== =========================================================================================== -matrix() - empty matrix of size 0 -matrix(size_t, size_t) - from the dimensions. Does NOT initialize the elements of matrix to 0 ! -matrix(const matrix &) - copy construction -matrix (PyObject * X) - Construct a new matrix from the data of the Python object X. - NB : X is a borrowed reference, matrix does affect its counting reference. - - it takes a **copy** of the data of X (or of numpy(X)) into C++. - - X is first transformed into a numpy by the python numpy lib itself - (as if you call numpy.array(X)) if X is not a numpy array or an array of the wrong type - (e.g. you construct an matrix from a numpy of int....). -matrix(const T & X) - Type T models the :ref:`HasImmutableArrayInterface` concept. - - X must have the appropriate domain dimension (checked at compile time). - - Constructs a new matrix of domain X.domain() and fills it with evaluation of X. -========================================== =========================================================================================== - -* Examples :: - - matrix Mc(2,2); - matrix Mf(2,2); - -.. warning:: - The constructor from the dimensions does **NOT** initialize the matrix to 0 - (because it may not be optimal). - If needed, do it explicitely by `a()=0`; - - - -Constructors of matrix_views ----------------------------------------------- - -Automatic construction -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -matrix_view are normally automatically constructed by making (partial) views, ref:`Slicing`, e.g. :: - - matrix A(2,2); - A(range(),2) ; // --> this makes a view... - A() ; // --> this makes a view over the full array - - -Explicit construction -^^^^^^^^^^^^^^^^^^^^^^^^ - -To explicitly make a view of an matrix, use make_view or () :: - - matrix A(2,2); - make_view(A) //-> a view... - make_view(A) = 13 ; // to assign e.g. - A() = 13; // idem - - - -====================================================================== =========================================================================================================== -Constructors of matrix_view Comments -====================================================================== =========================================================================================================== -matrix_view(const matrix_view &) - Copy construction (shallow copy) -matrix_view(const T & X) - `[Advanced]` T is any type such that X.indexmap() and X.storage() can be used to construct a view. -====================================================================== =========================================================================================================== - - diff --git a/doc/reference/c++/arrays_old/memory_management.rst b/doc/reference/c++/arrays_old/memory_management.rst deleted file mode 100644 index 46116973..00000000 --- a/doc/reference/c++/arrays_old/memory_management.rst +++ /dev/null @@ -1,40 +0,0 @@ -Memory management -======================== - -.. highlight:: c - -* `View classes` contains a **reference counting** system to the memory block they view. - - This guarantees that memory will not be dellocated before the destruction of the view. - - The memory block will be dellocated when its array and all array_view - pointing to it or to a portion of it will be destroyed, and only at that moment. - - Example:: - - array *p = new array (2,3); // create an array p - array_view B(*p); // making a view - delete p; // p is gone... - B(0,0) = 314; cout< res; - /// I compute res - return boost::python::object(array_view(res)); // Cf section for convertion... - } - - The interpreter then take the responsability of destroying the data when needed (meaning here, long after f has returned, - when the python object returned will be cleaned). - - diff --git a/doc/reference/c++/arrays_old/options.rst b/doc/reference/c++/arrays_old/options.rst deleted file mode 100644 index 7a677e89..00000000 --- a/doc/reference/c++/arrays_old/options.rst +++ /dev/null @@ -1,79 +0,0 @@ -.. highlight:: c - -.. _option_template: - -Options template parameters -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -The `Option::options` template can take several parameters that determine various aspects of `array` (resp. `matrix`, `vector`) at **compile time**. -These parameters can be given in any order in the `options` template, but only at most only (they all have a default, that can also be changed -using macros for maximum flexibility). -These options parameters determine : - -* the storage order of indices in memory at **compile time**. - - =============================== =========================================================================================== - Template parameter Meaning - =============================== =========================================================================================== - `Tag::C` `[default]` C-style order - `Tag::Fortran` Fortran-style order ( NB the base index is still 0, not 1 !) - memory_order The tuple of size rank determines the storage order of indices in memory - memory_order_p

P is a triqs::arrays::Permutations::permutation that determines the storage order in memory - P[0] is the fastest index, P[rank - 1] the slowest index - =============================== =========================================================================================== - - Example of custom order :: - - custom > // same as index_order_Fortran, the first index is the fastest - custom > // same as index_order_C, the last index is the fastest - custom > // if indices are (i,j,k) : index j is fastest, then k, then i - -* whether bound check are done while accessing array elements or slicing the arrays. - - =============================== ==================================================================================================== - Template parameter Meaning - =============================== ==================================================================================================== - `Tag::NoBoundCheck` `[default]` Do no bound check at all. This is the default, for speed ! - `Tag::BoundCheck` Do bound checks. Boundary errors throw triqs::arrays::key_error exception (See ??) - =============================== ==================================================================================================== - - The macro `TRIQS_ARRAYS_ENFORCE_BOUNDCHECK` change the default to `Tag::BoundCheck` (for debugging purposes). - - Rational : one may want for some array to activate bound check and catch the exception, while keeping faster unchecked arrays elsewhere. - -* whether arrays are initialized at construction. - - =============================== ==================================================================================================== - Template parameter Meaning - =============================== ==================================================================================================== - `Tag::no_init` `[default]` No initialization of the array at construction (fastest). - `Tag::default_init` Initialize with the default constructor of ValueType - `Tag::nan_inf_init` Initialize to nan for floating type and std::limits::infinity() for integer type - as in D e.g. mainly for detecting easily lack of initialization errors. - =============================== ==================================================================================================== - - Note that : - - * The macro `TRIQS_ARRAYS_ENFORCE_INIT_NAN_INF` change the default to `Tag::nan_inf_init`. - * The macro `TRIQS_ARRAYS_ENFORCE_INIT_DEFAULT` change the default to `Tag::default_init`. - -* Several simple aliases are defined in the namespace Option for the most current cases : - - =============================== =============================== - Option aliases - =============================== =============================== - Option::C Option::options<`Tag::C`> - Option::Default Option::options<`Tag::C`> - Option::Fortran Option::options<`Tag::Fortran>` - =============================== =============================== - - -* Examples :: - - array > > A0 (2,3,4); - array > > A1 (2,3,4); - array > > A2 (2,3,4); - array A4 (2,3,4); - - - diff --git a/doc/reference/c++/arrays_old/python_tools.rst b/doc/reference/c++/arrays_old/python_tools.rst deleted file mode 100644 index a3055493..00000000 --- a/doc/reference/c++/arrays_old/python_tools.rst +++ /dev/null @@ -1,18 +0,0 @@ -.. highlight:: c - -:: - - namespace Py_to_C { - template struct convert{ - static bool possible (boost::python::object x); // returns true iif the conversion of x into a T is possible - static T invoke(boost::python::object x); // do it ! - }; - } - - namespace C_to_Py { - template struct convert { - static boost::python::object invoke (T const & x); // convert T into a python object - }; - } - - diff --git a/doc/reference/c++/arrays_old/shape.rst b/doc/reference/c++/arrays_old/shape.rst deleted file mode 100644 index b4f6b074..00000000 --- a/doc/reference/c++/arrays_old/shape.rst +++ /dev/null @@ -1,33 +0,0 @@ - -.. _Shape: - -Shape, resize -================================== - -Shape --------------------- - -* array, matrix and vector have a method shape() that returns a `shape_type` object - i.e. a mini_vector - -* Tutorial:: - - array A(2,3); - A.shape()[0] == 2; - A.shape()[1] == 3; - -Resize --------- - -* The `value class` array, matrix and vector can be resized:: - - array A(2,3); - A.resize ( make_shape (5,5) ) - - matrix M; - M.resize( 3,3); - - vector V; - V.resize(10); - - diff --git a/doc/reference/c++/arrays_old/slicing.rst b/doc/reference/c++/arrays_old/slicing.rst deleted file mode 100644 index fc8cfbfb..00000000 --- a/doc/reference/c++/arrays_old/slicing.rst +++ /dev/null @@ -1,101 +0,0 @@ -.. highlight:: c - -.. _Slicing: - -Element access & partial views -================================== - -Element access --------------------- - -* The elements of the arrays can be naturally accessed using the ( ) operator:: - - array A(10,10); // defines an array - A(1,2) =10 // element access. - -Partial view, slices ------------------------- - -Various kind of partial views and slices can be made on arrays and matrices. - -* A `partial view` is defined as a view of a restricted portion of the array while - a `slice` is strictly speaking a partial view of a lower dimension of the original array, - e.g. a column of a matrix. - -* Partial views uses the ( ) operator, as the evaluation of the array:: - - array A(10,10); // defines an array - - A ( 1, range(0,2) ) // 1d slice - A ( 1, range() ) // 1d slice taking all the second dim - - A( range(0,10,2), range(0,10,2)) // a 2d slice viewing every each elements with even coordinates. - - array_view SL( A(0,range(0,3))); // naming the view - auto SL = A(0,range(0,3)); // C++0x with auto [check this] - -* **Return type** : - - * Partial views of array or array_view return an array_view. - * Partial views of vector or vector_view return an vector_view. - * 2d partial views of matrix or matrix_view return matrix_view. - * BUT : (1d) slices of matrix or matrix_view return vector_view. - * 0d slices of anything are converted to the `value_type` of the array. - -The `range` type -^^^^^^^^^^^^^^^^^^^^^ - - `range` mimics the python `range`. It can be constructed with : - - * no argument : it then takes the whole set of indices in the dimension (like `:` in python) :: - - A(range(), 0) // take the first column of A - - * two arguments to specify a range :: - - A(range (0,3), 0) // means A(0,0), A(1,0), A(2,0) - - .. warning:: - the second element is excluded : range(0,3) is 0,1,2, like in Python. - - * three arguments : a range with a step :: - - A(range(0,4,2), 0) // means A(0,0), A(2,0) - - -The `ellipsis` type -^^^^^^^^^^^^^^^^^^^^^^ - -* Ellipsis can be provided in place of `range`, as in python. The type `ellipsis` is similar to range - except that it is implicitely repeated to as much as necessary. - -* Example :: - - array B(2,3,4,5) ; - - B(0, ellipsis(), 3 ) ; // same as B(0, range(),range(), 3 ) - B(0, ellipsis(),2, 3 ) ; // same as B(0, range(), 2, 3 ) - B( ellipsis(),2, 3 ) ; // same as B( range(),range(), 2, 3 ) - -* NB : there can be at most one ellipsis per expression (otherwise it would be meaning less). - -* Example of usage : - - Ellipsis are useful to write generic algorithms. For example, imagine that you want to sum - arrays on their first index :: - - template - array_view sum0 (ArrayType const & A) { - array res = A(0,ellipsis()); - for (size_t u =1; u< A.shape()[0]; ++u) res += A(u,ellipsis()); - return res; - } - - ///..... - - array A(5,2); - array B(5,2,3); - std::cout<< sum0(A) << sum0(B) < class vector_view; - template class vector; - - ============================ ================================== ========================== ==================================================================== - Template parameter Accepted type Access in the class Meaning - ============================ ================================== ========================== ==================================================================== - ValueType normally a scalar, but any default value_type The type of the element of the vector - constructible type (?). - Opt Option::options< ...> opt_type Compile time options, see below. - ============================ ================================== ========================== ==================================================================== - -Options template parameters are described :ref:`here `. -The memory order is obviously useless here... - -.. _vector_constructors: - -Constructors ------------------ - -Intentionally, vector and vector_view have only a few constructors : - -========================================== =========================================================================================== -Constructors of vector Comments -========================================== =========================================================================================== -vector() - empty vector of size 0 -vector(size_t) - from the dimension. Does **NOT** initialize the elements of the vector to 0 ! -vector(const vector &) - copy construction -vector (PyObject * X) - Construct a new vector from the data of the Python object X. - NB : X is a borrowed reference, vector does affect its counting reference. - - it takes a **copy** of the data of X (or of numpy(X)) into C++. - - X is first transformed into a numpy by the python numpy lib itself - (as if you call numpy.array(X)) if X is not a numpy array or an array of the wrong type - (e.g. you construct an vector from a numpy of int....). -vector(const T & X) - Type T models the :ref:`HasImmutableArrayInterface` concept. - - X must have the appropriate domain (checked at compile time). - - Constructs a new vector of domain X.domain() and fills it with evaluation of X. -========================================== =========================================================================================== - -* Examples :: - - vector A(10,2); - A()=0; // init A to 0 - -.. warning:: - The constructor from the dimensions does **NOT** initialize the vector to 0 - (because it may not be optimal). - If needed, do it explicitely by `a()=0`; - -Constructors of vector_views ----------------------------------------------- - -Automatic construction -^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -vector_view are normally automatically constructed by making (partial) views, ref:`Slicing`, e.g. :: - - tqa::vector A(2,2); - A(range(),2) ; // --> this makes a view... - A() ; // --> this makes a view over the full array - - -Explicit construction -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -To explicitly make a view of an vector, use make_view :: - - vector A(2,2); - make_view(A) //-> a view... - make_view(A) = 13 ; // to assign e.g. - - - -====================================================================== =========================================================================================================== -Constructors of vector_view Comments -====================================================================== =========================================================================================================== -vector_view(const array_view &) - Copy construction (shallow copy) -vector_view(const T & X) - `[Advanced]` T is any type such that X.indexmap() and X.storage() can be used to construct a view. -====================================================================== =========================================================================================================== - - diff --git a/doc/reference/c++/arrays_old/view_or_not_view.rst b/doc/reference/c++/arrays_old/view_or_not_view.rst deleted file mode 100644 index 9239a8d3..00000000 --- a/doc/reference/c++/arrays_old/view_or_not_view.rst +++ /dev/null @@ -1,80 +0,0 @@ -.. highlight:: c - -A fundamental distinction : value classes vs view classes -================================================================= - -A fundamental distinction in TRIQS is the difference between the `value class` -and the corresponding `view classes`. - -* The **value classes** : e.g. array, matrix, vector, ... - - * They own their data. - * They have a `value semantic`. - * When copied, assigned, they make deep copy. - * They are default constructible. - * They resize themselves under assignment if necessary (like e.g. std::vector). - * They can be put in standard STL containers. - - -* The **view classes** : e.g. array_view, matrix_view, vector_view, ... - - * They are just partial views of the corresponding value class. - * They do not own their data. - * They have a `view semantic`, with reference counting to the data. - * They never make deep copy, either in copying, constructing, transforming to python, - * They can not be resized. - * They are useful to work on arrays, matrices, or on some part thereof, without making copies - (e.g. on a slice of an array, a column of a matrix). - - -Each class is therefore provided in two flavors, and it is crucial for the user -to learn and understand their difference in order to use them appropriately. - -This distinction extends to more complex objects in TRIQS, such as Green functions. - - -Behaviour comparison table between value and view classes ------------------------------------------------------------- - -A detailed comparison of the behaviour of the two type of classes is provided in the following table. - - - -=================== ======================================================================= ====================================================================================== -Topics value class (e.g. array, matrix, vector0 view (e.g. array_view, matrix_view, vector_view) -=================== ======================================================================= ====================================================================================== -Contructors - Constructors create a new fresh memory block. - Constructors only make another view of the data. -Default Contructor - Default constructor creates an empty array - No default constructor (what would it view ?). -Copy contructors - Create a new fresh memory block. - Make another view of the data (shallow copy). - - Make a *true* copy of the data. - *Never* copy data, so they are quite quick. -Assignment (=) - The assignment operator creates a new datablock if size mismatchs. - The assignment operator just copy data into the view. - - Hence, assignment never fail for size reason Behaviour is undefined if the size of the view is too small. - (unless there is a memory allocation exception of course) (it throws if the array has option `Tag::BoundCheck` or if - `TRIQS_ARRAYS_ENFORCE_BOUNDCHECK` is defined, Cf below). -Resizing - Can be resized, invalidating all references/pointers to the data. - Can be not be resized. -Invalidation - References/pointers to the data may become invalid after resize, - References/pointers to the data are still valid after assignment. - or assignment. -STL - Can be used with STL container and algorithms - Can not be put in STL containers, but have STL compliant iterators. -Data Contiguity - The data are contiguous in memory. - The data may not be contiguous in memory (e.g. a slice). -=================== ======================================================================= ====================================================================================== - -An example --------------- - -To be written.... - -A simple example :: - - array f(size_t n) { - array res(n,n); - res() = 0; - return res; - } - - array_view f(size_t n) { - array res(n,n); - res() = 0; - return res; - } - - diff --git a/doc/reference/c++/gf/gf_and_view.rst b/doc/reference/c++/gf/gf_and_view.rst index dcccccb4..a01aaf6a 100644 --- a/doc/reference/c++/gf/gf_and_view.rst +++ b/doc/reference/c++/gf/gf_and_view.rst @@ -7,52 +7,77 @@ gf and gf_view gf is the Green function container, with its view gf_view, defined as :: - template class gf; // the regular type - template class gf_view; // the view type + template class gf; // regular type + template class gf_view; // view type In this section, we describe all common properties of gf. In practice, one works with one of the specialization described below (:ref:`gf_special`), with their specific aspects (domains, interpolation methods, ....). - + Template parameters ---------------------------- -* Variable determines the domain of the Green function. It can be : +* **Variable** determines the domain of the Green function +* **Target** determines the value of the Green function +* **Opt** is currently not used [default to void]. It may be used to differentiate different Green functions + with the same Variable, Target but different behaviour (e.g. interpolation, ...). - ============================ ==================================================================== - Variable tag Meaning - ============================ ==================================================================== - :ref:`imfreq` Imaginary Matsubara frequency - imtime Imaginary Matsubara time - refreq Real frequency - retime Real time - legendre Legendre polynomial representation - block_index Block Green function - cartesian_product Cartesian product of gf ... functions. - ============================ ==================================================================== +They can have the following values : -* Target determines the value of the Green function ++--------------------------+--------------------------------------------+ +| Variable | Meaning | ++==========================+============================================+ +| gf_imfreq | Imaginary Matsubara frequency | ++--------------------------+--------------------------------------------+ +| imtime | Imaginary Matsubara time | ++--------------------------+--------------------------------------------+ +| refreq | Real frequency | ++--------------------------+--------------------------------------------+ +| retime | Real time | ++--------------------------+--------------------------------------------+ +| legendre | Legendre polynomial representation | ++--------------------------+--------------------------------------------+ +| block_index | Block Green function | ++--------------------------+--------------------------------------------+ +| cartesian_product | Cartesian product of gf ... functions. | ++--------------------------+--------------------------------------------+ - ============================ ==================================================================== - Target tag Meaning - ============================ ==================================================================== - matrix_valued [default] The function is matrix valued. - scalar_valued The function is scalar valued (double, complex...). - ============================ ==================================================================== +and -* Opt is currently not used [default to void]. ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ .. _gf_special: Specializations ------------------- ++-----------------+------------------+--------------------+----------------------------+ +| Variable\Target | matrix_valued | scalar_valued | G (another Green function) | ++=================+==================+====================+============================+ +| imfreq | :doc:`gf_imfreq` | :doc:`gf_imfreq_s` | | ++-----------------+------------------+--------------------+----------------------------+ +| imtime | :doc:`gf_imtime` | :doc:`gf_imtime_s` | | ++-----------------+------------------+--------------------+----------------------------+ +| refreq | :doc:`gf_refreq` | :doc:`gf_refreq_s` | | ++-----------------+------------------+--------------------+----------------------------+ +| retime | :doc:`gf_retime` | :doc:`gf_retime_s` | | ++-----------------+------------------+--------------------+----------------------------+ +| block_index | | | :doc:`block_gf` | ++-----------------+------------------+--------------------+----------------------------+ .. toctree:: + :hidden: :maxdepth: 1 gf_imfreq + gf_imfreq_s gf_refreq gf_imtime gf_retime diff --git a/doc/reference/c++/gf/gf_imfreq.rst b/doc/reference/c++/gf/gf_imfreq.rst index 88f2a398..b22ff37a 100644 --- a/doc/reference/c++/gf/gf_imfreq.rst +++ b/doc/reference/c++/gf/gf_imfreq.rst @@ -2,7 +2,7 @@ .. _gf_imfreq: -gf & gf +gf ========================================================== This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara frequencies. @@ -11,6 +11,7 @@ Domain & mesh ---------------- + Singularity ------------- @@ -22,7 +23,7 @@ Factories The factories are :: - make_gf(mesh m, matrix_shape_t shape, local::tail_view t = local::tail(shape) ) + make_gf(gf_mesh m, matrix_shape_t shape, local::tail_view t = local::tail(shape) ) make_gf(double beta, statistic_enum S, matrix_shape_t shape, size_t Nmax = 1025, local::tail_view t = local::tail(shape) ) diff --git a/doc/reference/c++/gf/gf_imfreq_s.rst b/doc/reference/c++/gf/gf_imfreq_s.rst new file mode 100644 index 00000000..0710c380 --- /dev/null +++ b/doc/reference/c++/gf/gf_imfreq_s.rst @@ -0,0 +1,79 @@ +.. highlight:: c + +.. _gf_imfreq_s: + +gf +========================================================== + +This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara frequencies. + +Domain & mesh +---------------- + + +Singularity +------------- + +:ref:`gf_tail`. + +Factories +------------- + + +The factories are :: + + make_gf(mesh m, local::tail_view t = local::tail({1,1}) ) + make_gf(double beta, statistic_enum S, size_t Nmax = 1025, local::tail_view t = local::tail({1,1}) ) + + +Interpolation method +--------------------- + +None + +Data storage +--------------- + +* `data_t` : 1d array (C ordered) of complex. + +* g.data()(i) is the value of g for the i-th point of the mesh. + +HDF5 storage convention +--------------------------- + +h5 tag : `ImFreq_s` + + +Examples +--------- + + +.. compileblock:: + + #include + using namespace triqs::gfs; + int main() { + double beta=1; // inverse temperature + size_t n_freq=5; // we will have 5 points including iw=0 and iw=beta + + auto GF = make_gf(beta, Fermion, n_freq); + }; + + +An alternative declaration with an explicit construction of the underlying mesh: + +.. compileblock:: + + + #include + using namespace triqs::gfs; + int main(){ + double beta=10; + int Nfreq =100; + + auto GF = make_gf(gf_mesh{beta,Fermion,Nfreq}); + // auto GF2 = make_gf({beta,Fermion,Nfreq}); + } + + + diff --git a/doc/reference/c++/gf/plan.rst b/doc/reference/c++/gf/plan.txt similarity index 100% rename from doc/reference/c++/gf/plan.rst rename to doc/reference/c++/gf/plan.txt diff --git a/doc/reference/c++/lectures/basic/C++11.rst b/doc/reference/c++/lectures/basic/C++11.rst deleted file mode 100644 index ee680d75..00000000 --- a/doc/reference/c++/lectures/basic/C++11.rst +++ /dev/null @@ -1,326 +0,0 @@ -.. highlight:: c - -The new C++11 standard -=============================== - -C++11 is the new standard of the C++, voted in 2011 by the International C++ Standard Committee. -It brings several improvements to the C++. -The goal of this page is to highlight a few constructs, which may be useful for us in scientific computing -in order to make C++ more accessible and ease library writing. - -A simple presentation of C++11 is given in `wikipedia `_. - -Implementation ------------------ - -The first thing is to know if/whether this standard is usable, i.e. implemented in compilers. -It turns out that recent versions of compilers (gcc, clang, icc) have already most of the features, -and the implementation is very rapidely progressing. - -A simple `page `_ present what is implemented in various compilers. -A detailed table of the implemented features with versions exists for `gcc `_ -and for `clang `_. - -.. note:: - To use C++11, pass the option -std=c++0x (all compilers), or -std=c++11 (gcc >= 4.7, clang >= 3.0). - -Rational for using C++11 ----------------------------- - -C++11 add some libraries (standardizing several boost libraries, e.g. shared_ptr), -but also and most importantly several new features in the language itself. -Some are basically *syntactic sugar*, some have a deep effect on the programming style. - -Some of these features allow a much cleaner programming with : - -* functionnal programming : introduction of lambda, in particular, also auto, decltype. -* metaprogramming (variadic template, auto, decltype, etc). - -C++11 allows : - -* Write simpler code for the User of libraries -* A lot of simplification in the writing of libraries themselves. - -In my opinion the gain in clarity that one can achieve with C++11 vastly outweigh the -cost in portability (to old versions of compilers). - -There are already many articles on C++11, and soon books, so I will only focuss on -a few examples and basic explanations. - -**Except when specified, the examples below compiles with gcc 4.6 (Ubuntu 12.04 LTS default) and clang 3.1 (current stable)** - -Rvalue references ------------------------ - -Probably one of the most useful feature is the *move semantics*, or rvalue references. - -The problem -^^^^^^^^^^^^^ - - *How to return efficiently a big object built by a function ?* - -Example : I have a function that creates a new std::vector (or any array, Green functions). -How do I return it ? - -A simple solution:: - - std::vector f ( my_parameters) { - std::vector res( vec_size); - // fill res - res[0]=1; //... - return res; - } - /// usage - std::vector v = f(1,2); // constructed a new v from the result of f call. - v = f(3,2); // reassign an already constructed v to the result of v - -Normally, this solution has a major defect : *res* is **copied** when returned from the function, -hence there is a (potentially huge) performance penalty in doing such simple code. -Most current compilers will already elide (remove) the copy in the first case, but it is not guaranteed. -And in the second case, there will be a copy. - -This situation has lead to several workarounds idioms, obscuring the coding style in C++, like :: - - // solution 1 : use pointers - std::vector * f (my_parameters) { // <---- return a pointer - std::vector res = new std::vector(vec_size); // <--- dynamical allocation - // fill res - (*res)[0] = 1; /// <---- not so nice to write - return res; - } - /// usage - std::vector * v = f(1,2); // <--- put this in shared_ptr ?? - *v = f(3,2); // <--- Oops : memory leak ! - -This has a big pb too : who is going to take care of deallocating the pointer (see shared_ptr section) ?? - -Another solution is :: - - // solution 2 : pass by reference - void f( std::vector & res, my_parameters) { // <---- pass a vector and change it - /// fill res.... - } - // usage - std::vector v; - f(v,1,2); // <------ not nice, not natural - - -The C++11 solution -^^^^^^^^^^^^^^^^^^^^^^ - -The **efficient** C++11 code is now **the simplest one** in this situation :: - - std::vector f (my_parameters) { - std::vector res(vec_size); - res[0]=1; // fill res ... - return res; - } - - /// usage - std::vector v = f(1,2); // no copy (copy elision) ! - v = f(2,3); // no copy (move operator =) !! - -In that case, there will be no useless copy. The data of the temporary will be "stolen" by v. -Since the temporary is ... temporary, this is a good strategy. Note that : - -* This is transparent in this code. The user of the vector object, - simply write the **simplest, natural code**. (it is a good idea still to understand the basic idea, see below). - -* The charge of implementing the move semantics rely on the **library writer, not the user of the class**. - -* Still, implementing this is not hard (see example below). - -* A detailed series of (readable) articles on copy elision, and move semantics by D. Abrahams - can be found `here `_. - -An example of a move constructor -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Let us now make a minimal vector class to illustrate the idea and how to implement move semantics. - -In a situation like above, the vector (or anything that provides a **move constructor**) will -alleviate the copy by *stealing the data from the returned temporary*. -The temporary will be discarded, so there is no pb in stealing its data. - -Typically a class with move semantics will provide : - -* a copy constructor, an operator =, as in C++03. -* a move constructor, a move operator= - -.. compileblock:: - - #include - #include - class my_vector { - double * data; // a piece of memory - size_t size_; // size of the piece - public : - my_vector(int s=0): data(s>0 ? new double[s] : NULL ), size_(s){} // basic constructor. Just allocate memory - ~my_vector() { if (data) delete[] data;} // free the pointer - - my_vector(my_vector const & x) { // copy contructor - std::cout << "copy construction"<< std::endl; - size_ = x.size(); data = new double[size_]; - memcpy(this->data, x.data, x.size()*sizeof(double)); // copy data in memory (low level C) - } - - my_vector & operator=( my_vector const & x) { // normal = operator - std::cout <<" normal ="<< std::endl; - size_ = x.size(); if (data) delete[] data; data = new double[size_]; - memcpy(this->data, x.data, x.size()*sizeof(double)); // copy data in memory (low level C) - return * this; - } - - my_vector(my_vector && V) { // C++11 move constructor - std::cout <<" move constructor"<< std::endl; - data = V.data; V.data = NULL; // swap the data with V and let V empty ! - size_ = V.size_; V.size_ = 0; - } - - my_vector & operator=( my_vector && V) { // C++11 move = operator - std::cout <<" Move ="<< std::endl; - std::swap(V.data, this->data); // swap data - std::swap(V.size_, this->size_); - return * this; - } - // ------ standard interface to vector (minimal) -------------------- - double const & operator[](size_t i) const { return (data)[i];} - double & operator[](size_t i) { return (data)[i];} - size_t size() const { return size_;} - friend std::ostream & operator <<( std::ostream & out, my_vector const & v) { - for (size_t i =0; i A(2,3), B(2,3); - auto expr = A + 2*B; - - -for range loop ------------------ - -For any object that have iterators (in fact it is more general ....), a simpler syntax to iterate on them is. -Example :: - - std::vector V = {1,2,3}; - for (auto & x : V) { - // do something with x ....(it is a reference) - } - - -Initializer list ------------------- - -Example :: - - std::vector V = {1,2,3}; - triqs::arrays::matrix M = {{1,2},{3,4}}; // not implemented yet.... - -Lambda ------------- - -As in python, C++11 now supports simple lambda functions :: - - auto f = [](int) { return x+1;} - int u = f(3); // use as a regular function... - -* It still allow inlining. - -* This is particularly useful with the STL algorithms, or similar iteration algorithms : - -.. - .. literalinclude:: count_if.cpp - -.. compileblock:: - - #include - #include - #include - - int main() - { - std::vector v = { -1, -2, -3, 4, -4, 3, -7, 8, 9, 10 }; - int num_items1 = std::count_if(v.begin(), v.end(), [](double i) {return i>0;}); - std::cout << "number of >0 : " << num_items1 << '\n'; - - int bound = 3; - num_items1 = std::count_if(v.begin(), v.end(), [&bound](double i) {return i>bound;}); - std::cout << "number of >bound : " << num_items1 << '\n'; - } - - -Variadic template ---------------------- - -Improvement of the standard library ---------------------------------------- - -Most notably : - -* std::shared_ptr, and std::unique_ptr (see Smart Pointers). -* std::unordered_map - -In most cases, this is a standardisation of some boost libraries. - diff --git a/doc/reference/c++/lectures/basic/compile_time_decision.rst b/doc/reference/c++/lectures/basic/compile_time_decision.rst deleted file mode 100644 index fb1a1114..00000000 --- a/doc/reference/c++/lectures/basic/compile_time_decision.rst +++ /dev/null @@ -1,13 +0,0 @@ -Various ways to take a decision at compile time -############################################### - -* Overloading : often forgotten. Use of tags - -* std::conditionnal or mpl::if\_ - -* std::enable_if - In combination with a template <> + 1 default ... - -* template technique : specialisation, traits ... - - diff --git a/doc/reference/c++/lectures/basic/contents.rst b/doc/reference/c++/lectures/basic/contents.rst deleted file mode 100644 index 3fa06ed3..00000000 --- a/doc/reference/c++/lectures/basic/contents.rst +++ /dev/null @@ -1,16 +0,0 @@ -General C++ techniques and idioms -************************************************************************************ - -.. highlight:: c - -.. toctree:: - :maxdepth: 1 - :numbered: - - C++11 - concepts - smart_ptr - functions - compile_time_decision - - diff --git a/doc/reference/c++/lectures/basic/functions.rst b/doc/reference/c++/lectures/basic/functions.rst deleted file mode 100644 index 13c73c1f..00000000 --- a/doc/reference/c++/lectures/basic/functions.rst +++ /dev/null @@ -1,14 +0,0 @@ -Functions becomes first class citizens ... -################################################# - -* functional programming : another paradigm ... - -* std::function : store a function in a an object -* currifying: std::bind, std::placeholders. -* Do not USE function pointers.... -* Function objects : wrap functions into an objects.... -* lambda : example of use, STL.... - - - - diff --git a/doc/reference/c++/lectures/basic/smart_ptr.rst b/doc/reference/c++/lectures/basic/smart_ptr.rst deleted file mode 100644 index 067969dc..00000000 --- a/doc/reference/c++/lectures/basic/smart_ptr.rst +++ /dev/null @@ -1,182 +0,0 @@ -.. highlight:: c - -Pointers : from the evil to the smart... -###################################################### - -Raw pointers -=============== - -Raw pointers are ordinary C/C++ pointers. - -Pointers allow to manipulate objects by reference (or indirectly). - -RefTo : stack/heap. Pointers in C++... - -They are used e.g. to allocate a new object on the heap :: - - { - A * p = new A( constructor_arguments); // allocated and constructed on the heap. - A a( constructor_arguments); // constructed on the stack - - delete p; // release memory - } - -Raw pointers are the source of many bugs in C++. Two main dangers are : - -* memory leak : when one forget to release the memory allocated to the pointer:: - - { - A * p = new A( constructor_arguments); // allocated and constructed on the heap. - // - // forgot to delete. but p is gone.... - } - -* a more subtle memory leak :: - - { - A * p = new A( constructor_arguments); // allocated and constructed on the heap. - // working with p... - delete p; - } - - What happens if there is an exception during the work ?? - The pointer is not deallocated (but all objects on stack would be destroyed). - -* dangling pointers : reusing a pointer after the object has been deleted :: - - { - A * p = new A( constructor_arguments); // allocated and constructed on the heap. - //... - A *q =p; a copy of the pointer is made somewhere - //... - delete p; - // ... - q->method(..) ; // DEAD : segfault. - } - -To conclude : - -* pointers are a low level object - -* pointers induce many dangers... - -The solution -==================== - -Smart pointer are a way to manipulate objects by reference with strong guarantee -that such common problems will **never** happen. - -The bottom line can be summarized as : - - *Smart pointers are pointers with strong guarantee to avoid dangling references and memory leaks.* - - -In order to avoid pointers problems, the style recommendation are : - -#. Try to avoid pointers ! Use object by value, and objects like std::vector - which hide the dynamic allocation of memory from the user in a safe manner. - -#. If you need a pointer to an object A which is pointed to by exactly one pointer - (a local object), then use std::unique_ptr (or boost:scoped_ptr if you insist on non C++11). - -#. If you need an object A shared by multiple objects, i.e. pointed to by several pointers - use a std::shared_ptr. - -#. Only when other options are not good, and in practice only at low level libraries, - use raw pointers to allocate objects, with new and delete. - -std::unique_ptr ---------------------- - -Example :: - - { - std::unique_ptr p( new A( constructor_arguments) ); // allocated and constructed on the heap and stored. - // - // working with p. p has the semantics of a pointer : - // p-> method(), p->data, etc ... - // NO DELETE : when the std::unique_ptr is getting out of scope and is destroyed the - // underlying raw pointer are freed. - } - -Note that : - -* a unique_ptr **can not be copied**. If you put it as a member of a class, this class can therefore not be copied. - -* a unique_ptr **can be moved** (see lecture on move semantics). - -* This writing is exception-safe. If an exception occurs in the treatment, then C++ guarantees - the all the local objects are destroyed, hence p and the underlying raw pointer is deleted. - -* unique_ptr has little if no overhead in critical code. - - -std::shared_ptr ---------------------- - -(or boost::shared_ptr < C++11). - -In the case where we want to **share** an object, i.e. that multiple pointers points to an object, -we have to use std::shared_ptr. - -A shared_ptr is no more than : - -* a pointer to the object - -* a reference counter, counting how many shared_ptr point to the objects. - - -When the shared_ptr are copied, assigned, the reference counter is updated accordingly. - -* When a shared_ptr is created, copied, ref_counter is increased by 1. -* When a shared_ptr is deleted, ref_counter is decreased by 1. - -When the last shared_ptr to an object is destroyed, the reference counter gets -to 0 and the pointed objects is deleted. - -Example :: - - class A; - - boost::shared_ptr p (construction_parameters); - - // work with p like a pointer - // p has pointer semantic. *p, p->method as for raw pointer - - // NO delete - - -Example :: - - class A; - - boost::shared_ptr p1 (construction_parameters); - boost::shared_ptr p2 (construction_parameters); - - // p1 ----> ( a1, ref_counter = 1 ) - // p2 ----> ( a2, ref_counter = 1 ) - - { // group start - boost::shared_ptr p3 = p2; - // p1 ----> ( a1, ref_counter = 1 ) - // p2 ----> ( a2, ref_counter = 2 ) - // p3 ----> ( a2, ref_counter = 2 ) - - } // group ends ... p3 is deleted - - // p1 ----> ( a1, ref_counter = 1 ) - // p2 ----> ( a2, ref_counter = 1 ) - - p1 = p2; - // p1 ----> ( a1, ref_counter = 2 ) - // p2 ----> ( a1, ref_counter = 2 ) - // a2 has been destroyed. - - - - - - - - - diff --git a/doc/reference/c++/lectures/contents.rst b/doc/reference/c++/lectures/contents.rst deleted file mode 100644 index 50a6520c..00000000 --- a/doc/reference/c++/lectures/contents.rst +++ /dev/null @@ -1,30 +0,0 @@ -Notes for TRIQS C++ contributor -********************************* - - *Our goal : write clear, readeable, high-level and yet efficient programs.* - -The goal of these notes is to be : - -* a starting point for learning modern C++ for scientific computation : - Traditionnally, it is often said that genericity, abstraction (i.e. high-level) programming - is not compatible with the efficency that we need in scientific computations. - The goal of these lectures is to show that we can achieve these goals with today C++ and libraries. - -* a style guide for the TRIQS project (C++ part). - -Most of the content is more general (and independant of the TRIQS project). -They are a few notes on modern C++ for scientific computations. - -The intended audience is physicists with a basic knowledge of C++ (??). -These lectures do not focuss on basic syntax, but on useful **idioms**. - - -.. highlight:: c - -.. toctree:: - :maxdepth: 1 - :numbered: - - basic/contents - patterns/contents - concepts/contents diff --git a/doc/reference/c++/lectures/goal.rst b/doc/reference/c++/lectures/goal.rst deleted file mode 100644 index 9c1b9197..00000000 --- a/doc/reference/c++/lectures/goal.rst +++ /dev/null @@ -1,25 +0,0 @@ -Introduction -============ - -*Our goal : write clear, readeable, high-level and yet efficient programs.* - -This is the motto of the lectures. -Traditionnally, it is often said that genericity, abstraction (i.e. high-level) programming -is not compatible with the efficency that we need in scientific computations. -The goal of these lectures is to show that we can achieve these goals with today C++ and libraries. - -Why is it important? --------------------- - -* Our codes are changing all the times. We want to have code reuse. - -* more to come here - - -Goal of the lectures --------------------- - -* Present some patterns and techniques which are useful to write high level and efficient codes. -* Fix the "pattern" design of TRIQS. One problem : one standardised solution. - - diff --git a/doc/reference/c++/lectures/intro.rst b/doc/reference/c++/lectures/intro.rst deleted file mode 100644 index c2306552..00000000 --- a/doc/reference/c++/lectures/intro.rst +++ /dev/null @@ -1,31 +0,0 @@ - -User Manual : the complete (?) reference -******************************************* - -.. warning:: - This library is beta. - - This manual is not finished: work in progress. - - -.. highlight:: c - -.. toctree:: - :maxdepth: 2 - :numbered: - - view_or_not_view - basic_classes - slicing - shape - debug - assignment - algebras - foreach - functional - STL - H5 - Interop_Python - IO - blas_lapack - FAQ diff --git a/doc/reference/c++/lectures/patterns/assign_delegation.rst b/doc/reference/c++/lectures/patterns/assign_delegation.rst deleted file mode 100644 index b243a4e2..00000000 --- a/doc/reference/c++/lectures/patterns/assign_delegation.rst +++ /dev/null @@ -1,133 +0,0 @@ -.. highlight:: c - - -Assignment Delegation -=============================== - -Problem ------------- - -Given a class A which defines an `operator =`, we would like to -be able to implement various optimisation depending on the RHS (right hand side). - -Normally, we can simply define various operator in class A :: - - class A { - // ... - template A & operator = (RHS rhs); - } - -maybe with some enabler to choose for different types. - -Normally RHS are expected to model some precise concept (e.g. array, green function, etc...). -In other words, RHS are looking "like an A". -It is then easy to implement a generic version of operator = that uses just the concept. - -We would like to be able to implement various optimisations depending on the RHS (right hand side). -This optimisation are decided are compile time and maybe also at run time. - -Of course, we could split the definition of operator = into cases with enable_if. -But we want a non-intrusive solution, i.e. be able to add an optimisation -for a new type of RHS, without modifying the class A. - -Real life examples -^^^^^^^^^^^^^^^^^^^^^^^^ - -* array, matrix, vector, with RHS a generic proto expression. For example in triqs::arrays :: - - M = inverse(M); // (1) - M = 2*inverse(M); // (2) - - inverse returns an inverse_lazy object. We would like that the inversion is done in place, *without - temporary* in this (1), but obviously not for a general expression like (2). - -* Green function, for example in the case :: - - gf = fourier (gt) ; - - The Fourier transformation could be done directly in the gf data, without a temporary, - if this becomes a bottleneck the code (?). - - The point here is also that we can first implement a generic algorithm. If we need, - we know that we will be able to optimise this temporary *without any change in the user code*. - -Solution ------------------- - -A simple solution is to delegate the assignement and simply use overload e.g. :: - - template - void assign_delegation (LHS & lhs, RHS const & rhs) { - // implement the generic algorithm - } - - class A { - // ... - template A & operator = (RHS const & rhs) { assign_delegation(*this, rhs);} - } - -And then later when defining MyRHS_with_optimisation :: - - class MyRHS_with_optimisation { - // .... implement the object ... - - template - friend void assign_delegation (LHS & lhs, MyRHS_with_optimisation const & rhs) { - // implement the special algorithm - } - }; - -**NB** : Adding the function as a friend in the class has 2 advantages over defining it outside: - - * It is friend, hence it can access the data of rhs. - * It will be injected in the namespace where the class MyRHS_with_optimisation, and therefore it will be found by ADL. - -Variant -^^^^^^^^^^ - -With the same technique, one can also delegate the compound assignments : +=, -=, *=, /=. - -Let us label these operators by a char :A, S, M, D. We can then differentiate between the various -case by adding an additionnal dummy parameter of mpl::char_ type:: - - template - void compound_assign_delegation (LHS & lhs, RHS const & rhs, mpl::char_) { - // implement the generic algorithm - } - - class A { - - template A & operator += (RHS const & rhs) { - compound_assign_delegation(*this, rhs, mpl::char_<'A'>() ); - } - - // and so on... - } - -And then later when defining MyRHS_with_optimisation :: - - class MyRHS_with_optimisation { - // .... implement the object ... - - template - friend void compound_assign_delegation (LHS & lhs, MyRHS_with_optimisation const & rhs,mpl::char_<'E'>) { - // implement the special algorithm - } - }; - -**NB** : - -* One can then overload some or all of the operations for a given class. - -* mpl::char_<'A'> is a different **type** than mpl::char_<'S'>, so the standard overload mechanism of C++ will - choose the right function at compile time (the `mpl::char_` object is tiny, and will be optimized away). - It better than passing a simple char and make a switch at run time (because there is no compile time choice in this case). - -* One could a priori also add a char OP in the template parameter list, but then in the MyRHS_with_optimisation class, we would need - to specialize a function, which results in less readeable code, IMHO. - -Examples ------------ - -Used in triqs::arrays. - diff --git a/doc/reference/c++/lectures/patterns/contents.rst b/doc/reference/c++/lectures/patterns/contents.rst deleted file mode 100644 index 817cb031..00000000 --- a/doc/reference/c++/lectures/patterns/contents.rst +++ /dev/null @@ -1,24 +0,0 @@ -.. highlight:: c - -Some useful idioms and patterns -####################################### - -The goal of this list of patterns is to : - -* illustrate some useful (?) patterns in C++. -* document TRIQS code and standardize the way things are implemented in TRIQS - - -.. toctree:: - :maxdepth: 1 - - h5_read_write - lazy_objects - type_erasure - python - assign_delegation - proto - preproc - iterator - - diff --git a/doc/reference/c++/lectures/patterns/h5_read_write.rst b/doc/reference/c++/lectures/patterns/h5_read_write.rst deleted file mode 100644 index f2695e3c..00000000 --- a/doc/reference/c++/lectures/patterns/h5_read_write.rst +++ /dev/null @@ -1,148 +0,0 @@ -.. highlight:: c - -HDF5 write/read convention -================================ - - -All TRIQS objects that can be h5-serialized (i.e. put into an hdf5 file and reconstructed from it) -follow the following convention. - -Convention --------------- - -For each h5-serializable TRIQS object, one defines **free functions** *h5_read*, *h5_write* with the synopsis :: - - // Write the object into the group or file F, either directly (for array, basic type) - // or by opening a subgroup and writing its components into this subgroup. - h5_read ( h5::group_or_file F, std::string const & Name, MyNiceObject & obj); - - // The inverse operation. It fills the object from the file. The object must exists before. - h5_write ( h5::group_or_file F, std::string const & Name, MyNiceObject const & obj); | - -* h5::group_or_file is a little opaque object that can be transparently constructed - from the two types H5::Group and H5::H5File (respectively the Group and the File in the HDF5 library). - -* Recommended practice : **define these functions as friend functions of the class**. - - *Rational* : they will access the data and by resolved by ADL. - -Example : ----------- - -Suppose I have an object with one array A and a string: - -.. compileblock:: - - #include - #include - namespace tqa= triqs::arrays; - - class MyNiceObject { - tqa::array A; - double x; - std::string s; - public : - MyNiceObject(): A(2,3), x(3), s("wonderful") { A() = 3;} - - friend void h5_write (triqs::h5::group fg, std::string subgroup_name, MyNiceObject const & obj) { - auto g = fg.create_group(subgroup_name); - h5_write(g,"A",obj.A); - h5_write(g,"x",obj.x); - h5_write(g,"s",obj.s); - } - - friend void h5_read (triqs::h5::group fg, std::string subgroup_name, MyNiceObject & obj){ - auto g = fg.open_group(subgroup_name); - h5_read(g,"A",obj.A); - h5_read(g,"x",obj.x); - h5_read(g,"s",obj.s); - } - }; - - int main() { - H5::H5File file("ess.h5", H5F_ACC_TRUNC ); - MyNiceObject a; - h5_write(file, "a1", a); - h5_write(file, "a2", a); - } - -The generated file looks like:: - - > h5ls -r ess.h5 - / Group - /a1 Group - /a1/A Dataset {2, 3} - /a1/s Dataset {SCALAR} - /a1/x Dataset {SCALAR} - /a2 Group - /a2/A Dataset {2, 3} - /a2/s Dataset {SCALAR} - /a2/x Dataset {SCALAR} - -and the full dump as :: - - HDF5 "ess.h5" { - GROUP "/" { - GROUP "a1" { - DATASET "A" { - DATATYPE H5T_IEEE_F64LE - DATASPACE SIMPLE { ( 2, 3 ) / ( 2, 3 ) } - DATA { - (0,0): 3, 3, 3, - (1,0): 3, 3, 3 - } - } - DATASET "s" { - DATATYPE H5T_STRING { - STRSIZE 9; - STRPAD H5T_STR_NULLTERM; - CSET H5T_CSET_ASCII; - CTYPE H5T_C_S1; - } - DATASPACE SCALAR - DATA { - (0): "wonderful" - } - } - DATASET "x" { - DATATYPE H5T_IEEE_F64LE - DATASPACE SCALAR - DATA { - (0): 3 - } - } - } - GROUP "a2" { - DATASET "A" { - DATATYPE H5T_IEEE_F64LE - DATASPACE SIMPLE { ( 2, 3 ) / ( 2, 3 ) } - DATA { - (0,0): 3, 3, 3, - (1,0): 3, 3, 3 - } - } - DATASET "s" { - DATATYPE H5T_STRING { - STRSIZE 9; - STRPAD H5T_STR_NULLTERM; - CSET H5T_CSET_ASCII; - CTYPE H5T_C_S1; - } - DATASPACE SCALAR - DATA { - (0): "wonderful" - } - } - DATASET "x" { - DATATYPE H5T_IEEE_F64LE - DATASPACE SCALAR - DATA { - (0): 3 - } - } - } - } - } - - - diff --git a/doc/reference/c++/lectures/patterns/lazy_objects.rst b/doc/reference/c++/lectures/patterns/lazy_objects.rst deleted file mode 100644 index 6a0ae9d7..00000000 --- a/doc/reference/c++/lectures/patterns/lazy_objects.rst +++ /dev/null @@ -1,47 +0,0 @@ -.. highlight:: c - -Lazy Objects -==================== - -PreCompute in Lazy Object ------------------------------ - -* why a lazy object -* example - -* idea : keep a view or a copy and then implement the concept. -* add possibly some optimisation - -* When precomputation are needed. - --> LazyPreCompute pattern with the shared_ptr - - Ex :: - - class MyLazyObject { - - // .. - private: - - struct internal_data { - // contains the data for the precomputation - // e.g. std::vector data; - internal_data(MyLazyObject const & P) { - // do the computation and fill the data - // data = ...; - } - }; - friend struct internal_data; // internal_data can read the parent object P - mutable boost::shared_ptr _id; // holds the data - void activate() const { if (!_id) _id= boost::make_shared(*this);} - - } - -Comments : - -* The data is hold in a specialized struct `internal_data`, which contains the precomputation. -* The data is hold by a shared_ptr, _id : - - * The copy of MyLazyObject is then cheap and safe with shared_ptr. - * If the computation is not necessary, _id is empty, there is no memory usage. - * activate() activates the calculation of the internal data. - * _id is mutable : the precomputation is private, so even a const object should be able to activate safely. diff --git a/doc/reference/c++/lectures/patterns/python.rst b/doc/reference/c++/lectures/patterns/python.rst deleted file mode 100644 index 93b5b447..00000000 --- a/doc/reference/c++/lectures/patterns/python.rst +++ /dev/null @@ -1,57 +0,0 @@ -Python bindings -##################### - -All the idioms that we use from boost.python .... - -Introduction to boost.python -============================== - -* wrapping versus converting. - -Example of wrapping a function, a class, with an array. - -Using a more complex converters. - -Wrapping enable a converter. - -* converting : principle. - Ref to doc. - The simple triqs way to build a converter. - Example of range ! - -* Passing a pyobject. - -* Difference between PyObject* and bpy::object - Reference counting. - -* More sophisticated data structure. - -* THE TRIQS tools : a thin layer (converter) for recursive conversion - (e.g. vector , unordered_map, ....) - - -boost.python : a C++ API for python -====================================== - -Code injection -================= - -The problem. - -The solution : injection and documentation. - -Examples -========== - -Passing a dict for parameters in a routine. - -Passing a list of gf_view into a vector of gf_view... ----> make it a const vector to avoid changing it.... - -Other options -============== - -* Swig -* Cython - - diff --git a/doc/reference/c++/lectures/patterns/type_erasure.rst b/doc/reference/c++/lectures/patterns/type_erasure.rst deleted file mode 100644 index 270c74db..00000000 --- a/doc/reference/c++/lectures/patterns/type_erasure.rst +++ /dev/null @@ -1,55 +0,0 @@ -Type erasure : from static to dynamical polymorphism -############################################################ - -* What is the issue ? -* A basic example : std::function -* Example : -* A real life example : QMC class. -* The classical C++ solution with virtual functions - - - -Polymorphism -================ - -Static polymorphism (template, generic programming) ---------------------------------------------------------- - -Dynamical polymorphism. - -The crossing : why ? How ? -Concepts. When is the decision taken : at compile time or a run time ? - -A simple example : std::function -Can handle any function or object which ahs the concept double(int,int) e.g. -It **erases** the type of the object. - -The pb : find an envelop for an object with : - * the concept - * one type, that can handle any type. ----> Add Abrahams definitions ( contruct + =). - -An example : - -Pb is twofold : -* keep the undelying alive -* call its methods wihout remembering its type. - ---> show a solution -Use std::function + std::bind variant without lambdas. -Use shared_ptr -* shared_ptr the universal handle - ---> Hum : pb : restore the value semantics ! - -* The traditionnal C++ way using 3 classes. ugly... - -The envelop as a concept checker ??? ---> It does not compile it T does not model the concept. - -Python : a general type eraser : PyObject * or boost::python::object. - -2 frontiers : template, C++ polymorphism, Python : in easy and in speed. - - - diff --git a/doc/reference/c++/lectures/patterns/view_value.rst b/doc/reference/c++/lectures/patterns/view_value.rst deleted file mode 100644 index 8a4fcb02..00000000 --- a/doc/reference/c++/lectures/patterns/view_value.rst +++ /dev/null @@ -1,18 +0,0 @@ -ViewValue pattern in TRIQS -################################## - -* for array, matrix, vector, tail, gf of any kind. -* always of the same form. - -* A value class with value semantics -* A view class to take partial views & slices, with view semantics. - -How to implement ? - -* X_impl - --> X and X_view by derivation : reimplement only construction (PATTERN) and =. - -* Rational : - do NOT derive X from X_view ! - - diff --git a/doc/reference/c++/lectures/patterns/writing_an_iterator.rst b/doc/reference/c++/lectures/patterns/writing_an_iterator.rst deleted file mode 100644 index 865089d1..00000000 --- a/doc/reference/c++/lectures/patterns/writing_an_iterator.rst +++ /dev/null @@ -1,42 +0,0 @@ -How to write quickly and efficently an iterator -====================================================== - -* Do not use too much iterators, usually for loops are more efficient. - -* Iterators are mainly useful in combination with STL. - -* The best way to have an STL compliant iterator is to use the boost::iterator_facade documented - `here `_ - -Example (stupid since std::vector already has an iterator) :: - - #include - - template - class my_iterator : public boost::iterator_facade< my_iterator , const std::ptrdiff_t , boost::forward_traversal_tag > { - public: - typedef CuboidMap indexmap_type; - typedef typename indexmap_type::domain_type domain_type; - typedef IterationOrder iteration_order; - typedef typename domain_type::index_value_type indices_type; - typedef const std::ptrdiff_t return_type; - - my_iterator (): Parent(NULL), pos(0),atend(true) {} - my_iterator (const indexmap_type & P, bool atEnd=false): Parent(&P), pos(Parent->start_shift()),atend(atEnd) {} - - indices_type const & indices() const { return indices_tuple; } - operator bool() const { return !atend;} - - private: - friend class boost::iterator_core_access; - void increment(); // below - bool equal(my_iterator const & other) const {return ((other.Parent==Parent)&&(other.atend==atend)&&(other.pos==pos));} - return_type & dereference() const { assert (!atend); return pos; } - - const indexmap_type * Parent; - indices_type indices_tuple; - std::ptrdiff_t pos; - bool atend; - }; - - diff --git a/doc/reference/c++/lectures/smart_pointers.cpp b/doc/reference/c++/lectures/smart_pointers.cpp deleted file mode 100644 index f4b95197..00000000 --- a/doc/reference/c++/lectures/smart_pointers.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include -#include - - -class myObject{ - - int i; - - public: - - myObject(){} - myObject(int i_):i(i_){} - int getValue() const {return i;} -}; - - - -struct A{ - - -}; - -struct B{ - int data; - boost::shared_ptr p; -}; - -struct C{ - boost::shared_ptr p; -}; - - -int main(){ - - /* WHY POINTERS? */ - - - - - myObject A(0); // object is allocated on the stack - - /* stack VS heap */ - // small variables on the stack - // dynamical memory on the heap == "new" - - //when running, the OS has allocated the stack size... - //when allocatinh by new ==> HEAP, one asks the kernel for space in the system ... takes time! - //RULE: in critical loops, do not allocate using new ! - - //stack: only for objects whose size is known at compile time - - //array lib : minivector on stack, to avoid allocation problems (ulimit... - - - std::cout << A.getValue()<getValue() << std::endl; - std::cout << (*p1).getValue() << std::endl; - - - - /* THE ISSUE */ - delete(p1); // deletes the pointer but not the underlying data - - - - /* SHARED POINTER */ - boost::shared_ptr p1_sh(new myObject(1)); - - std::cout << p1_sh->getValue()< p2_sh = p1_sh; - - std::cout << p2_sh->getValue()<