diff --git a/cmake/Doxyfile.in b/cmake/Doxyfile.in index e69e928f..32b9c0f2 100644 --- a/cmake/Doxyfile.in +++ b/cmake/Doxyfile.in @@ -521,7 +521,7 @@ EXCLUDE_PATTERNS = # The symbol name can be a fully qualified name, a word, or if the wildcard * is used, # a substring. Examples: ANamespace, AClass, AClass::ANamespace, ANamespace::*Test -EXCLUDE_SYMBOLS = *::*details::* +EXCLUDE_SYMBOLS = *::*details::*, _aux*, *::gfs_implementation::* # The EXAMPLE_PATH tag can be used to specify one or more files or # directories that contain example code fragments that are included (see diff --git a/doc/index.rst b/doc/index.rst index 444bfb7b..e3819ad9 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -19,9 +19,9 @@ TRIQS (**T**\oolbox for **R**\esearch on **I**\nteracting **Q**\uantum **S**\yst is a scientific project providing a set of C++ and Python libraries to develop new tools for the study of interacting quantum systems. -The goal of this toolkit is to provide condensed matter physicists with +The goal of this toolkit is to provide our team with some high level, efficient and simple to use libraries in C++ and Python, -and to promote the use of modern programming techniques in our field. +and to promote the use of modern programming techniques. TRIQS is free software (GPL). @@ -29,7 +29,7 @@ TRIQS applications ----------------------- Based on the TRIQS toolkit, several :ref:`full-fledged applications ` -are also maintained by the TRIQS collaboration. They allow for example to +are also available. They allow for example to solve a generic quantum impurity model or to run a complete LDA+DMFT calculation. diff --git a/doc/installation/osx_lion.rst b/doc/installation/osx_lion.rst index 04755a23..aeb0a3ab 100644 --- a/doc/installation/osx_lion.rst +++ b/doc/installation/osx_lion.rst @@ -4,28 +4,49 @@ .. _install_on_osx_lion: -Installing required libraries on Mac OS X -========================================= +Installing required libraries on Mac OS X [EXPERIMENTAL] +============================================================== -This is an installation guide for Os X Mountain Lion. -It may work for older -versions of Mac OS X but previous versions of Mac OS X are not supported. +Disclaimer +------------- -NB: The installation of TRIQS under previous versions of OS X requires installing clang (via Xcode). -(On Mountain Lion, clang (llvm) replaces gcc as the default C++ compiler). +We provide here some instructions to install and use TRIQS on OS X. -We strongly recommend the following installation procedure, which provides a clean way to set up all dependencies, so that all -of them are compatible with each other. Only the installation via homebrew is supported for the Mac. +While the installation of TRIQS itself is as straightforward as on Linux systems, +the installation of the standard scientific libraries used by TRIQS +(mpi, hdf5, boost, fftw, ipython, ....) is not as simple as in e.g. Debian/Ubuntu, +where they are packaged with the distribution itself (i.e. "part of the system"). -.. warning:: +This general problem of scientific computing on OS X is clearly +illustrated by the large number of third-party attempts +to provide "easy" installation of scientific libraries : +brew, macports, fink, Enthought Python distribution. + +In our experience, none of this solution is perfect, nor complete or stable : they are still +very far from the quality and stability of a Linux distribution like Debian/Ubuntu. +It *almost* works, but at the end, there are some issues, sometimes severe, +in the way scientific librairies are installed. +(e.g. currently the default version of mpi and hdf5 installed by brew are in conflict : +the simple mpi "Hello World" crashes when linked with hdf5_cpp). + +Moreover, because there is no notion of "distribution" (except in Enthought, which unfortunately is incomplete), +the versions of the libraries are always changing e.g. in brew. +As a result, the installation instructions may work on one day, and suddenly stop to work +the day after. + +While this has a priori **nothing to do with TRIQS** and its applications, it clearly impacts its installation and usage. +We are looking for a more robust solution to this OS X installation mess; **help welcome !**. + +In the following, we describe an installation procedure which worked (at least one day), +on 10.8 and 10.9 (at least on the Mac of one of the developer !). - Because brew evolve with time, there is no notion of distribution on Mac, like e.g. Ubuntu. - So, while the procedure worked at some point, there can be no guarantee that it still does. - Installation of the dependencies -------------------------------- +We describe an installation procedure which is known to have worked at least one day, +on 10.8 and 10.9, (at least on the Mac of one of the developer !). + 1. Install `homebrew `_. Run ``brew doctor`` and resolve potential conflicts before continuing. @@ -33,7 +54,8 @@ Installation of the dependencies 2. Install XCode (directly from the Mac store). In Preferences/Downloads, install "Command Line tools". 3. Install several packages which are needed: :: - + + brew tap homebrew/science brew install cmake brew install gfortran brew install --enable-cxx hdf5 @@ -42,11 +64,9 @@ Installation of the dependencies brew install open-mpi brew install zmq brew install python - brew install doxygen - - #brew formula has been repaired. Temporary using our own - #until this is back in the master. - #When 1.55 is out, the regular brew formula should work again ... + + #brew formula has been repaired, since boost installation of mpi.python is a complete mess + #which needs to be fixed manually (except in Debian/Ubuntu where it is correct). ### brew install boost --without-single --with-mpi --with-c++11 brew install http://ipht.cea.fr/triqs/formulas/boost.rb --without-single --with-mpi --with-c++11 -v @@ -69,22 +89,32 @@ Installation of the dependencies 6. If you wish to compile the documentation locally, install sphinx, its dependencies and mathjax: :: + brew install doxygen pip install sphinx easy_install pyparsing==1.5.7 git clone https://github.com/mathjax/MathJax.git MathJax NB : you need pyparsing <=1.5.7 since apparently v.2.0 works only for python 3. +7. If you wish to build the documentation locally, + configure TRIQS with the option -DPython_use_mpi4py=ON (workaround boost.mpi.python bug). + +8. **Set up** the environment variable, e.g. in your ~/.bash_profile (workaround for issue #43) :: + + export HDF5_DEBUG="all" + + or your code will crash when launched without mpirun + (due to a bug in hdf5 C++/ openmpi, nothing to do with TRIQS, so we can not fix it). Possible issues --------------- -If you encounter the following error: :: +* If you encounter the following error: :: /usr/local/include/ft2build.h:56:38: error: freetype/config/ftheader.h: No such file or directory -in the installation of matplotlib, you need to pass the proper include path. Locate the freetype directory -with the header file and pass the include path through ``CPPFLAGS``: :: + in the installation of matplotlib, you need to pass the proper include path. Locate the freetype directory + with the header file and pass the include path through ``CPPFLAGS``: :: CPPFLAGS=-I/usr/X11/include/freetype2/ pip install git+https://github.com/matplotlib/matplotlib.git#egg=matplotlib-dev diff --git a/doc/reference/c++/arrays/CMakeLists.txt b/doc/reference/c++/arrays/CMakeLists.txt index fc43772c..12643128 100644 --- a/doc/reference/c++/arrays/CMakeLists.txt +++ b/doc/reference/c++/arrays/CMakeLists.txt @@ -1,10 +1,10 @@ # Doxygen sources set_property(GLOBAL APPEND PROPERTY DOXYGEN_SOURCES - ${TRIQS_SOURCE_DIR}/triqs/arrays/h5/simple_read_write.hpp + #${TRIQS_SOURCE_DIR}/triqs/arrays/h5/simple_read_write.hpp ${TRIQS_SOURCE_DIR}/triqs/arrays/h5/array_stack.hpp ${TRIQS_SOURCE_DIR}/triqs/h5/group.hpp - ${TRIQS_SOURCE_DIR}/triqs/arrays/array.hpp - ${TRIQS_SOURCE_DIR}/triqs/arrays/matrix.hpp - ${TRIQS_SOURCE_DIR}/triqs/arrays/vector.hpp + # ${TRIQS_SOURCE_DIR}/triqs/arrays/array.hpp + #${TRIQS_SOURCE_DIR}/triqs/arrays/matrix.hpp + #${TRIQS_SOURCE_DIR}/triqs/arrays/vector.hpp ) diff --git a/doc/reference/c++/gf/CMakeLists.txt b/doc/reference/c++/gf/CMakeLists.txt index 88e36173..97327f47 100644 --- a/doc/reference/c++/gf/CMakeLists.txt +++ b/doc/reference/c++/gf/CMakeLists.txt @@ -1,3 +1,10 @@ # Doxygen sources -set_property(GLOBAL APPEND PROPERTY DOXYGEN_SOURCES ${TRIQS_SOURCE_DIR}/triqs/gfs/gf.hpp) - +SET(TRIQS_GFS_SRC_DIR ${TRIQS_SOURCE_DIR}/triqs/gfs) +set_property(GLOBAL APPEND PROPERTY DOXYGEN_SOURCES + ${TRIQS_GFS_SRC_DIR}/meshes/linear.hpp + ${TRIQS_GFS_SRC_DIR}/meshes/matsubara_freq.hpp + ${TRIQS_GFS_SRC_DIR}/meshes/matsubara_time.hpp + ${TRIQS_GFS_SRC_DIR}/meshes/product.hpp + ${TRIQS_GFS_SRC_DIR}/refreq.hpp + ${TRIQS_GFS_SRC_DIR}/domains/matsubara.hpp + ) diff --git a/doc/reference/c++/gf/clef.rst b/doc/reference/c++/gf/clef.rst index f7ca3faa..f71346fd 100644 --- a/doc/reference/c++/gf/clef.rst +++ b/doc/reference/c++/gf/clef.rst @@ -5,18 +5,17 @@ Interaction with CLEF expressions ============================================ -* The gf containers and their view classes can be used with the :doc:`../clef/contents` library : +The gf containers and their view classes can be used with the :doc:`../clef/contents` library : - * They can be called with CLEF expressions. - * :doc:`Automatic assignment<../clef/assign>` has been set up. +* They can be called with CLEF expressions. +* :doc:`Automatic assignment<../clef/assign>` has been set up. -* Using the CLEF library offers a quick and efficient way to fill an array with multiple advantages : +Using the CLEF library offers a quick and efficient way to fill an array with multiple advantages : - * It is simpler and more readeable than a series of for loops. - * It is usually more optimal since the for loops are automatically written in the TraversalOrder of the array. +* It is simpler and more readeable than a series of for loops. +* It is more optimal since the loops are automatically written in the best order for memory traversal. - -* **Example** : +**Example** : .. compileblock:: @@ -36,9 +35,6 @@ Interaction with CLEF expressions .. note:: - The syntax uses a <<, not = since the array is not assigned to an expression - but filled by the evaluation thereof. - The LHS uses () and not brackets, even though it is on the mesh, because of the strange C++ limitation that [] can not be overloaded for several variables... diff --git a/doc/reference/c++/gf/concepts.rst b/doc/reference/c++/gf/concepts.rst index bf68da5b..af13b50e 100644 --- a/doc/reference/c++/gf/concepts.rst +++ b/doc/reference/c++/gf/concepts.rst @@ -37,22 +37,21 @@ Domain * **Purpose** : The domain of definition of a function. It is a mathematical definition of the domain, and does not contain any mesh, or details on its representation in a computer. -* **Refines** : RegularType, BoostSerializable, H5-serializable. +* **Refines** : RegularType. * **Definition** : -+----------------------------------------------------------------------------+---------------------------------------------------------------------+ -| Elements | Comment | -+============================================================================+=====================================================================+ -| point_t | Type of element in the domain (int, int, double, k_vector, ...) as | -| | in the call of a function over this domain. In particular, in | -| | Matsubara, it is a complex. | -+----------------------------------------------------------------------------+---------------------------------------------------------------------+ ++----------+--------------------------------------------------------------------+ +| Elements | Comment | ++==========+====================================================================+ +| point_t | Type of element in the domain (int, int, double, k_vector, ...) as | +| | in the call of a function over this domain. | ++----------+--------------------------------------------------------------------+ * **Examples** : - * Matsubara frequencies (boson/fermion) * Matsubara time + * Matsubara frequencies (boson/fermion) : in this case, point_t is `matsubara_freq`, a simple type containing (n, beta, statistics). * Real frequencies * Real time * Brillouin zone @@ -75,9 +74,7 @@ PureFunctionOnDomain +--------------------------------------+----------------------------------------------------------+ | Elements | Comment | +======================================+==========================================================+ -| domain_t | Type of the Domain represented, modelling Domain concept | -+--------------------------------------+----------------------------------------------------------+ -| domain_t const & domain() const | Returns the domain | +| domain_t const & domain() const | Returns the domain (deduced as domain_t) | +--------------------------------------+----------------------------------------------------------+ | operator (domain_t::point_t) const | Calling for all elements of the Domain (including infty | | | if it is in the domain... | @@ -99,38 +96,44 @@ Mesh It does not really need to be a mesh : e.g. if the function is represented on a polynomial basis, it is the parameters of this representation (max number of coordinates, e.g.) -* **Refines** : RegularType, HasConstIterator BoostSerializable, H5-serializable, Printable. +* **Refines** : RegularType, H5-serializable, Printable. * **Definition** : -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| Elements | Comment | -+==============================================================+===============================================================================+ -| domain_t | Type of the Domain represented, modeling the Domain concept | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| domain_t const & domain() const | Access to the domain | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| index_t | Type of indices of a point on the grid. Typically a tuple of long or a long | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| size_t size() const | The number of points in the mesh. | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| domain_t::point_t index_to_point(index_t) const | From the index of a mesh point, compute the corresponding point in the domain | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| size_t index_to_linear(index_t const &) const | Flattening the index of the mesh into a contiguous linear index | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| mesh_point_t | A type modelling MeshPoint concept (see below). | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| mesh_point_t operator[](index_t const & index ) const | From an index, return a mesh_point_t containing this a ref to this mesh and | -| | the index. | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| free function | | -| foreach ( mesh_t, lambda) | ??????????????????????????????????? | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| mesh_pt_generator iterator | A generator of all the mesh points. | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ -| const_iterator begin()/end() const | Standard access to iterator on the mesh | -| const_iterator cbegin()/cend() const | Standard access to iterator on the mesh | -+--------------------------------------------------------------+-------------------------------------------------------------------------------+ ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| Elements | Comment | ++=======================================================+===============================================================================+ +| domain_t | Type of the Domain represented, modeling the Domain concept | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| domain_t const & domain() const | Access to the domain | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| index_t | Type of indices of a point on the grid. Typically a tuple of long or a long | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| long size() const | The number of points in the mesh. | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| domain_t::point_t index_to_point(index_t) const | From the index of a mesh point, compute the corresponding point in the domain | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| long index_to_linear(index_t const &) const | Flattening the index of the mesh into a contiguous linear index | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| mesh_point_t | A type modeling MeshPoint concept (see below). | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| mesh_point_t operator[](index_t const & index ) const | From an index, return a mesh_point_t containing this a ref to this mesh and | +| | the index. | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| mesh_pt_generator const_iterator | A generator of all the mesh points. | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ +| const_iterator begin()/end() const | Standard access to iterator on the mesh Standard access to iterator on the | +| cbegin()/cend() const | mesh | ++-------------------------------------------------------+-------------------------------------------------------------------------------+ + + ++---------------------------+-----------------------------------------------------------+ +| Free functions | Comment | ++===========================+===========================================================+ +| void foreach (mesh_t, F) | If F is a function of synopsis | +| | auto F( mesh_t::mesh_point_t) | +| | it calls F for each point on the mesh, in arbitrary order | ++---------------------------+-----------------------------------------------------------+ .. _Concept_MeshPoint: @@ -193,4 +196,3 @@ one can write naturally :: // In this expression, w is casted to the domain_t::point_t, here a complex // which allows to evaluate the function - diff --git a/doc/reference/c++/gf/contents.rst b/doc/reference/c++/gf/contents.rst index 7378c68b..90008f1c 100644 --- a/doc/reference/c++/gf/contents.rst +++ b/doc/reference/c++/gf/contents.rst @@ -3,7 +3,7 @@ Green functions ================= -The TRIQS library provides a generic container `gf` and its view `gf_view`, to store and manipulate +The TRIQS library provides a generic container `gf` and its views `gf_view` and `gf_const_view`, to store and manipulate various Green functions. .. warning:: @@ -15,11 +15,11 @@ various Green functions. intro gf_and_view - gf_part_eval_curry - slicing - concepts - meshes + gf_special + clef tail fourier - gf_misc + concepts implementation_notes + + diff --git a/doc/reference/c++/gf/fourier.rst b/doc/reference/c++/gf/fourier.rst index aeaff1f4..af0bd99e 100644 --- a/doc/reference/c++/gf/fourier.rst +++ b/doc/reference/c++/gf/fourier.rst @@ -6,22 +6,64 @@ Fourier transforms The Fourier transforms from real and imaginary frequencies to times, and inverse, are currently implemented, along with the analogous transformation from the Legendre expansion to imaginary time and frequencies. -Synopsis and example -====================== +Synopsis +========= **Synopsis** :: - lazy_object fourier (gf const &) - lazy_object fourier (gf_const_view const &) + auto fourier (gf const &) + auto fourier (gf_view const &) + auto fourier (gf_const_view const &) - lazy_object inverse_fourier (gf const &) - lazy_object inverse_fourier (gf_const_view const &) + auto inverse_fourier (gf const &) + auto inverse_fourier (gf_view const &) + auto inverse_fourier (gf_const_view const &) + + gf make_gf_from_fourier(gf const&); + gf make_gf_from_fourier(gf_view const&); + gf make_gf_from_fourier(gf_const_view const&); + + gf make_gf_from_inverse_fourier(gf const&); + gf make_gf_from_inverse_fourier(gf_view const&); + gf make_gf_from_inverse_fourier(gf_const_view const&); + + +**fourier, inverse_fourier** The fourier/inverse_fourier functions do **not** perform the Fourier transformation, but returns a small lazy object (basically saying "Fourier Transform of XXX"), which is then used in an assignment of a *view* of a gf. +The reason is the following: when putting e.g. a Fourier transform of a function in time, say gt, +into a Green function in frequencies, say gw, we want to say something like:: + + gw = fourier(gt); // ??? (1) + +However, if the fourier function performs the transformation, how could it know the details +of the mesh of gw ? That information is not available when calling *fourier*. + +Since *fourier* returns a small lazy object, the library can then rewrite (1) internally into something like :: + + call_the_fourier_implementation(gt, gw); + +where all the information about the mesh of gw is now available to the implementation. + +Moreover, since fourier(gt) does not possess a domain (for the same reason), (1) +makes no sense : RHS of gf assignment requires a domain (cf concepts). +We therefore use *a view* as LHS:: + + gw() = fourier(gt); // correct usage. + +**make_gf_from_fourier, make_gf_from_inverse_fourier** + +In the case where we want to create a *new* container from the fourier transform of gt, +we can use the function make_gf_from_fourier, in which choice is explicitly made to generate a new gf with the same number of points in the mesh. +(Cf example below). + +DOC TO BE FINISHED. + Example +========= .. compileblock:: @@ -36,21 +78,14 @@ Example triqs::clef::placeholder<0> om_; gw (om_) << 1/(om_-a); - gt() = inverse_fourier(gw); // fills the *View* with the contents of the FFT. - // NB : the mesh must have the same size. + // fills a full *view* of gt with the contents of the FFT + // NB : the mesh of gt *must* have the same size as the mesh of gw. + gt() = inverse_fourier(gw); // make a new fresh gf, with the same size mesh, from the FFT of gt auto gw2 = make_gf_from_fourier(gt); } -Note that : - -* the LHS of the = must be a view, since the RHS can not compute the domain of the function - (how many points to use ?). - -* In the make_gf_from_fourier function, choice is explicitly made to generate a new gf with the same number of points in the mesh. - - Convention =========== @@ -69,7 +104,7 @@ For Matsubara (imaginary) time/frequency: The :math:`\omega_n`'s are :math:`\frac{(2n+1)\pi}{\beta}` for fermions, :math:`\frac{2n\pi}{\beta}` for bosons (as :math:`G(\tau+\beta)=-G(\tau)` for fermions, :math:`G(\tau)` for bosons). -* + .. toctree:: :maxdepth: 1 diff --git a/doc/reference/c++/gf/gf_and_view.rst b/doc/reference/c++/gf/gf_and_view.rst index 079fd1b6..9980a3a9 100644 --- a/doc/reference/c++/gf/gf_and_view.rst +++ b/doc/reference/c++/gf/gf_and_view.rst @@ -2,8 +2,8 @@ .. _gf_and_view: -gf and gf_view -================= +gf and views +============== **Synopsis**: @@ -30,6 +30,9 @@ Template parameters | Opt | void | option_t | Currently unused | +-----------------------------------------+-------------------------------+-------------------------------+--------------------------------------+ + +The various specializations of the container and its views are decribed in in the :ref:`specializations` page. + The *Variable* template parameter can take the following values : +--------------------------+--------------------------------------------+ @@ -62,39 +65,6 @@ The *Target* template parameter can take the following values : | tensor_valued | The function is tensor valued with rank | +-------------------------+-----------------------------------------------------+ -.. _gf_special: - -Specializations -------------------- - -+-------------------------+-------------------------------+-------------------------------+---------------------------+ -| Variable/Target | matrix_valued | scalar_valued | gf_valued | -+=========================+===============================+===============================+===========================+ -| 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` | -+-------------------------+-------------------------------+-------------------------------+---------------------------+ -| cartesian_product | :doc:`gf_product` | :doc:`gf_product` | | -+-------------------------+-------------------------------+-------------------------------+---------------------------+ - -.. toctree:: - :hidden: - :maxdepth: 1 - - gf_imfreq - gf_imfreq_s - gf_refreq - gf_imtime - gf_retime - gf_block - gf_product - Member types -------------------------------------- @@ -165,13 +135,15 @@ Non-member functions ------------------------ -+---------------------------------+-------------------------------------------+ -| Member function | Meaning | -+=================================+===========================================+ -| :ref:`swap` | Swap of 2 containers | -+---------------------------------+-------------------------------------------+ -| :ref:`operator\<\<` | Writing to stream | -+---------------------------------+-------------------------------------------+ ++----------------------------------------------------------------------+------------------------------------------------------------+ +| Member function | Meaning | ++======================================================================+============================================================+ +| :ref:`swap` | Swap of 2 containers | ++----------------------------------------------------------------------+------------------------------------------------------------+ +| :ref:`operator\<\<` | Writing to stream | ++----------------------------------------------------------------------+------------------------------------------------------------+ +| :ref:`reinterpret_scalar_valued_gf_as_matrix_valued` | Reinterpret a scalar_valued gf as a 1x1 matrix_valued view | ++----------------------------------------------------------------------+------------------------------------------------------------+ .. toctree:: @@ -180,16 +152,5 @@ Non-member functions stream swap -Interaction with CLEF expressions -------------------------------------- - -The gf containers and their view classes can be :doc:`used with the CLEF Library `. - -.. toctree:: - :hidden: - - clef - - diff --git a/doc/reference/c++/gf/gf_block.rst b/doc/reference/c++/gf/gf_block.rst index b168a41d..43eaa5d5 100644 --- a/doc/reference/c++/gf/gf_block.rst +++ b/doc/reference/c++/gf/gf_block.rst @@ -2,15 +2,15 @@ .. _gf_block: -block_gf (alias of gf) -=================================================== +Block Green functions +======================= This is a specialisation of :ref:`gf_and_view` for block functions. A block Green function is nothing but Green function on a discrete domain representing the -block indices, and whose value is a Green function itself. +block indices, and whose value is itself a Green function. -For convenience, for following aliases are provided :: +For convenience, the following aliases are provided :: template using block_gf = gf >; template using block_gf_view = gf_view >; diff --git a/doc/reference/c++/gf/gf_imfreq.rst b/doc/reference/c++/gf/gf_imfreq.rst index 2fc32a06..1fe8b2ca 100644 --- a/doc/reference/c++/gf/gf_imfreq.rst +++ b/doc/reference/c++/gf/gf_imfreq.rst @@ -2,39 +2,81 @@ .. _gf_imfreq: -gf +Matsubara frequencies ========================================================== +This is a specialisation of :ref:`gf` for imaginary Matsubara frequencies. -This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara frequencies. +Synopsis +------------ + +.. code:: + + gf + +The *Target* template parameter can take the following values : ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ + Domain & mesh ---------------- +The domain is :doxy:`matsubara_freq_domain`. + +The Matsubara frequencies are :math:`\omega_n=\frac{(2n+1)\pi}{\beta}` for fermions and :math:`\omega_n=\frac{2n\pi}{\beta}` for bosons. + +The mesh is :doxy:`matsubara_freq_mesh`. + Singularity ------------- -:ref:`gf_tail`. +The singularity is a high frequency expansion, :ref:`gf_tail`. -Interpolation method +Evaluation method --------------------- -None +* No interpolation. + +* Return type : + + * If Target==scalar_valued : a complex + * If Target==matrix_valued : an object modeling ImmutableMatrix concept. + +* When the point is outside of the mesh, the evaluation of the gf returns : + + * the evaluation of the high frequency tail if no_tail is not set. + * 0 otherwise + Data storage --------------- -* `data_t` : 3d array (C ordered) of complex. +* If Target==scalar_valued : + + * `data_t` : 1d array of complex. -* g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. + * g.data()(i) is the value of g for the i-th point of the mesh. + +* If Target==matrix_valued : + + * `data_t` : 3d array (C ordered) of complex. + + * g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. HDF5 storage convention --------------------------- h5 tag : `ImFreq` +TODO : DECIDE if we have 2 tag, one for scalar, one for matrix.... Examples --------- @@ -43,25 +85,37 @@ Examples .. compileblock:: #include - using namespace triqs::gfs; using triqs::clef::placeholder; + using namespace triqs::gfs; int main(){ double beta=10; int Nfreq =100; + // --- first a matrix_valued function ------------ + // First give information to build the mesh, second to build the target - auto GF1 = gf { {beta,Fermion,Nfreq}, {1,1} }; + auto g1 = gf { {beta,Fermion,Nfreq}, {1,1} }; // or a more verbose/explicit form ... - auto GF2 = gf { gf_mesh{beta,Fermion,Nfreq}, make_shape(1,1) }; + auto g2 = gf { gf_mesh{beta,Fermion,Nfreq}, make_shape(1,1) }; // Filling the gf with something... - placeholder<0> wn_; - GF1(wn_) << 1/ (wn_ + 2); + triqs::clef::placeholder<0> wn_; + g1(wn_) << 1/ (wn_ + 2); // evaluation at n=3 - std::cout << GF1(3) << " == "<< 1/ ( 1_j * std::acos(-1) / beta * (2*3+1) + 2) << std::endl; + std::cout << g1(3) << " == "<< 1/ ( 1_j * M_PI / beta * (2*3+1) + 2) << std::endl; // the high frequency expansion was automatically computed. - std::cout << GF1.singularity() << std::endl; + //std::cout << g1.singularity() << std::endl; // a bit verbose.. + + // --- a scalar_valued function ------------ + + // same a before, but without the same of the target space ... + auto g3 = gf { {beta,Fermion,Nfreq} }; + auto g4 = gf { gf_mesh{beta,Fermion,Nfreq} }; + + g3(wn_) << 1/ (wn_ + 2); + + // evaluation at n=3. + std::cout << g3(3) << " == "<< 1/ ( 1_j * std::acos(-1) / beta * (2*3+1) + 2) << std::endl; } - diff --git a/doc/reference/c++/gf/gf_imfreq_s.rst b/doc/reference/c++/gf/gf_imfreq_s.rst deleted file mode 100644 index 81650bc2..00000000 --- a/doc/reference/c++/gf/gf_imfreq_s.rst +++ /dev/null @@ -1,66 +0,0 @@ -.. 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`. - - -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; using triqs::clef::placeholder; - int main(){ - double beta=10; - int Nfreq =100; - - // First give information to build the mesh, second to build the target - auto GF1 = gf { {beta,Fermion,Nfreq} }; - // or a more verbose/explicit form ... - auto GF2 = gf { gf_mesh{beta,Fermion,Nfreq} }; - - // Filling the gf with something... - placeholder<0> wn_; - GF1(wn_) << 1/ (wn_ + 2); - - // evaluation at n=3 - std::cout << GF1(3) << " == "<< 1/ ( 1_j * std::acos(-1) / beta * (2*3+1) + 2) << std::endl; - // the high frequency expansion was automatically computed. - std::cout << GF1.singularity() << std::endl; - } - - - diff --git a/doc/reference/c++/gf/gf_imtime.rst b/doc/reference/c++/gf/gf_imtime.rst index ed92edab..ceb99123 100644 --- a/doc/reference/c++/gf/gf_imtime.rst +++ b/doc/reference/c++/gf/gf_imtime.rst @@ -2,29 +2,81 @@ .. _gf_imtime: -gf -=================================================== +Matsubara imaginary time +========================================================== -This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara time. +This is a specialisation of :ref:`gf` for imaginary Matsubara time. + +Synopsis +------------ + +.. code:: + + gf + +The *Target* template parameter can take the following values : ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ + Domain & mesh ---------------- +The domain is the set of real numbers between 0 and :math:`\beta` +since the function is periodic (resp. antiperiodic) for bosons (resp. fermions), i.e. + +* :math:`G(\tau+\beta)=-G(\tau)` for fermions +* :math:`G(\tau+\beta)=G(\tau)` for bosons. + +The domain is implemented in the class :doxy:`matsubara_time_domain`. + +The mesh is :doxy:`matsubara_time_mesh`. + + Singularity ------------- -Interpolation method +The singularity is a high frequency expansion, :ref:`gf_tail`. + + +Evaluation method --------------------- -None. +* Use a linear interpolation between the two closest point of the mesh. + +* Return type : + + * If Target==scalar_valued : a complex + * If Target==matrix_valued : an object modeling ImmutableMatrix concept. + +* When the point is outside of the mesh, the evaluation of the gf returns : + + * the evaluation of the high frequency tail if no_tail is not set. + * 0 otherwise + Data storage --------------- -* `data_t` : 3d array (C ordered) of complex. +* If Target==scalar_valued : + + * `data_t` : 1d array of complex. -* g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. + * g.data()(i) is the value of g for the i-th point of the mesh. +* If Target==matrix_valued : + + * `data_t` : 3d array (C ordered) of complex. + + * g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. + + +TO DO : complex OR DOUBLE : FIX and document !! HDF5 storage convention --------------------------- @@ -38,24 +90,38 @@ Examples .. compileblock:: #include - using namespace triqs::gfs; using triqs::clef::placeholder; + using namespace triqs::gfs; int main(){ double beta=10, a = 1; - int n_times=5; + int n_times=1000; + + // --- first a matrix_valued function ------------ // First give information to build the mesh, second to build the target - auto GF1 = gf { {beta,Fermion,n_times}, {1,1} }; + auto g1 = gf { {beta,Fermion,n_times}, {1,1} }; + // or a more verbose/explicit form ... - auto GF2 = gf { gf_mesh{beta,Fermion,n_times}, make_shape(1,1) }; - - // Filling the gf with something... - placeholder<0> tau_; - //GF1(tau_) << exp ( - a * tau_) / (1 + exp(- beta * a)); + auto g2 = gf { gf_mesh{beta,Fermion,n_times}, make_shape(1,1) }; - // evaluation at n=3 - std::cout << make_matrix(GF1(3)) << " == "<< exp ( - a * 3) / (1 + exp(- beta * a)) << std::endl; + // Filling the gf with something... COMMENT HERE : ok only because of no_tail + triqs::clef::placeholder<0> tau_; + g1(tau_) << exp ( - a * tau_) / (1 + exp(- beta * a)); + + // evaluation at tau=3.2 + std::cout << triqs::arrays::make_matrix(g1(3.2)) << " == "<< exp ( - a * 3.2) / (1 + exp(- beta * a)) << std::endl; + + // --- a scalar_valued function ------------ + + // same a before, but without the same of the target space ... + auto g3 = gf { {beta,Fermion,n_times} }; + + g3(tau_) << exp ( - a * tau_) / (1 + exp(- beta * a)); + + // evaluation at tau=3.2 + std::cout << g3(3.2) << " == "<< exp ( - a * 3.2) / (1 + exp(- beta * a)) << std::endl; } + diff --git a/doc/reference/c++/gf/gf_legendre.rst b/doc/reference/c++/gf/gf_legendre.rst new file mode 100644 index 00000000..a553ee82 --- /dev/null +++ b/doc/reference/c++/gf/gf_legendre.rst @@ -0,0 +1,66 @@ +.. highlight:: c + +.. _gf_legendre: + +Legendre representation +========================================================== + +This is a specialisation of :ref:`gf` for Legendre polynomial expansion. + +Synopsis +------------ + +.. code:: + + gf + +The *Target* template parameter can take the following values : + ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ + +Domain & mesh +---------------- + +TO BE WRITTEN + +Singularity +------------- + +None + + +Evaluation method +--------------------- + +TO BE WRITTEN + +Data storage +--------------- + +TO BE WRITTEN + +HDF5 storage convention +--------------------------- + +h5 tag : `Legendre` + +Examples +--------- + +.. compileblock:: + + #include + using namespace triqs::gfs; + + int main() { + // We want a 2x2 matrix valued function on this mesh... + //auto g = gf { {wmin, wmax, n_freq}, {2,2} }; + }; + + diff --git a/doc/reference/c++/gf/gf_product.rst b/doc/reference/c++/gf/gf_product.rst index e7a6b204..1366c90b 100644 --- a/doc/reference/c++/gf/gf_product.rst +++ b/doc/reference/c++/gf/gf_product.rst @@ -2,8 +2,8 @@ .. _gf_product: -Green functions on cartesian product domains -=================================================== +Multiple variables +=================== .. warning:: @@ -41,7 +41,18 @@ HDF5 storage convention For convenience, in hdf5 files, the arrays has higher dimension, so that the first indices are *not* flatten. EXPLAIN. +Functional techniques +------------------------ +See : + +.. toctree:: + + :maxdepth: 1 + + gf_part_eval_curry + + Examples --------- .. compileblock:: @@ -76,14 +87,5 @@ Examples } -See also : - -.. toctree:: - - :maxdepth: 1 - - gf_part_eval_curry - - diff --git a/doc/reference/c++/gf/gf_refreq.rst b/doc/reference/c++/gf/gf_refreq.rst index afc28d5e..a236cb4d 100644 --- a/doc/reference/c++/gf/gf_refreq.rst +++ b/doc/reference/c++/gf/gf_refreq.rst @@ -2,29 +2,75 @@ .. _gf_refreq: -gf -=================================================== +Real frequencies +========================================================== -This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara frequencies. +This is a specialisation of :ref:`gf` for real frequencies. + + +Synopsis +------------ + +.. code:: + + gf + +The *Target* template parameter can take the following values : ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ + Domain & mesh ---------------- +CORRECT THIS !! PUT OUT THE C++ impl for doxy doc + +The domain is :doxy:`matsubara_freq_domain`. + +The mesh is :doxy:`matsubara_freq_mesh`. + Singularity ------------- -Interpolation method +The singularity is a high frequency expansion, :ref:`gf_tail`. + + +Evaluation method --------------------- -Linear interpolation on the mesh. +* Linear interpolation on the mesh. + +* Return type : + + * If Target==scalar_valued : a complex + * If Target==matrix_valued : an object modeling ImmutableMatrix concept. + +* When the point is outside of the mesh, the evaluation of the gf returns : + + * the evaluation of the high frequency tail if no_tail is not set. + * 0 otherwise + Data storage --------------- -* `data_t` : 3d array (C ordered) of complex. +* If Target==scalar_valued : + + * `data_t` : 1d array of complex. -* g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. + * g.data()(i) is the value of g for the i-th point of the mesh. + +* If Target==matrix_valued : + + * `data_t` : 3d array (C ordered) of complex. + + * g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. HDF5 storage convention --------------------------- @@ -36,7 +82,7 @@ Examples .. compileblock:: - #include + #include using namespace triqs::gfs; int main() { diff --git a/doc/reference/c++/gf/gf_reinterpret.rst b/doc/reference/c++/gf/gf_reinterpret.rst new file mode 100644 index 00000000..73b1007a --- /dev/null +++ b/doc/reference/c++/gf/gf_reinterpret.rst @@ -0,0 +1,24 @@ +.. highlight:: c + +.. _gf_reinterpret: + +Target reinterpretation +========================= + +**Synopsis** :: + + gf_view + reinterpret_scalar_valued_gf_as_matrix_valued(gf_view); + + gf_view + reinterpret_scalar_valued_gf_as_matrix_valued(gf &); + + gf_const view + reinterpret_scalar_valued_gf_as_matrix_valued(gf_const_view); + + gf_const_view + reinterpret_scalar_valued_gf_as_matrix_valued(gf const &); + + +Given a gf or a view, scalar_valued, it returns a view of a 1x1 matrix_valued function on the same data. + diff --git a/doc/reference/c++/gf/gf_retime.rst b/doc/reference/c++/gf/gf_retime.rst index 57a0900b..3819ea12 100644 --- a/doc/reference/c++/gf/gf_retime.rst +++ b/doc/reference/c++/gf/gf_retime.rst @@ -2,29 +2,76 @@ .. _gf_retime: -gf +Real time =================================================== This is a specialisation of :ref:`gf_and_view` for imaginary Matsubara frequencies. + + +Synopsis +------------ + +.. code:: + + gf + +The *Target* template parameter can take the following values : + ++-------------------------+-----------------------------------------------------+ +| Target | Meaning | ++=========================+=====================================================+ +| scalar_valued | The function is scalar valued (double, complex...). | ++-------------------------+-----------------------------------------------------+ +| matrix_valued [default] | The function is matrix valued. | ++-------------------------+-----------------------------------------------------+ + Domain & mesh ---------------- +CORRECT THIS !! PUT OUT THE C++ impl for doxy doc + +The domain is :doxy:`matsubara_freq_domain`. + +The mesh is :doxy:`matsubara_freq_mesh`. + Singularity ------------- -Interpolation method +The singularity is a high frequency expansion, :ref:`gf_tail`. + + +Evaluation method --------------------- -Linear interpolation on the mesh. +* Linear interpolation on the mesh. + +* Return type : + + * If Target==scalar_valued : a complex + * If Target==matrix_valued : an object modeling ImmutableMatrix concept. + +* When the point is outside of the mesh, the evaluation of the gf returns : + + * the evaluation of the high frequency tail if no_tail is not set. + * 0 otherwise + Data storage --------------- -* `data_t` : 3d array (C ordered) of complex. +* If Target==scalar_valued : + + * `data_t` : 1d array of complex. -* g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. + * g.data()(i) is the value of g for the i-th point of the mesh. + +* If Target==matrix_valued : + + * `data_t` : 3d array (C ordered) of complex. + + * g.data()(i, range(), range()) is the value of g for the i-th point of the mesh. HDF5 storage convention --------------------------- @@ -36,7 +83,7 @@ Examples .. compileblock:: - #include + #include using namespace triqs::gfs; int main() { diff --git a/doc/reference/c++/gf/gf_special.rst b/doc/reference/c++/gf/gf_special.rst new file mode 100644 index 00000000..fe7e9f85 --- /dev/null +++ b/doc/reference/c++/gf/gf_special.rst @@ -0,0 +1,21 @@ +.. _gf_special: + +Specializations +------------------- + +In this section, we describe the various specializations of the generic container *gf*. +We concentrate only on specific properties, since every specialization has the members and functions described +for :ref:`gf container`. + +.. toctree:: + :maxdepth: 1 + + gf_imfreq + gf_imtime + gf_refreq + gf_retime + gf_block + gf_legendre + gf_product + + diff --git a/doc/reference/c++/gf/meshes.rst b/doc/reference/c++/gf/meshes.rst.old similarity index 85% rename from doc/reference/c++/gf/meshes.rst rename to doc/reference/c++/gf/meshes.rst.old index 418150ec..568cfc78 100644 --- a/doc/reference/c++/gf/meshes.rst +++ b/doc/reference/c++/gf/meshes.rst.old @@ -1,9 +1,10 @@ .. highlight:: c +Domains & Meshes +################## -Meshes -####### +The :doxy:`full C++ documentation` is available here. The linear meshes @@ -100,23 +101,6 @@ The domain is the set of real numbers. By default, the mesh kind is ``full_bins``. -Matsubara time ---------------- - -The domain is (approximatively) the set of real numbers between 0 and :math:`\beta`. -In fact, other points are also in the domain, but the values at these points are given by the values on this restricted domain. -:math:`G(\tau+\beta)=-G(\tau)` for fermions, :math:`G(\tau+\beta)=G(\tau)` for bosons. - -The limits from above or below at these both points can be different. -Depending on what one needs, we can choose ``full_bins``, ``half_bins`` or ``without_last``. - - -Matsubara frequency --------------------- - -The domain is discrete. The Matsubara frequencies are :math:`\omega_n=\frac{(2n+1)\pi}{\beta}` for fermions and :math:`\omega_n=\frac{2n\pi}{\beta}` for bosons. - - Products of meshes =================== diff --git a/triqs/arrays/h5/array_stack.hpp b/triqs/arrays/h5/array_stack.hpp index 349c0eaf..71b4cfb6 100644 --- a/triqs/arrays/h5/array_stack.hpp +++ b/triqs/arrays/h5/array_stack.hpp @@ -27,9 +27,7 @@ namespace triqs { namespace arrays { - /** - * The implementation class - */ + /// The implementation class template class array_stack_impl { static const size_t dim = R; static const bool base_is_array = dim > 0; diff --git a/triqs/gfs/deprecated/re_im_time.hpp b/triqs/gfs/deprecated/re_im_time.hpp index 4d313d53..ce8cf1f4 100644 --- a/triqs/gfs/deprecated/re_im_time.hpp +++ b/triqs/gfs/deprecated/re_im_time.hpp @@ -25,6 +25,7 @@ #include "../retime.hpp" #include "../imtime.hpp" #include "../meshes/product.hpp" +#include "./tool.hpp" namespace triqs { namespace gfs { diff --git a/triqs/gfs/deprecated/tool.hpp b/triqs/gfs/deprecated/tool.hpp new file mode 100644 index 00000000..b99a4c01 --- /dev/null +++ b/triqs/gfs/deprecated/tool.hpp @@ -0,0 +1,43 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * + * TRIQS is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * TRIQS. If not, see . + * + ******************************************************************************/ +#pragma once +#include "../gf.hpp" + +namespace triqs { namespace gfs { + + // make_gf and make_gf_view forward any args to them + template + gf make_gf(gf_mesh m, U &&... x) { + return gfs_implementation::factories::make_gf(std::move(m), std::forward(x)...); + } + + template + gf make_gf(U &&... x) { + return gfs_implementation::factories::make_gf(std::forward(x)...); + } + + template + gf_view make_gf_view(U &&... x) { + return gfs_implementation::factories::make_gf_view(std::forward(x)...); + } + +}} + diff --git a/triqs/gfs/deprecated/two_real_times.hpp b/triqs/gfs/deprecated/two_real_times.hpp index 3cad1917..309f9124 100644 --- a/triqs/gfs/deprecated/two_real_times.hpp +++ b/triqs/gfs/deprecated/two_real_times.hpp @@ -24,6 +24,7 @@ #include "../gf.hpp" #include "../retime.hpp" #include "../meshes/product.hpp" +#include "./tool.hpp" namespace triqs { namespace gfs { diff --git a/triqs/gfs/domains/R.hpp b/triqs/gfs/domains/R.hpp index a6f90d3c..a775a9fc 100644 --- a/triqs/gfs/domains/R.hpp +++ b/triqs/gfs/domains/R.hpp @@ -18,23 +18,22 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_DOMAIN_R_H -#define TRIQS_GF_DOMAIN_R_H +#pragma once #include "../tools.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { /// The domain struct R_domain { - typedef double point_t; - bool operator == (R_domain const & D) const { return true; } - friend void h5_write (h5::group fg, std::string subgroup_name, R_domain const & d) {} - friend void h5_read (h5::group fg, std::string subgroup_name, R_domain & d){ } + using point_t = double; + bool operator==(R_domain const& D) const { return true; } + friend void h5_write(h5::group fg, std::string subgroup_name, R_domain const& d) {} + friend void h5_read(h5::group fg, std::string subgroup_name, R_domain& d) {} friend class boost::serialization::access; - template void serialize(Archive & ar, const unsigned int version) {} + template void serialize(Archive& ar, const unsigned int version) {} }; - -}} -#endif +} +} diff --git a/triqs/gfs/domains/discrete.hpp b/triqs/gfs/domains/discrete.hpp index febc5548..a4480ec1 100644 --- a/triqs/gfs/domains/discrete.hpp +++ b/triqs/gfs/domains/discrete.hpp @@ -1,4 +1,3 @@ - /******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems @@ -19,61 +18,67 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_DISCRETE_DOMAIN_H -#define TRIQS_GF_DISCRETE_DOMAIN_H +#pragma once #include "../tools.hpp" #include -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - /// The domain class discrete_domain { - size_t Nmax; - std::vector _names;// name of the points (e.g. for block) - std::map _inv_names; - void init_inv() { - for (size_t i =0; i _names; // optional name of the points (e.g. for block) + std::map _inv_names; // table inverting the previous vector + void init_inv() { + int i = 0; + for (auto const& x : _names) _inv_names[x] = i++; + } + public: - typedef long point_t; - size_t size() const { return Nmax;}; - discrete_domain (size_t Nmax_=1) : Nmax(Nmax_) { - for (size_t i =0; i && Names) : Nmax(Names.size()), _names(Names) {init_inv(); } - discrete_domain (std::vector const & Names) : Nmax(Names.size()), _names(Names) { init_inv();} - discrete_domain (std::initializer_list const & Names) : Nmax(Names.size()), _names(Names) { init_inv();} + discrete_domain(std::vector Names) : Nmax(Names.size()), _names(std::move(Names)) { init_inv(); } + discrete_domain(std::initializer_list const& Names) : Nmax(Names.size()), _names(Names) { init_inv(); } - std::vector const & names() const { return _names;} - int index_from_name(std::string const & s) const { return _inv_names.at(s);} - - bool operator == (discrete_domain const & D) const { return (Nmax == D.Nmax);} + std::vector const& names() const { return _names; } + int index_from_name(std::string const& s) const { return _inv_names.at(s); } + + bool operator==(discrete_domain const& D) const { return (Nmax == D.Nmax); } /// Write into HDF5 - friend void h5_write (h5::group fg, std::string subgroup_name, discrete_domain const & d) { - h5::group gr = fg.create_group(subgroup_name); - h5_write(gr,"n_max",d.Nmax); + friend void h5_write(h5::group fg, std::string subgroup_name, discrete_domain const& d) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "n_max", d.Nmax); + // THE NAME ARE MISSING } /// Read from HDF5 - friend void h5_read (h5::group fg, std::string subgroup_name, discrete_domain & d){ + friend void h5_read(h5::group fg, std::string subgroup_name, discrete_domain& d) { h5::group gr = fg.open_group(subgroup_name); - long n; - h5_read(gr,"n_max",n); + int n; + h5_read(gr, "n_max", n); d = discrete_domain(n); + // NAME ARE MISSING } // BOOST Serialization friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("n_max",Nmax); - ar & boost::serialization::make_nvp("names",_names); - } - + template void serialize(Archive& ar, const unsigned int version) { + ar& boost::serialization::make_nvp("n_max", Nmax); + ar& boost::serialization::make_nvp("names", _names); + ar& boost::serialization::make_nvp("names_inv", _inv_names); + } }; - -}} -#endif +} +} diff --git a/triqs/gfs/domains/legendre.hpp b/triqs/gfs/domains/legendre.hpp index 09f7da23..dd606ef1 100644 --- a/triqs/gfs/domains/legendre.hpp +++ b/triqs/gfs/domains/legendre.hpp @@ -1,4 +1,3 @@ - /******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems @@ -19,60 +18,59 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_LEGENDRE_DOMAIN_H -#define TRIQS_GF_LEGENDRE_DOMAIN_H - +#pragma once #include "../tools.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - /// The domain class legendre_domain { public: - double beta; statistic_enum statistic; size_t Nmax; - typedef long point_t; - size_t size() const { return Nmax;}; + using point_t = long; + size_t size() const { + return Nmax; + }; - legendre_domain (double beta_ = 1, statistic_enum stat_ = Fermion, size_t Nmax_=1) : beta(beta_), statistic(stat_), Nmax(Nmax_) {} + legendre_domain(double beta_ = 1, statistic_enum stat_ = Fermion, size_t Nmax_ = 1) + : beta(beta_), statistic(stat_), Nmax(Nmax_) {} - bool operator == (legendre_domain const & D) const { - return ((std::abs(beta - D.beta)<1.e-15) && (statistic == D.statistic) && (Nmax == D.Nmax)); + bool operator==(legendre_domain const& D) const { + return ((std::abs(beta - D.beta) < 1.e-15) && (statistic == D.statistic) && (Nmax == D.Nmax)); } /// Write into HDF5 - friend void h5_write (h5::group fg, std::string subgroup_name, legendre_domain const & d) { - h5::group gr = fg.create_group(subgroup_name); - h5_write(gr,"n_max",d.Nmax); - h5_write(gr,"beta",d.beta); - h5_write(gr,"statistic",(d.statistic==Fermion ? "F" : "B")); + friend void h5_write(h5::group fg, std::string subgroup_name, legendre_domain const& d) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "n_max", d.Nmax); + h5_write(gr, "beta", d.beta); + h5_write(gr, "statistic", (d.statistic == Fermion ? "F" : "B")); } /// Read from HDF5 - friend void h5_read (h5::group fg, std::string subgroup_name, legendre_domain & d){ + friend void h5_read(h5::group fg, std::string subgroup_name, legendre_domain& d) { h5::group gr = fg.open_group(subgroup_name); - long n; double beta; std::string statistic; - h5_read(gr,"n_max",n); - h5_read(gr,"beta",beta); - h5_read(gr,"statistic",statistic); - d = legendre_domain(beta,(statistic=="F" ? Fermion : Boson),n); + long n; + double beta; + std::string statistic; + h5_read(gr, "n_max", n); + h5_read(gr, "beta", beta); + h5_read(gr, "statistic", statistic); + d = legendre_domain(beta, (statistic == "F" ? Fermion : Boson), n); } // BOOST Serialization friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("n_max",Nmax); - ar & boost::serialization::make_nvp("beta",beta); - ar & boost::serialization::make_nvp("statistic",statistic); + template void serialize(Archive& ar, const unsigned int version) { + ar& boost::serialization::make_nvp("n_max", Nmax); + ar& boost::serialization::make_nvp("beta", beta); + ar& boost::serialization::make_nvp("statistic", statistic); } - }; - -}} -#endif +} +} diff --git a/triqs/gfs/domains/matsubara.hpp b/triqs/gfs/domains/matsubara.hpp index 9157f78e..2250b5fc 100644 --- a/triqs/gfs/domains/matsubara.hpp +++ b/triqs/gfs/domains/matsubara.hpp @@ -18,26 +18,33 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_DOMAIN_MATSUBARA_H -#define TRIQS_GF_DOMAIN_MATSUBARA_H +#pragma once #include "../tools.hpp" #include namespace triqs { namespace gfs { - // --------------- One matsubara frequency, with its arithmetics ------------------------- - // all operations are done by casting to complex, except addition and substraction of matsubara_freq - + /** + * A matsubara frequency, i.e. + * * n : int, the index + * * beta : double, the temperature inverse + * * statistic : Fermion or Boson + * + * * Can be casted into a complex. + * + * * Every operations is done by casting to complex, except addition and substraction of matsubara_freq, which return matsubara_freq + * and work on the index + **/ struct matsubara_freq : public utility::arithmetic_ops_by_cast_disable_same_type> { int n; double beta; statistic_enum statistic; matsubara_freq() : n(0), beta(1), statistic(Fermion) {} - matsubara_freq(int const &n_, double beta_, statistic_enum stat_) : n(n_), beta(beta_), statistic(stat_) {} + matsubara_freq(int n_, double beta_, statistic_enum stat_) : n(n_), beta(beta_), statistic(stat_) {} using cast_t = std::complex; operator cast_t() const { - return {0, M_PI * (2 * n + (statistic == Fermion ? 1 : 0)) / beta}; + return {0, M_PI * (2 * n + statistic) / beta}; } }; @@ -50,7 +57,7 @@ namespace gfs { } inline matsubara_freq operator-(matsubara_freq const &mp) { - return {-(mp.n + mp.statistic==Fermion ? 1: 0), mp.beta, mp.statistic}; + return {-(mp.n + mp.statistic == Fermion ? 1 : 0), mp.beta, mp.statistic}; } //--------------------------------------------------------------------------------------------------------- @@ -59,12 +66,12 @@ namespace gfs { using point_t = typename std::conditional, double>::type; double beta; statistic_enum statistic; - matsubara_domain() = default; - matsubara_domain(double Beta, statistic_enum s) : beta(Beta), statistic(s) { + matsubara_domain(double beta, statistic_enum s) : beta(beta), statistic(s) { if (beta < 0) TRIQS_RUNTIME_ERROR << "Matsubara domain construction : beta <0 : beta =" << beta << "\n"; } + matsubara_domain() : matsubara_domain(1, Fermion) {} matsubara_domain(matsubara_domain const &) = default; - matsubara_domain(matsubara_domain const &x) : beta(x.beta), statistic(x.statistic) {} + matsubara_domain(matsubara_domain const &x) : matsubara_domain(x.beta, x.statistic) {} bool operator==(matsubara_domain const &D) const { return ((std::abs(beta - D.beta) < 1.e-15) && (statistic == D.statistic)); } /// Write into HDF5 @@ -91,7 +98,9 @@ namespace gfs { ar &boost::serialization::make_nvp("statistic", statistic); } }; -} -} -#endif + + using matsubara_freq_domain = matsubara_domain; + using matsubara_time_domain = matsubara_domain; +} +} diff --git a/triqs/gfs/domains/product.hpp b/triqs/gfs/domains/product.hpp index 2e4d38fd..dccc366b 100644 --- a/triqs/gfs/domains/product.hpp +++ b/triqs/gfs/domains/product.hpp @@ -18,19 +18,18 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_DOM_PRODUCT_H -#define TRIQS_GF_DOM_PRODUCT_H +#pragma once #include "../tools.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - template - struct domain_product { - typedef std::tuple point_t; - std::tuple domains; - domain_product() = default; - domain_product(Domains const & ... doms) : domains(doms...) {} - friend bool operator == (domain_product const & D1, domain_product const & D2) { return D1.domains == D2.domains;} - }; - -}} -#endif + template struct domain_product { + using point_t = std::tuple; + std::tuple domains; + domain_product() = default; + domain_product(Domains const&... doms) : domains(doms...) {} + friend bool operator==(domain_product const& D1, domain_product const& D2) { return D1.domains == D2.domains; } + // implement boost serializable, hdf5 if needed... (done at the mesh level). + }; +} +} diff --git a/triqs/gfs/evaluators.hpp b/triqs/gfs/evaluators.hpp index 2de24cc2..3bb24487 100644 --- a/triqs/gfs/evaluators.hpp +++ b/triqs/gfs/evaluators.hpp @@ -1,8 +1,8 @@ /******************************************************************************* - * + * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2013 by O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,57 +18,62 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_EVALUATOR_H -#define TRIQS_GF_EVALUATOR_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" -namespace triqs { namespace gfs { - namespace gfs_implementation { +namespace triqs { +namespace gfs { + namespace gfs_implementation { - // simple evaluation : take the point on the grid... - struct evaluator_grid_simple { - long n; - evaluator_grid_simple() = default; - - template - evaluator_grid_simple (MeshType const & m, PointType const & p) { n=p; } - template auto operator()(F const & f) const DECL_AND_RETURN(f (n)); - }; - - // a linear interpolation - struct evaluator_grid_linear_interpolation { - double w1, w2; size_t n1, n2; - - evaluator_grid_linear_interpolation() = default; - - template - evaluator_grid_linear_interpolation (MeshType const & m, PointType const & p, double prefactor=1) { - bool in; double w; - std::tie(in, n1, w) = windowing(m,p); - if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds"; - w1 = prefactor * (1-w); w2 = prefactor * w; n2 = n1 +1; - } - - template auto operator()(F const & f) const DECL_AND_RETURN(w1 * f(n1) + w2 * f (n2)); - }; - - // the evaluator for various types. - template struct evaluator_fnt_on_mesh; - - // can not use inherited constructors, too recent... -#define TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(NEWCLASS,CLASS) : CLASS { template NEWCLASS(T &&... t) : CLASS(std::forward(t)...){};}; - - // - template - struct evaluator_one_var { - public : - static constexpr int arity = 1; - evaluator_one_var() = default; - template auto operator()(G const * g, double x) const DECL_AND_RETURN(evaluator_fnt_on_mesh (g->mesh(),x)(on_mesh(*g))); - template typename G::singularity_t const & operator()(G const * g,freq_infty const &) const {return g->singularity();} + // simple evaluation : take the point on the grid... + struct evaluator_grid_simple { + long n; + evaluator_grid_simple() = default; + template evaluator_grid_simple(MeshType const &m, PointType const &p) { n = p; } + template auto operator()(F const &f) const DECL_AND_RETURN(f(n)); }; + // a linear interpolation + struct evaluator_grid_linear_interpolation { + double w1, w2; + size_t n1, n2; + + evaluator_grid_linear_interpolation() = default; + + template + evaluator_grid_linear_interpolation(MeshType const &m, PointType const &p, double prefactor = 1) { + bool in; + double w; + std::tie(in, n1, w) = windowing(m, p); + if (!in) TRIQS_RUNTIME_ERROR << " Evaluation out of bounds"; + w1 = prefactor * (1 - w); + w2 = prefactor * w; + n2 = n1 + 1; + } + + template auto operator()(F const &f) const DECL_AND_RETURN(w1 *f(n1) + w2 *f(n2)); + }; + + // the evaluator for various types. + template struct evaluator_fnt_on_mesh; + +// can not use inherited constructors, too recent... +#define TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(NEWCLASS, CLASS) : CLASS { \ + template NEWCLASS(T &&... t) : CLASS(std::forward(t)...) {}; \ + }; + + // + template struct evaluator_one_var { + public: + static constexpr int arity = 1; + evaluator_one_var() = default; + template + auto operator()(G const *g, double x) const DECL_AND_RETURN(evaluator_fnt_on_mesh(g -> mesh(), x)(on_mesh(*g))); + template typename G::singularity_t const &operator()(G const *g, freq_infty const &) const { + return g->singularity(); + } + }; } -}} -#endif +} +} diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index 5c1128e6..589e1b68 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -18,8 +18,7 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_GFBASECLASS_H -#define TRIQS_GF_GFBASECLASS_H +#pragma once #include #include #include @@ -29,7 +28,8 @@ #include "./tools.hpp" #include "./data_proxies.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { using utility::factory; using arrays::make_shape; using triqs::make_clone; @@ -48,48 +48,43 @@ namespace triqs { namespace gfs { using gf_const_view = gf_view; // the implementation class - template class gf_impl; + template class gf_impl; // various implementation traits namespace gfs_implementation { // never use using of this... // evaluator regroup functions to evaluate the function. - template struct evaluator{ static constexpr int arity = 0;}; + template struct evaluator { + static constexpr int arity = 0; + }; // closest_point mechanism - template struct get_closest_point; + template struct get_closest_point; // singularity - template struct singularity { typedef nothing type;}; + template struct singularity { + using type = nothing; + }; // symmetry - template struct symmetry { typedef nothing type;}; + template struct symmetry { + using type = nothing; + }; // data_proxy contains function to manipulate the data array, but no data itself. - // this is used to specialize this part of the code to array of dim 3 (matrix gf), dim 1 (scalar gf) and vector (e.g. block gf, ...) - template struct data_proxy; + // this is used to specialize this part of the code to array of dim 3 (matrix gf), dim 1 (scalar gf) and vector (e.g. block gf, + // ...) + template struct data_proxy; // Traits to read/write in hdf5 files. Can be specialized for some case (Cf block). Defined below - template struct h5_name; // value is a const char + template struct h5_name; // value is a const char template struct h5_rw; - + // factories regroup all factories (constructors..) for all types of gf. Defaults implemented below. template struct factories; } // gfs_implementation - // OBSOLETE : kept for backward compatibility only. Do not document. - // make_gf and make_gf_view forward any args to them - template - gf make_gf(gf_mesh m, U && ... x) - { return gfs_implementation::factories::make_gf(std::move(m),std::forward(x)...);} - - template - gf make_gf(U && ... x) { return gfs_implementation::factories::make_gf(std::forward(x)...);} - - template - gf_view make_gf_view(U && ... x) { return gfs_implementation::factories::make_gf_view(std::forward(x)...);} - // The trait that "marks" the Green function TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(ImmutableGreenFunction); @@ -102,288 +97,297 @@ namespace triqs { namespace gfs { // overload get_shape for a vector to simplify code below in gf block case. template long get_shape(std::vector const &x) { return x.size(); } - + /// A common implementation class for gf and gf_view. They will only redefine contructor and = ... template class gf_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) { static_assert(!(!IsView && IsConst), "Internal error"); - public : - typedef gf_view mutable_view_type; - typedef gf_const_view const_view_type; - typedef typename std::conditional ::type view_type; - typedef gf regular_type; + public: + using mutable_view_type = gf_view; + using const_view_type = gf_const_view; + using view_type = typename std::conditional::type; + using regular_type = gf; - typedef Variable variable_t; - typedef Target target_t; - typedef Opt option_t; + using variable_t = Variable; + using target_t = Target; + using option_t = Opt; - typedef gf_mesh mesh_t; - typedef typename mesh_t::domain_t domain_t; - typedef typename mesh_t::mesh_point_t mesh_point_t; - typedef typename mesh_t::index_t mesh_index_t; - typedef typename gfs_implementation::symmetry::type symmetry_t; - typedef gfs_implementation::evaluator evaluator_t; + using mesh_t = gf_mesh; + using domain_t = typename mesh_t::domain_t; + using mesh_point_t = typename mesh_t::mesh_point_t; + using mesh_index_t = typename mesh_t::index_t; + using symmetry_t = typename gfs_implementation::symmetry::type; + using evaluator_t = gfs_implementation::evaluator; - typedef gfs_implementation::data_proxy data_proxy_t; - typedef typename data_proxy_t::storage_t data_regular_t; - typedef typename data_proxy_t::storage_view_t data_view_t; - typedef typename data_proxy_t::storage_const_view_t data_const_view_t; - typedef typename std::conditional::type, - data_regular_t>::type data_t; + using data_proxy_t = gfs_implementation::data_proxy; + using data_regular_t = typename data_proxy_t::storage_t; + using data_view_t = typename data_proxy_t::storage_view_t; + using data_const_view_t = typename data_proxy_t::storage_const_view_t; + using data_t = typename std::conditional::type, + data_regular_t>::type; - typedef typename gfs_implementation::singularity::type singularity_non_view_t; - typedef typename view_type_if_exists_else_type::type singularity_view_t; - typedef typename std::conditional::type singularity_t; + using singularity_non_view_t = typename gfs_implementation::singularity::type; + using singularity_view_t = typename view_type_if_exists_else_type::type; + using singularity_t = typename std::conditional::type; - mesh_t const & mesh() const { return _mesh;} - domain_t const & domain() const { return _mesh.domain();} - data_t & data() { return _data;} - data_t const & data() const { return _data;} - singularity_t & singularity() { return _singularity;} - singularity_t const & singularity() const { return _singularity;} - symmetry_t const & symmetry() const { return _symmetry;} - evaluator_t const & get_evaluator() const { return _evaluator;} + mesh_t const &mesh() const { return _mesh; } + domain_t const &domain() const { return _mesh.domain(); } + data_t &data() { return _data; } + data_t const &data() const { return _data; } + singularity_t &singularity() { return _singularity; } + singularity_t const &singularity() const { return _singularity; } + symmetry_t const &symmetry() const { return _symmetry; } + evaluator_t const &get_evaluator() const { return _evaluator; } - auto get_data_shape() const DECL_AND_RETURN(get_shape(this->data())); + auto get_data_shape() const DECL_AND_RETURN(get_shape(this -> data())); - protected: - mesh_t _mesh; - data_t _data; - singularity_t _singularity; - symmetry_t _symmetry; - evaluator_t _evaluator; - data_proxy_t _data_proxy; + protected: + mesh_t _mesh; + data_t _data; + singularity_t _singularity; + symmetry_t _symmetry; + evaluator_t _evaluator; + data_proxy_t _data_proxy; - // --------------------------------Constructors ----------------------------------------------- - // all protected but one, this is an implementation class, see gf/gf_view later for public one - gf_impl() {} // all arrays of zero size (empty) + // --------------------------------Constructors ----------------------------------------------- + // all protected but one, this is an implementation class, see gf/gf_view later for public one + gf_impl() {} // all arrays of zero size (empty) - public : //everyone can make a copy and a move (for clef lib in particular, this one needs to be public) - gf_impl(gf_impl const &x) - : _mesh(x.mesh()), - _data(factory(x.data())), - _singularity(factory(x.singularity())), - _symmetry(x.symmetry()), - _evaluator(x._evaluator) {} + public: // everyone can make a copy and a move (for clef lib in particular, this one needs to be public) + gf_impl(gf_impl const &x) + : _mesh(x.mesh()), + _data(factory(x.data())), + _singularity(factory(x.singularity())), + _symmetry(x.symmetry()), + _evaluator(x._evaluator) {} - gf_impl(gf_impl &&) = default; + gf_impl(gf_impl &&) = default; - protected: - template - gf_impl(G &&x, bool) // bool to disambiguate + protected: + template + gf_impl(G &&x, bool) // bool to disambiguate : _mesh(x.mesh()), _data(factory(x.data())), _singularity(factory(x.singularity())), _symmetry(x.symmetry()), _evaluator(x.get_evaluator()) {} - template - gf_impl(M &&m, D &&dat, S &&sing, SY &&sy, EV &&ev) - : _mesh(std::forward(m)), - _data(std::forward(dat)), - _singularity(std::forward(sing)), - _symmetry(std::forward(sy)), - _evaluator(std::forward(ev)) {} + template + gf_impl(M &&m, D &&dat, S &&sing, SY &&sy, EV &&ev) + : _mesh(std::forward(m)), + _data(std::forward(dat)), + _singularity(std::forward(sing)), + _symmetry(std::forward(sy)), + _evaluator(std::forward(ev)) {} - void operator = (gf_impl const & rhs) = delete; // done in derived class. + void operator=(gf_impl const &rhs) = delete; // done in derived class. - void swap_impl (gf_impl & b) noexcept { - using std::swap; - swap(this->_mesh, b._mesh); - swap(this->_data, b._data); - swap(this->_singularity, b._singularity); - swap(this->_symmetry, b._symmetry); - swap(this->_evaluator, b._evaluator); - } + void swap_impl(gf_impl &b) noexcept { + using std::swap; + swap(this->_mesh, b._mesh); + swap(this->_data, b._data); + swap(this->_singularity, b._singularity); + swap(this->_symmetry, b._symmetry); + swap(this->_evaluator, b._evaluator); + } - public: + public: + // ------------- All the call operators ----------------------------- - // ------------- All the call operators ----------------------------- + // First, a simple () returns a view, like for an array... + const_view_type operator()() const { return *this; } + view_type operator()() { return *this; } - // First, a simple () returns a view, like for an array... - const_view_type operator()() const { return *this; } - view_type operator()() { return *this; } + /// Calls are (perfectly) forwarded to the evaluator::operator(), except mesh_point_t and when + /// there is at least one lazy argument ... + template // match any argument list, picking out the first type : () is not permitted + typename std::add_const::value || + ((sizeof...(Args) != evaluator_t::arity) && (evaluator_t::arity != -1)) // if -1 : no check + , + std::result_of // what is the result type of call + >::type // end of lazy_disable_if + >::type // end of add_Const + operator()(Args &&... args) const { + return _evaluator(this, std::forward(args)...); + } - /// Calls are (perfectly) forwarded to the evaluator::operator(), except mesh_point_t and when - /// there is at least one lazy argument ... - template // match any argument list, picking out the first type : () is not permitted - typename std::add_const::value || - ((sizeof...(Args) != evaluator_t::arity) && (evaluator_t::arity != -1)) // if -1 : no check - , - std::result_of // what is the result type of call - >::type // end of lazy_disable_if - >::type // end of add_Const - operator()(Args &&... args) const { - return _evaluator(this, std::forward(args)...); - } + template typename clef::_result_of::make_expr_call::type operator()(Args &&... args) & { + return clef::make_expr_call(*this, std::forward(args)...); + } - template - typename clef::_result_of::make_expr_call::type operator()(Args &&... args) & { - return clef::make_expr_call(*this, std::forward(args)...); - } + template + typename clef::_result_of::make_expr_call::type operator()(Args &&... args) const &{ + return clef::make_expr_call(*this, std::forward(args)...); + } - template - typename clef::_result_of::make_expr_call::type operator()(Args &&... args) const &{ - return clef::make_expr_call(*this, std::forward(args)...); - } + template typename clef::_result_of::make_expr_call::type operator()(Args &&... args) && { + return clef::make_expr_call(std::move(*this), std::forward(args)...); + } - template typename clef::_result_of::make_expr_call::type operator()(Args &&... args) && { - return clef::make_expr_call(std::move(*this), std::forward(args)...); - } + // ------------- All the [] operators ----------------------------- + // [] and access to the grid point + using r_type = typename std::result_of::type; + using cr_type = typename std::result_of::type; - // ------------- All the [] operators ----------------------------- - // [] and access to the grid point - typedef typename std::result_of::type r_type; - typedef typename std::result_of::type cr_type; + r_type operator[](mesh_index_t const &arg) { return _data_proxy(_data, _mesh.index_to_linear(arg)); } + cr_type operator[](mesh_index_t const &arg) const { return _data_proxy(_data, _mesh.index_to_linear(arg)); } - r_type operator[] (mesh_index_t const & arg) { return _data_proxy(_data,_mesh.index_to_linear(arg));} - cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));} + r_type operator[](mesh_point_t const &x) { return _data_proxy(_data, x.linear_index()); } + cr_type operator[](mesh_point_t const &x) const { return _data_proxy(_data, x.linear_index()); } - r_type operator[] (mesh_point_t const & x) { return _data_proxy(_data, x.linear_index());} - cr_type operator[] (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());} + template r_type operator[](closest_pt_wrap const &p) { + return _data_proxy(_data, + _mesh.index_to_linear(gfs_implementation::get_closest_point::invoke(this, p))); + } + template cr_type operator[](closest_pt_wrap const &p) const { + return _data_proxy(_data, + _mesh.index_to_linear(gfs_implementation::get_closest_point::invoke(this, p))); + } - template - r_type operator[] (closest_pt_wrap const & p) { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point::invoke(this,p)));} - template - cr_type operator[] (closest_pt_wrap const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point::invoke(this,p)));} + template + typename clef::_result_of::make_expr_subscript::type operator[](Arg &&arg) const &{ + return clef::make_expr_subscript(*this, std::forward(arg)); + } - template - typename clef::_result_of::make_expr_subscript::type - operator[](Arg && arg) const & { return clef::make_expr_subscript(*this,std::forward(arg));} + template typename clef::_result_of::make_expr_subscript::type operator[](Arg &&arg) & { + return clef::make_expr_subscript(*this, std::forward(arg)); + } - template - typename clef::_result_of::make_expr_subscript::type - operator[](Arg && arg) & { return clef::make_expr_subscript(*this,std::forward(arg));} + template typename clef::_result_of::make_expr_subscript::type operator[](Arg &&arg) && { + return clef::make_expr_subscript(std::move(*this), std::forward(arg)); + } - template - typename clef::_result_of::make_expr_subscript::type - operator[](Arg && arg) && { return clef::make_expr_subscript(std::move(*this),std::forward(arg));} + /// --------------------- A direct access to the grid point -------------------------- + template r_type on_mesh(Args &&... args) { + return _data_proxy(_data, _mesh.index_to_linear(mesh_index_t(std::forward(args)...))); + } - /// --------------------- A direct access to the grid point -------------------------- - template - r_type on_mesh (Args&&... args) { return _data_proxy(_data,_mesh.index_to_linear(mesh_index_t(std::forward(args)...)));} + template cr_type on_mesh(Args &&... args) const { + return _data_proxy(_data, _mesh.index_to_linear(mesh_index_t(std::forward(args)...))); + } - template - cr_type on_mesh (Args&&... args) const { return _data_proxy(_data,_mesh.index_to_linear(mesh_index_t(std::forward(args)...)));} + // The on_mesh little adaptor .... + private: + struct _on_mesh_wrapper_const { + gf_impl const &f; + template cr_type operator()(Args &&... args) const { return f.on_mesh(std::forward(args)...); } + }; + struct _on_mesh_wrapper { + gf_impl &f; + template r_type operator()(Args &&... args) const { return f.on_mesh(std::forward(args)...); } + }; - // The on_mesh little adaptor .... - private: - struct _on_mesh_wrapper_const { - gf_impl const &f; - template cr_type operator()(Args &&... args) const { return f.on_mesh(std::forward(args)...); } - }; - struct _on_mesh_wrapper { - gf_impl &f; - template r_type operator()(Args &&... args) const { return f.on_mesh(std::forward(args)...); } - }; + public: + _on_mesh_wrapper_const friend on_mesh(gf_impl const &f) { + return {f}; + } + _on_mesh_wrapper friend on_mesh(gf_impl &f) { + return {f}; + } - public: - _on_mesh_wrapper_const friend on_mesh(gf_impl const &f) { - return {f}; - } - _on_mesh_wrapper friend on_mesh(gf_impl &f) { - return {f}; - } + //----------------------------- HDF5 ----------------------------- - //----------------------------- HDF5 ----------------------------- + friend std::string get_triqs_hdf5_data_scheme(gf_impl const &g) { + return "Gf" + gfs_implementation::h5_name::invoke(); + } - friend std::string get_triqs_hdf5_data_scheme(gf_impl const & g) { return "Gf" + gfs_implementation::h5_name::invoke();} + friend class gfs_implementation::h5_rw; - friend class gfs_implementation::h5_rw; + /// Write into HDF5 + friend void h5_write(h5::group fg, std::string subgroup_name, gf_impl const &g) { + auto gr = fg.create_group(subgroup_name); + gr.write_triqs_hdf5_data_scheme(g); + gfs_implementation::h5_rw::write(gr, g); + } - /// Write into HDF5 - friend void h5_write (h5::group fg, std::string subgroup_name, gf_impl const & g) { - auto gr = fg.create_group(subgroup_name); - gr.write_triqs_hdf5_data_scheme(g); - gfs_implementation::h5_rw::write(gr, g); - } + /// Read from HDF5 + friend void h5_read(h5::group fg, std::string subgroup_name, gf_impl &g) { + auto gr = fg.open_group(subgroup_name); + // Check the attribute or throw + auto tag_file = gr.read_triqs_hdf5_data_scheme(); + auto tag_expected = get_triqs_hdf5_data_scheme(g); + if (tag_file != tag_expected) + TRIQS_RUNTIME_ERROR << "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found " << tag_file + << " while I expected " << tag_expected; + gfs_implementation::h5_rw::read(gr, g); + } - /// Read from HDF5 - friend void h5_read (h5::group fg, std::string subgroup_name, gf_impl & g) { - auto gr = fg.open_group(subgroup_name); - // Check the attribute or throw - auto tag_file = gr.read_triqs_hdf5_data_scheme(); - auto tag_expected= get_triqs_hdf5_data_scheme(g); - if (tag_file != tag_expected) - TRIQS_RUNTIME_ERROR<< "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found "<::read(gr, g); - } - - //----------------------------- BOOST Serialization ----------------------------- - friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("data",_data); - ar & boost::serialization::make_nvp("singularity",_singularity); - ar & boost::serialization::make_nvp("mesh",_mesh); - ar & boost::serialization::make_nvp("symmetry",_symmetry); - } - - /// print - friend std::ostream & operator << (std::ostream & out, gf_impl const & x) { return out<<(IsView ? "gf_view": "gf");} - friend std::ostream & triqs_nvl_formal_print(std::ostream & out, gf_impl const & x) { return out<<(IsView ? "gf_view": "gf");} + //----------------------------- BOOST Serialization ----------------------------- + friend class boost::serialization::access; + template void serialize(Archive &ar, const unsigned int version) { + ar &boost::serialization::make_nvp("data", _data); + ar &boost::serialization::make_nvp("singularity", _singularity); + ar &boost::serialization::make_nvp("mesh", _mesh); + ar &boost::serialization::make_nvp("symmetry", _symmetry); + } + /// print + friend std::ostream &operator<<(std::ostream &out, gf_impl const &x) { return out << (IsView ? "gf_view" : "gf"); } + friend std::ostream &triqs_nvl_formal_print(std::ostream &out, gf_impl const &x) { return out << (IsView ? "gf_view" : "gf"); } }; - // -------------------------Interaction with the CLEF library : auto assignement implemnetation-------------------------------------- + // -------------------------Interaction with the CLEF library : auto assignement implementation ----------------- // auto assignment of the gf (gf(om_) << expression fills the functions by evaluation of expression) - template - void triqs_clef_auto_assign(gf_impl &g, RHS const &rhs) { - triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of>::type()); - assign_from_expression(g.singularity(), rhs); - // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !) - // if f is an expression, replace the placeholder with a simple tail. If f is a function callable on freq_infty, - // it uses the fact that tail_non_view_t can be casted into freq_infty - } + template + void triqs_clef_auto_assign(gf_impl &g, RHS const &rhs) { + triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of>::type()); + assign_from_expression(g.singularity(), rhs); + // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !) + // if f is an expression, replace the placeholder with a simple tail. If f is a function callable on freq_infty, + // it uses the fact that tail_non_view_t can be casted into freq_infty + } - // enable the writing g[om_] << .... also - template - void triqs_clef_auto_assign_subscript(gf_impl &g, RHS const &rhs) { - triqs_clef_auto_assign(g, rhs); - } + // enable the writing g[om_] << .... also + template + void triqs_clef_auto_assign_subscript(gf_impl &g, RHS const &rhs) { + triqs_clef_auto_assign(g, rhs); + } - template - void triqs_gf_clef_auto_assign_impl_aux_assign(G &&g, RHS &&rhs, std::integral_constant) { - std::forward(g) = std::forward(rhs); - } + template + void triqs_gf_clef_auto_assign_impl_aux_assign(G &&g, RHS &&rhs, std::integral_constant) { + std::forward(g) = std::forward(rhs); + } - template - void triqs_gf_clef_auto_assign_impl_aux_assign(G &&g, clef::make_fun_impl &&rhs, std::integral_constant) { - triqs_clef_auto_assign_impl(std::forward(g), std::forward>(rhs), std::integral_constant()); - } + template + void triqs_gf_clef_auto_assign_impl_aux_assign(G &&g, clef::make_fun_impl &&rhs, std::integral_constant) { + triqs_clef_auto_assign_impl(std::forward(g), std::forward>(rhs), + std::integral_constant()); + } - template - void triqs_clef_auto_assign_impl(gf_impl &g, RHS const &rhs, - std::integral_constant) { - for (auto const &w : g.mesh()) { - triqs_gf_clef_auto_assign_impl_aux_assign(g[w], rhs(w), std::integral_constant()); - //(*this)[w] = rhs(w); - } + template + void triqs_clef_auto_assign_impl(gf_impl &g, RHS const &rhs, + std::integral_constant) { + for (auto const &w : g.mesh()) { + triqs_gf_clef_auto_assign_impl_aux_assign(g[w], rhs(w), std::integral_constant()); + //(*this)[w] = rhs(w); } + } - template - void triqs_clef_auto_assign_impl(gf_impl &g, RHS const &rhs, - std::integral_constant) { - for (auto const &w : g.mesh()) { - triqs_gf_clef_auto_assign_impl_aux_assign(g[w], triqs::tuple::apply(rhs, w.components_tuple()), std::integral_constant()); - //(*this)[w] = triqs::tuple::apply(rhs, w.components_tuple()); - } + template + void triqs_clef_auto_assign_impl(gf_impl &g, RHS const &rhs, + std::integral_constant) { + for (auto const &w : g.mesh()) { + triqs_gf_clef_auto_assign_impl_aux_assign(g[w], triqs::tuple::apply(rhs, w.components_tuple()), + std::integral_constant()); + //(*this)[w] = triqs::tuple::apply(rhs, w.components_tuple()); } + } // -------------------------The regular class of GF -------------------------------------------------------- - - template class gf : public gf_impl { - typedef gf_impl B; - typedef gfs_implementation::factories factory; - public : + + template class gf : public gf_impl { + using B = gf_impl; + using factory = gfs_implementation::factories; + + public: gf() : B() {} gf(gf const &g) : B(g) {} gf(gf &&g) noexcept : B(std::move(g)) {} - gf(gf_view const &g) : B(g, bool{}) {} - gf(gf_const_view const &g) : B(g, bool{}) {} + gf(gf_view const &g) : B(g, bool {}) {} + gf(gf_const_view const &g) : B(g, bool {}) {} template gf(GfType const &x, typename std::enable_if::value>::type *dummy = 0) @@ -394,50 +398,64 @@ namespace triqs { namespace gfs { gf(typename B::mesh_t m, typename B::data_t dat, typename B::singularity_view_t const &si, typename B::symmetry_t const &s) : B(std::move(m), std::move(dat), si, s, typename B::evaluator_t{}) {} - typedef typename factory::target_shape_t target_shape_t; + using target_shape_t = typename factory::target_shape_t; - gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{} ): - B(std::move(m), factory::make_data(m,shape), factory::make_singularity(m,shape), typename B::symmetry_t {}, typename B::evaluator_t{}) {} + gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}) + : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, + typename B::evaluator_t{}) {} - friend void swap (gf & a, gf & b) noexcept { a.swap_impl (b);} + friend void swap(gf &a, gf &b) noexcept { a.swap_impl(b); } - gf & operator = (gf const & rhs) { *this = gf(rhs); return *this;} // use move = - gf & operator = (gf & rhs) { *this = gf(rhs); return *this;} // use move = - gf & operator = (gf && rhs) noexcept { swap(*this,rhs); return *this;} - - template void operator = (RHS && rhs) { - this->_mesh = rhs.mesh(); - this->_data.resize(get_gf_data_shape(rhs)); - for (auto const & w: this->mesh()) { (*this)[w] = rhs[w]; } - this->_singularity = rhs.singularity(); - // to be implemented : there is none in the gf_expr in particular.... - //this->_symmetry = rhs.symmetry(); + gf &operator=(gf const &rhs) { + *this = gf(rhs); + return *this; + } // use move = + gf &operator=(gf &rhs) { + *this = gf(rhs); + return *this; + } // use move = + gf &operator=(gf &&rhs) noexcept { + swap(*this, rhs); + return *this; } + template void operator=(RHS &&rhs) { + this->_mesh = rhs.mesh(); + this->_data.resize(get_gf_data_shape(rhs)); + for (auto const &w : this->mesh()) { + (*this)[w] = rhs[w]; + } + this->_singularity = rhs.singularity(); + // to be implemented : there is none in the gf_expr in particular.... + // this->_symmetry = rhs.symmetry(); + } }; // --------------------------The const View class of GF ------------------------------------------------------- - - template class gf_view : public gf_impl { - typedef gf_impl B; - public : + + template + class gf_view : public gf_impl { + using B = gf_impl; + + public: gf_view() = delete; gf_view(gf_view const &g) : B(g) {} gf_view(gf_view &&g) noexcept : B(std::move(g)) {} - gf_view(gf_impl const &g) : B(g, bool{}) {} // from a const_view - gf_view(gf_impl const &g) : B(g, bool{}) {} // from a view - gf_view(gf_impl const &g) : B(g, bool{}) {} // from a const gf - gf_view(gf_impl &g) : B(g, bool{}) {} // from a gf & - gf_view(gf_impl &&g) noexcept : B(std::move(g), bool{}) {} // from a gf && + gf_view(gf_impl const &g) : B(g, bool {}) {} // from a const_view + gf_view(gf_impl const &g) : B(g, bool {}) {} // from a view + gf_view(gf_impl const &g) : B(g, bool {}) {} // from a const gf + gf_view(gf_impl &g) : B(g, bool {}) {} // from a gf & + gf_view(gf_impl &&g) noexcept : B(std::move(g), bool {}) {} // from a gf && template gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s) : B(m, factory(dat), t, s, typename B::evaluator_t{}) {} void rebind(gf_view const &X) noexcept { - this->_mesh = X._mesh; this->_symmetry = X._symmetry; - this->_data_proxy.rebind(this->_data,X); + this->_mesh = X._mesh; + this->_symmetry = X._symmetry; + this->_data_proxy.rebind(this->_data, X); this->_singularity.rebind(X._singularity); } @@ -445,33 +463,36 @@ namespace triqs { namespace gfs { rebind(gf_view{X}); } - gf_view & operator = (gf_view const & ) = delete; + gf_view &operator=(gf_view const &) = delete; }; // class gf_const_view // ------------------------- The View class of GF ------------------------------------------------------- - - template class gf_view : public gf_impl { - typedef gf_impl B; - public : + + template + class gf_view : public gf_impl { + using B = gf_impl; + + public: gf_view() = delete; gf_view(gf_view const &g) : B(g) {} gf_view(gf_view &&g) noexcept : B(std::move(g)) {} - gf_view(gf_impl const &g) = delete; // from a const view : impossible - gf_view(gf_impl const &g) : B(g, bool{}) {} // from a view - gf_view(gf_impl const &g) = delete; // from a const gf : impossible - gf_view(gf_impl &g) : B(g, bool{}) {} // from a gf & - gf_view(gf_impl &&g) noexcept : B(std::move(g), bool{}) {} // from a gf && + gf_view(gf_impl const &g) = delete; // from a const view : impossible + gf_view(gf_impl const &g) : B(g, bool {}) {} // from a view + gf_view(gf_impl const &g) = delete; // from a const gf : impossible + gf_view(gf_impl &g) : B(g, bool {}) {} // from a gf & + gf_view(gf_impl &&g) noexcept : B(std::move(g), bool {}) {} // from a gf && template gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s) : B(m, factory(dat), t, s, typename B::evaluator_t{}) {} - friend void swap (gf_view & a, gf_view & b) noexcept { a.swap_impl (b);} + friend void swap(gf_view &a, gf_view &b) noexcept { a.swap_impl(b); } void rebind(gf_view const &X) noexcept { - this->_mesh = X._mesh; this->_symmetry = X._symmetry; - this->_data_proxy.rebind(this->_data,X); + this->_mesh = X._mesh; + this->_symmetry = X._symmetry; + this->_data_proxy.rebind(this->_data, X); this->_singularity.rebind(X._singularity); } @@ -487,89 +508,93 @@ namespace triqs { namespace gfs { }; // class gf_view // delegate = so that I can overload it for specific RHS... - template - DISABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation( gf_view g, RHS const & rhs) { - if (!(g.mesh() == rhs.mesh())) TRIQS_RUNTIME_ERROR<<"Gf Assignment in View : incompatible mesh"< + DISABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation(gf_view g, RHS const &rhs) { + if (!(g.mesh() == rhs.mesh())) + TRIQS_RUNTIME_ERROR << "Gf Assignment in View : incompatible mesh" << g.mesh() << " vs " << rhs.mesh(); + for (auto const &w : g.mesh()) g[w] = rhs[w]; + g.singularity() = rhs.singularity(); + } - template - ENABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation( gf_view g, T const & x) { - gf_view::data_proxy_t::assign_to_scalar(g.data(), x); - g.singularity() = x; - } + template + ENABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation(gf_view g, T const &x) { + gf_view::data_proxy_t::assign_to_scalar(g.data(), x); + g.singularity() = x; + } // tool for lazy transformation - template struct gf_keeper { - gf_const_view g; - }; + template struct gf_keeper { + gf_const_view g; + }; // ---------------------------------- slicing ------------------------------------ - // slice - template - gf_view slice_target(gf_view g, Args &&... args) { - static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); - using arrays::range; - return {g.mesh(), g.data()(range(), std::forward(args)...), - slice_target(g.singularity(), std::forward(args)...), g.symmetry()}; - } + // slice + template + gf_view slice_target(gf_view g, Args &&... args) { + static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); + using arrays::range; + return {g.mesh(), g.data()(range(), std::forward(args)...), + slice_target(g.singularity(), std::forward(args)...), g.symmetry()}; + } - template - gf_view slice_target(gf &g, Args &&... args) { - return slice_target(g(), std::forward(args)...); - } + template + gf_view slice_target(gf &g, Args &&... args) { + return slice_target(g(), std::forward(args)...); + } - template - gf_const_view slice_target(gf const &g, Args &&... args) { - return slice_target(g(), std::forward(args)...); - } + template + gf_const_view slice_target(gf const &g, Args &&... args) { + return slice_target(g(), std::forward(args)...); + } - // slice to scalar - template - gf_view slice_target_to_scalar(gf_view g, - Args &&... args) { - static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); - using arrays::range; - return {g.mesh(), g.data()(range(), std::forward(args)...), - slice_target(g.singularity(), range(args, args + 1)...), g.symmetry()}; - } + // slice to scalar + template + gf_view slice_target_to_scalar(gf_view g, + Args &&... args) { + static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); + using arrays::range; + return {g.mesh(), g.data()(range(), std::forward(args)...), + slice_target(g.singularity(), range(args, args + 1)...), g.symmetry()}; + } - template - gf_view slice_target_to_scalar(gf &g, Args &&... args) { - return slice_target_to_scalar(g(), std::forward(args)...); - } + template + gf_view slice_target_to_scalar(gf &g, Args &&... args) { + return slice_target_to_scalar(g(), std::forward(args)...); + } - template - gf_const_view slice_target_to_scalar(gf const &g, Args &&... args) { - return slice_target_to_scalar(g(), std::forward(args)...); - } + template + gf_const_view slice_target_to_scalar(gf const &g, Args &&... args) { + return slice_target_to_scalar(g(), std::forward(args)...); + } - // a scalar_valued gf can be viewed as a 1x1 matrix - template - gf_view - reinterpret_scalar_valued_gf_as_matrix_valued(gf_view g) { - typedef typename gf_view::data_view_t a_t; - auto a = a_t{typename a_t::indexmap_type(join(g.data().shape(), make_shape(1, 1))), g.data().storage()}; - return {g.mesh(), a, g.singularity(), g.symmetry()}; - } + // ---------------------------------- target reinterpretation ------------------------------------ + + // a scalar_valued gf can be viewed as a 1x1 matrix + template + gf_view + reinterpret_scalar_valued_gf_as_matrix_valued(gf_view g) { + using a_t = typename gf_view::data_view_t; + auto a = a_t{typename a_t::indexmap_type(join(g.data().shape(), make_shape(1, 1))), g.data().storage()}; + return {g.mesh(), a, g.singularity(), g.symmetry()}; + } - template - gf_view reinterpret_scalar_valued_gf_as_matrix_valued(gf &g) { - return reinterpret_scalar_valued_gf_as_matrix_valued(g()); - } + template + gf_view reinterpret_scalar_valued_gf_as_matrix_valued(gf &g) { + return reinterpret_scalar_valued_gf_as_matrix_valued(g()); + } - template - gf_const_view - reinterpret_scalar_valued_gf_as_matrix_valued(gf const &g) { - return reinterpret_scalar_valued_gf_as_matrix_valued(g()); - } + template + gf_const_view + reinterpret_scalar_valued_gf_as_matrix_valued(gf const &g) { + return reinterpret_scalar_valued_gf_as_matrix_valued(g()); + } /* template gf_view slice_mesh (gf_impl const & g, Args... args) { - return gf_view(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry()); + return gf_view(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), + g.singularity(), g.symmetry()); }*/ namespace gfs_implementation { // implement some default traits @@ -578,9 +603,9 @@ namespace triqs { namespace gfs { // ----- tensor_valued template struct factories, Opt> { - typedef gf, Opt> gf_t; - typedef tqa::mini_vector target_shape_t; - typedef typename gf_t::mesh_t mesh_t; + using gf_t = gf, Opt>; + using target_shape_t = arrays::mini_vector; + using mesh_t = typename gf_t::mesh_t; static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) { typename gf_t::data_t A(shape.front_append(m.size())); @@ -595,9 +620,9 @@ namespace triqs { namespace gfs { // ----- matrix_valued template struct factories { - typedef gf gf_t; - typedef tqa::mini_vector target_shape_t; - typedef typename gf_t::mesh_t mesh_t; + using gf_t = gf; + using target_shape_t = arrays::mini_vector; + using mesh_t = typename gf_t::mesh_t; static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) { typename gf_t::data_t A(shape.front_append(m.size())); @@ -612,7 +637,7 @@ namespace triqs { namespace gfs { // ----- scalar_valued template struct factories { - typedef gf gf_t; + using gf_t = gf; struct target_shape_t { target_shape_t front_pop() const { // this make the get_target_shape function works in this case... return {}; @@ -620,8 +645,8 @@ namespace triqs { namespace gfs { target_shape_t() = default; template target_shape_t(utility::mini_vector) {} }; - - typedef typename gf_t::mesh_t mesh_t; + + using mesh_t = typename gf_t::mesh_t; static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) { typename gf_t::data_t A(m.size()); @@ -633,7 +658,7 @@ namespace triqs { namespace gfs { } }; - // --------------------- hdf5 --------------------------------------- + // --------------------- hdf5 --------------------------------------- // scalar function : just add a _s template struct h5_name { static std::string invoke() { return h5_name::invoke() + "_s"; } @@ -657,14 +682,12 @@ namespace triqs { namespace gfs { } }; } // gfs_implementation - - -}} +} +} // same as for arrays : views can not be swapped by the std::swap. Delete it namespace std { - template - void swap(triqs::gfs::gf_view &a, triqs::gfs::gf_view &b) = delete; +template +void swap(triqs::gfs::gf_view &a, triqs::gfs::gf_view &b) = delete; } #include "./gf_expr.hpp" -#endif diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index ecc1f725..604f4e04 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,8 +18,7 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_MATSUBARA_FREQ_H -#define TRIQS_GF_MATSUBARA_FREQ_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" @@ -32,22 +31,18 @@ namespace gfs { struct imfreq {}; template struct gf_mesh : matsubara_freq_mesh { - using B = matsubara_freq_mesh; - static double m1(double beta) { return std::acos(-1) / beta; } - gf_mesh() = default; - gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt - gf_mesh(typename B::domain_t const &d, int Nmax = 1025, bool positive_only = true) : B(d, Nmax, positive_only) {} - gf_mesh(double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta, S}, Nmax) {} + template gf_mesh(T &&... x) : matsubara_freq_mesh(std::forward(x)...) {} + //using matsubara_freq_mesh::matsubara_freq_mesh; }; namespace gfs_implementation { // singularity template <> struct singularity { - typedef local::tail type; + using type = local::tail; }; template <> struct singularity { - typedef local::tail type; + using type = local::tail; }; // h5 name @@ -124,10 +119,13 @@ namespace gfs { #ifdef __clang__ // to generate a clearer error message ? . Only ok on clang ? - template struct error { static_assert(n>0, "Green function can not be evaluated on a complex number !");}; + template struct error { + static_assert(n > 0, "Green function can not be evaluated on a complex number !"); + }; - template - error<0> operator()(G const *g, std::complex) const { return {};} + template error<0> operator()(G const *g, std::complex) const { + return {}; + } #endif template typename G::singularity_t const &operator()(G const *g, freq_infty const &) const { @@ -142,4 +140,3 @@ namespace gfs { } // gfs_implementation } } -#endif diff --git a/triqs/gfs/imtime.hpp b/triqs/gfs/imtime.hpp index a6447aa9..a8bab6b7 100644 --- a/triqs/gfs/imtime.hpp +++ b/triqs/gfs/imtime.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,83 +18,89 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_MATSUBARA_TIME_H -#define TRIQS_GF_MATSUBARA_TIME_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" #include "./local/no_tail.hpp" #include "./domains/matsubara.hpp" -#include "./meshes/linear.hpp" +#include "./meshes/matsubara_time.hpp" #include "./evaluators.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { struct imtime {}; - // gf_mesh type and its factories - template struct gf_mesh : linear_mesh> { - typedef linear_mesh> B; - gf_mesh() = default; - gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt - gf_mesh(typename B::domain_t d, int n_time_slices, mesh_kind mk = half_bins) : B(d, 0, d.beta, n_time_slices, mk) {} - gf_mesh(double beta, statistic_enum S, int n_time_slices, mesh_kind mk = half_bins) : gf_mesh({beta, S}, n_time_slices, mk) {} + template struct gf_mesh : matsubara_time_mesh { + template gf_mesh(T &&... x) : matsubara_time_mesh(std::forward(x)...) {} + // using matsubara_time_mesh::matsubara_time_mesh; }; namespace gfs_implementation { // singularity. If no_tail is given, then it is the default (nothing) - template<> struct singularity { typedef local::tail type;}; - template<> struct singularity { typedef local::tail type;}; + template <> struct singularity { + using type = local::tail; + }; + template <> struct singularity { + using type = local::tail; + }; // h5 name - template struct h5_name { static std::string invoke(){ return "ImTime";}}; + template struct h5_name { + static std::string invoke() { return "ImTime"; } + }; /// --------------------------- data access --------------------------------- - template struct data_proxy : data_proxy_array {}; - template struct data_proxy : data_proxy_array {}; + template struct data_proxy : data_proxy_array {}; + template struct data_proxy : data_proxy_array {}; /// --------------------------- closest mesh point on the grid --------------------------------- - template - struct get_closest_point { - // index_t is int - template - static int invoke(G const * g, closest_pt_wrap const & p) { - double x = (g->mesh().kind()==half_bins ? double(p.value) : double(p.value)+ 0.5*g->mesh().delta()); - int n = std::floor(x/g->mesh().delta()); - return n; - } - }; + template struct get_closest_point { + // index_t is int + template static int invoke(G const *g, closest_pt_wrap const &p) { + double x = (g->mesh().kind() == half_bins ? double(p.value) : double(p.value) + 0.5 * g->mesh().delta()); + int n = std::floor(x / g->mesh().delta()); + return n; + } + }; /// --------------------------- evaluator --------------------------------- // this one is specific because of the beta-antiperiodicity for fermions - template<> - struct evaluator_fnt_on_mesh { - double w1, w2; long n; + template <> struct evaluator_fnt_on_mesh { + double w1, w2; + long n; - evaluator_fnt_on_mesh() = default; + evaluator_fnt_on_mesh() = default; - evaluator_fnt_on_mesh (gf_mesh const & m, double tau) { - double beta = m.domain().beta; - int p = std::floor(tau/beta); - tau -= p*beta; - double w; bool in; - std::tie(in, n, w) = windowing(m,tau); - if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds"; - if ((m.domain().statistic == Fermion) && (p%2==1)) {w2 = -w; w1 = w-1;} else { w2 = w; w1 = 1-w;} + evaluator_fnt_on_mesh(gf_mesh const &m, double tau) { + double beta = m.domain().beta; + int p = std::floor(tau / beta); + tau -= p * beta; + double w; + bool in; + std::tie(in, n, w) = windowing(m, tau); + if (!in) TRIQS_RUNTIME_ERROR << " Evaluation out of bounds"; + if ((m.domain().statistic == Fermion) && (p % 2 == 1)) { + w2 = -w; + w1 = w - 1; + } else { + w2 = w; + w1 = 1 - w; } + } - template auto operator()(F const & f) const DECL_AND_RETURN(w1 * f(n) + w2 * f (n+1)); - }; + template auto operator()(F const &f) const DECL_AND_RETURN(w1 *f(n) + w2 *f(n + 1)); + }; - // now evaluator - template struct evaluator : evaluator_one_var{}; + // now evaluator + template struct evaluator : evaluator_one_var {}; } // gfs_implementation. - -}} -#endif +} +} diff --git a/triqs/gfs/local/fourier_matsubara.cpp b/triqs/gfs/local/fourier_matsubara.cpp index e51a0962..eb141322 100644 --- a/triqs/gfs/local/fourier_matsubara.cpp +++ b/triqs/gfs/local/fourier_matsubara.cpp @@ -32,7 +32,7 @@ namespace gfs { struct impl_worker { - tqa::vector g_in, g_out; + arrays::vector g_in, g_out; dcomplex oneFermion(dcomplex a, double b, double tau, double beta) { return -a * (b >= 0 ? exp(-b * tau) / (1 + exp(-beta * b)) : exp(b * (beta - tau)) / (1 + exp(beta * b))); diff --git a/triqs/gfs/local/fourier_matsubara.hpp b/triqs/gfs/local/fourier_matsubara.hpp index bbc20b08..a4c12790 100644 --- a/triqs/gfs/local/fourier_matsubara.hpp +++ b/triqs/gfs/local/fourier_matsubara.hpp @@ -44,8 +44,10 @@ namespace gfs { void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); // The version without tail : only possible in one direction - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); + void triqs_gf_view_assign_delegation(gf_view g, + gf_keeper const& L); + void triqs_gf_view_assign_delegation(gf_view g, + gf_keeper const& L); template gf_mesh make_mesh_fourier_compatible(gf_mesh const& m) { int L = m.size() - (m.kind() == full_bins ? 1 : 0); @@ -59,23 +61,20 @@ namespace gfs { } template - gf_view make_gf_from_fourier(gf_impl const& gt) { + gf make_gf_from_fourier(gf_impl const& gt) { auto gw = gf{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)}; gw() = fourier(gt); return gw; } template - gf_view make_gf_from_inverse_fourier(gf_impl const& gw, - mesh_kind mk = full_bins) { + gf make_gf_from_inverse_fourier(gf_impl const& gw, mesh_kind mk = full_bins) { auto gt = gf{make_mesh_fourier_compatible(gw.mesh(), mk), get_target_shape(gw)}; gt() = inverse_fourier(gw); return gt; } } -} -namespace triqs { namespace clef { TRIQS_CLEF_MAKE_FNT_LAZY(fourier); TRIQS_CLEF_MAKE_FNT_LAZY(inverse_fourier); diff --git a/triqs/gfs/local/fourier_real.cpp b/triqs/gfs/local/fourier_real.cpp index 6ebe06c7..aa600eb0 100644 --- a/triqs/gfs/local/fourier_real.cpp +++ b/triqs/gfs/local/fourier_real.cpp @@ -35,7 +35,7 @@ namespace triqs { namespace gfs { struct impl_worker { - tqa::vector g_in, g_out; + arrays::vector g_in, g_out; void direct (gf_view gw, gf_const_view gt){ @@ -83,7 +83,7 @@ namespace triqs { namespace gfs { const double a = gw.mesh().delta() * sqrt( double(L) ); auto ta = gw(freq_infty()); - tqa::vector g_in(L), g_out(L); + arrays::vector g_in(L), g_out(L); dcomplex t1 = ta(1)(0,0), t2 = ta.get_or_zero(2)(0,0); dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.; diff --git a/triqs/gfs/local/functions.cpp b/triqs/gfs/local/functions.cpp index c1279af7..239e89d2 100644 --- a/triqs/gfs/local/functions.cpp +++ b/triqs/gfs/local/functions.cpp @@ -24,20 +24,20 @@ namespace triqs { namespace gfs { dcomplex F(dcomplex a,double b,double Beta) {return -a/(1+exp(-Beta*b));} - using tqa::array; + using arrays::array; //------------------------------------------------------- // For Imaginary Matsubara Frequency functions // ------------------------------------------------------ - tqa::matrix density( gf_view const & G) { + arrays::matrix density( gf_view const & G) { dcomplex I(0,1); auto sh = G.data().shape().front_pop(); auto Beta = G.domain().beta; local::tail_view t = G(freq_infty()); if (!t.is_decreasing_at_infinity()) TRIQS_RUNTIME_ERROR<<" density computation : Green Function is not as 1/omega or less !!!"; const size_t N1=sh[0], N2 = sh[1]; - tqa::array dens_part(sh), dens_tail(sh), dens(sh); - tqa::matrix res(sh); + arrays::array dens_part(sh), dens_tail(sh), dens(sh); + arrays::matrix res(sh); dens_part()=0;dens()=0;dens_tail()=0; for (size_t n1=0; n1 density( gf_view const & gl) { + arrays::matrix density( gf_view const & gl) { auto sh = gl.data().shape().front_pop(); - tqa::matrix res(sh); + arrays::matrix res(sh); res() = 0.0; for (auto l : gl.mesh()) { @@ -102,21 +102,21 @@ namespace triqs { namespace gfs { } // Impose a discontinuity G(\tau=0)-G(\tau=\beta) - void enforce_discontinuity(gf_view & gl, tqa::array_view disc) { + void enforce_discontinuity(gf_view & gl, arrays::array_view disc) { double norm = 0.0; - tqa::vector t(gl.data().shape()[0]); + arrays::vector t(gl.data().shape()[0]); for (int i=0; i corr(disc.shape()); corr() = 0; + arrays::array corr(disc.shape()); corr() = 0; for (auto l : gl.mesh()) { corr += t(l.index()) * gl[l]; } - tqa::range R; + arrays::range R; for (auto l : gl.mesh()) { gl.data()(l.index(),R,R) += (disc - corr) * t(l.index()) / norm; } diff --git a/triqs/gfs/local/functions.hpp b/triqs/gfs/local/functions.hpp index 284660c2..d336ae45 100644 --- a/triqs/gfs/local/functions.hpp +++ b/triqs/gfs/local/functions.hpp @@ -32,9 +32,9 @@ namespace triqs { // For Imaginary Matsubara Frequency functions // ------------------------------------------------------ - tqa::matrix density(gf_view const & g); + arrays::matrix density(gf_view const & g); - tqa::matrix density(gf_view const & g); + arrays::matrix density(gf_view const & g); local::tail_view get_tail(gf_const_view gl, int size, int omin); @@ -43,7 +43,7 @@ namespace triqs { // For anything that has the ImmutableGfMatsubaraFreq concept, create such a function and compute // Here I choose to create G and call the function to avoid creating one code for each expression... //template - //TYPE_ENABLE_IF (tqa::matrix, ImmutableGfMatsubaraFreq) + //TYPE_ENABLE_IF (arrays::matrix, ImmutableGfMatsubaraFreq) //density( GfType const & G) { return density( gf_view(G));} } diff --git a/triqs/gfs/meshes/discrete.hpp b/triqs/gfs/meshes/discrete.hpp index e71e31c3..4bf9a5a6 100644 --- a/triqs/gfs/meshes/discrete.hpp +++ b/triqs/gfs/meshes/discrete.hpp @@ -18,86 +18,84 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_MESH_DISCRETE_H -#define TRIQS_GF_MESH_DISCRETE_H +#pragma once #include "./mesh_tools.hpp" #include "../domains/discrete.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - template - struct discrete_mesh { + template struct discrete_mesh { - typedef Domain domain_t; - typedef size_t index_t; + using domain_t = Domain; + using index_t = long; - discrete_mesh(domain_t dom) : _dom(std::move(dom)) {} - discrete_mesh() = default; + discrete_mesh(domain_t dom) : _dom(std::move(dom)) {} + discrete_mesh() = default; - domain_t const &domain() const { return _dom; } - size_t size() const {return _dom.size();} + domain_t const &domain() const { return _dom; } + long size() const { return _dom.size(); } - /// Conversions point <-> index <-> discrete_index - size_t index_to_point (index_t ind) const {return ind;} - size_t index_to_linear(index_t ind) const {return ind;} + /// Conversions point <-> index <-> discrete_index + long index_to_point(index_t ind) const { return ind; } + long index_to_linear(index_t ind) const { return ind; } - /// The wrapper for the mesh point - class mesh_point_t : tag::mesh_point,public utility::arithmetic_ops_by_cast { - discrete_mesh const * m; - index_t _index; - public: - mesh_point_t() = default; - mesh_point_t( discrete_mesh const & mesh, index_t const & index_): m(&mesh), _index(index_) {} - void advance() { ++_index;} - typedef size_t cast_t; - operator cast_t() const { return m->index_to_point(_index);} - size_t linear_index() const { return _index;} - size_t index() const { return _index;} - bool at_end() const { return (_index == m->size());} - void reset() {_index =0;} - }; - - /// Accessing a point of the mesh - mesh_point_t operator[](index_t i) const { return mesh_point_t (*this,i);} - mesh_point_t operator[](std::string const & s) const { return mesh_point_t (*this,_dom.index_from_name(s));} - - /// Iterating on all the points... - typedef mesh_pt_generator const_iterator; - const_iterator begin() const { return const_iterator (this);} - const_iterator end() const { return const_iterator (this, true);} - const_iterator cbegin() const { return const_iterator (this);} - const_iterator cend() const { return const_iterator (this, true);} - - /// Mesh comparison - bool operator == (discrete_mesh const & M) const { return (_dom == M._dom) ;} - bool operator != (discrete_mesh const & M) const { return !(operator==(M)); } - - /// Write into HDF5 - friend void h5_write (h5::group fg, std::string subgroup_name, discrete_mesh const & m) { - h5::group gr = fg.create_group(subgroup_name); - h5_write(gr,"domain",m.domain()); - } - - /// Read from HDF5 - friend void h5_read (h5::group fg, std::string subgroup_name, discrete_mesh & m){ - h5::group gr = fg.open_group(subgroup_name); - typename discrete_mesh::domain_t dom; - h5_read(gr,"domain",dom); - m = discrete_mesh(std::move(dom)); - } - - // BOOST Serialization - friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("domain",_dom); - } - - friend std::ostream &operator <<(std::ostream &sout, discrete_mesh const & m){return sout << "Discrete Mesh"; } - - private: - domain_t _dom; + /// The wrapper for the mesh point + class mesh_point_t : tag::mesh_point, public utility::arithmetic_ops_by_cast { + discrete_mesh const *m; + index_t _index; + public: + mesh_point_t() = default; + mesh_point_t(discrete_mesh const &mesh, index_t const &index_) : m(&mesh), _index(index_) {} + void advance() { ++_index; } + using cast_t = long; + operator cast_t() const { return m->index_to_point(_index); } + long linear_index() const { return _index; } + long index() const { return _index; } + bool at_end() const { return (_index == m->size()); } + void reset() { _index = 0; } }; -}} -#endif + + /// Accessing a point of the mesh + mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); } + mesh_point_t operator[](std::string const &s) const { return mesh_point_t(*this, _dom.index_from_name(s)); } + + /// Iterating on all the points... + using const_iterator = mesh_pt_generator; + const_iterator begin() const { return const_iterator(this); } + const_iterator end() const { return const_iterator(this, true); } + const_iterator cbegin() const { return const_iterator(this); } + const_iterator cend() const { return const_iterator(this, true); } + + /// Mesh comparison + bool operator==(discrete_mesh const &M) const { return (_dom == M._dom); } + bool operator!=(discrete_mesh const &M) const { return !(operator==(M)); } + + /// Write into HDF5 + friend void h5_write(h5::group fg, std::string subgroup_name, discrete_mesh const &m) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "domain", m.domain()); + } + + /// Read from HDF5 + friend void h5_read(h5::group fg, std::string subgroup_name, discrete_mesh &m) { + h5::group gr = fg.open_group(subgroup_name); + typename discrete_mesh::domain_t dom; + h5_read(gr, "domain", dom); + m = discrete_mesh(std::move(dom)); + } + + // BOOST Serialization + friend class boost::serialization::access; + template void serialize(Archive &ar, const unsigned int version) { + ar &boost::serialization::make_nvp("domain", _dom); + } + + friend std::ostream &operator<<(std::ostream &sout, discrete_mesh const &m) { return sout << "Discrete Mesh"; } + + private: + domain_t _dom; + }; +} +} diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp index 8aaa90b9..89835b78 100644 --- a/triqs/gfs/meshes/linear.hpp +++ b/triqs/gfs/meshes/linear.hpp @@ -18,8 +18,7 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_MESH_LINEAR_H -#define TRIQS_GF_MESH_LINEAR_H +#pragma once #include "./mesh_tools.hpp" namespace triqs { namespace gfs { @@ -31,15 +30,21 @@ namespace gfs { without_last }; + /** + * Linear mesh + */ template struct linear_mesh { - typedef Domain domain_t; - typedef size_t index_t; - typedef typename domain_t::point_t domain_pt_t; + using domain_t = Domain; + using index_t = long; + using domain_pt_t = typename domain_t::point_t; + + static_assert(!std::is_base_of, domain_pt_t>::value, + "Internal error : can not use Linear Mesh in this case"); linear_mesh() : _dom(), L(0), a_pt(0), b_pt(0), xmin(0), xmax(0), del(0), meshk(half_bins) {} - explicit linear_mesh(domain_t dom, double a, double b, size_t n_pts, mesh_kind mk) + explicit linear_mesh(domain_t dom, double a, double b, long n_pts, mesh_kind mk) : _dom(std::move(dom)), L(n_pts), a_pt(a), b_pt(b), meshk(mk) { switch (mk) { case half_bins: @@ -59,23 +64,17 @@ namespace gfs { } domain_t const &domain() const { return _dom; } - size_t size() const { return L; } + long size() const { return L; } double delta() const { return del; } double x_max() const { return xmax; } double x_min() const { return xmin; } mesh_kind kind() const { return meshk; } /// Conversions point <-> index <-> linear_index - domain_pt_t index_to_point(index_t ind) const { - return embed(xmin + ind * del, std::integral_constant, domain_pt_t>::value>()); - } - - private: // multiply by I is the type is a complex .... - domain_pt_t embed(double x, std::false_type) const { return x; } - domain_pt_t embed(double x, std::true_type) const { return std::complex(0, x); } + domain_pt_t index_to_point(index_t ind) const { return xmin + ind * del; } public: - size_t index_to_linear(index_t ind) const { return ind; } + long index_to_linear(index_t ind) const { return ind; } /// The wrapper for the mesh point class mesh_point_t : tag::mesh_point, public utility::arithmetic_ops_by_cast { @@ -86,24 +85,22 @@ namespace gfs { mesh_point_t() : m(nullptr) {} mesh_point_t(linear_mesh const &mesh, index_t const &index_) : m(&mesh), _index(index_) {} void advance() { ++_index; } - typedef domain_pt_t cast_t; + using cast_t = domain_pt_t; operator cast_t() const { return m->index_to_point(_index); } - size_t linear_index() const { return _index; } - size_t index() const { return _index; } + long linear_index() const { return _index; } + long index() const { return _index; } bool at_end() const { return (_index == m->size()); } void reset() { _index = 0; } }; /// Accessing a point of the mesh - mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); } - - private: - static double real_or_imag(domain_pt_t x, std::false_type) { return x; } - static double real_or_imag(domain_pt_t x, std::true_type) { return imag(x); } + mesh_point_t operator[](index_t i) const { + return {*this, i}; + } public: /// Iterating on all the points... - typedef mesh_pt_generator const_iterator; + using const_iterator = mesh_pt_generator; const_iterator begin() const { return const_iterator(this); } const_iterator end() const { return const_iterator(this, true); } const_iterator cbegin() const { return const_iterator(this); } @@ -142,7 +139,7 @@ namespace gfs { h5::group gr = fg.open_group(subgroup_name); typename linear_mesh::domain_t dom; double a, b; - size_t L; + long L; int k; mesh_kind mk; h5_read(gr, "domain", dom); @@ -181,17 +178,15 @@ namespace gfs { private: domain_t _dom; - size_t L; - double a_pt, b_pt; - double xmin, xmax; - double del; + long L; + double a_pt, b_pt, xmin, xmax, del; mesh_kind meshk; }; - // UNUSED + // UNUSED, only in deprecated code... /// Simple approximation of a point of the domain by a mesh point. No check - template size_t get_closest_mesh_pt_index(linear_mesh const &mesh, typename D::point_t const &x) { + template long get_closest_mesh_pt_index(linear_mesh const &mesh, typename D::point_t const &x) { double a = (x - mesh.x_min()) / mesh.delta(); return std::floor(a); } @@ -212,5 +207,4 @@ namespace gfs { } } } -#endif diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index b7d0f031..3b7e9f5a 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -25,47 +25,75 @@ namespace triqs { namespace gfs { - //----------------------------------- The mesh ----------------------------------------- struct matsubara_freq_mesh { + /// using domain_t = matsubara_domain; - using index_t = int; + /// + using index_t = long; + /// using domain_pt_t = typename domain_t::point_t; + /// Constructor matsubara_freq_mesh() : _dom(), n_max(0), _positive_only(true) {} - explicit matsubara_freq_mesh(domain_t dom, int n_pts, bool positive_only = true) + /// Constructor + matsubara_freq_mesh(domain_t dom, int n_pts=1025, bool positive_only = true) : _dom(std::move(dom)), n_max(n_pts), _positive_only(positive_only) {} + /// Constructor + matsubara_freq_mesh(double beta, statistic_enum S, int Nmax = 1025, bool positive_only = true) + : matsubara_freq_mesh({beta, S}, Nmax, positive_only) {} + + /// Copy constructor + matsubara_freq_mesh(matsubara_freq_mesh const &) = default; + + /// The corresponding domain domain_t const &domain() const { return _dom; } + + /** \brief First value of the index + * + * 0 if positive_only is true + * else : + * For fermions : -Nmax +1 + * For Bosons : -Nmax + **/ int index_start() const { return -(_positive_only ? 0 : n_max + (_dom.statistic == Fermion)); } - size_t size() const { return (_positive_only ? n_max : 2 * n_max + (_dom.statistic == Boson ? 1 : 2)); } - /// Conversions point <-> index <-> linear_index - domain_pt_t index_to_point(index_t ind) const { - return 1_j * M_PI * (2 * ind + (_dom.statistic == Fermion ? 1 : 0)) / _dom.beta; - } + /// Size (linear) of the mesh + long size() const { return (_positive_only ? n_max : 2 * n_max + (_dom.statistic == Boson ? 1 : 2)); } - public: - int index_to_linear(index_t ind) const { return ind + index_start(); } + /// From an index of a point in the mesh, returns the corresponding point in the domain + domain_pt_t index_to_point(index_t ind) const { return 1_j * M_PI * (2 * ind + _dom.statistic == Fermion) / _dom.beta; } - /// The mesh point + /// Flatten the index in the positive linear index for memory storage (almost trivial here). + long index_to_linear(index_t ind) const { return ind + index_start(); } + + /** + * The mesh point + * + * * NB : the mesh point is also in this case a matsubara_freq. + **/ struct mesh_point_t : tag::mesh_point, matsubara_freq { - index_t index_start, index_stop; mesh_point_t() = default; mesh_point_t(matsubara_freq_mesh const &mesh, index_t const &index_) : matsubara_freq(index_, mesh.domain().beta, mesh.domain().statistic), index_start(mesh.index_start()), index_stop(mesh.index_start() + mesh.size() - 1) {} void advance() { ++n; } - size_t linear_index() const { return n; } - size_t index() const { return n; } + long linear_index() const { return n; } + long index() const { return n; } bool at_end() const { return (n == index_stop); } void reset() { n = index_start; } + + private: + index_t index_start, index_stop; }; - /// Accessing a point of the mesh - mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); } + /// Accessing a point of the mesh from its index + mesh_point_t operator[](index_t i) const { + return {*this, i}; + } /// Iterating on all the points... using const_iterator = mesh_pt_generator; @@ -74,7 +102,6 @@ namespace gfs { const_iterator cbegin() const { return const_iterator(this); } const_iterator cend() const { return const_iterator(this, true); } - /// Mesh comparison bool operator==(matsubara_freq_mesh const &M) const { return (std::make_tuple(_dom, n_max, _positive_only, n_max) == std::make_tuple(M._dom, M.n_max, M._positive_only, n_max)); } @@ -86,8 +113,7 @@ namespace gfs { h5_write(gr, "domain", m.domain()); h5_write(gr, "size", m.size()); if (m._positive_only) { - // h5_write(gr, "start_at_0", m._positive_only); - // kept for backward compatibility + // kept ONLY for backward compatibility of archives auto beta = m.domain().beta; h5_write(gr, "min", Fermion ? M_PI / beta : 0); h5_write(gr, "max", ((Fermion ? 1 : 0) + 2 * m.size()) * M_PI / beta); @@ -109,14 +135,15 @@ namespace gfs { m = matsubara_freq_mesh{std::move(dom), L, s}; } - // BOOST Serialization friend class boost::serialization::access; + /// BOOST Serialization template void serialize(Archive &ar, const unsigned int version) { ar &boost::serialization::make_nvp("domain", _dom); ar &boost::serialization::make_nvp("size", n_max); ar &boost::serialization::make_nvp("kind", _positive_only); } + /// Simple print (just blabla and the size) friend std::ostream &operator<<(std::ostream &sout, matsubara_freq_mesh const &m) { return sout << "Matsubara Freq Mesh of size " << m.size(); } @@ -126,6 +153,19 @@ namespace gfs { int n_max; bool _positive_only; }; + + //------------------------------------------------------- + + /** \brief foreach for this mesh + * + * @param m : a mesh + * @param F : a function of synopsis auto F (matsubara_freq_mesh::mesh_point_t) + * + * Calls F on each point of the mesh, in arbitrary order. + **/ + template void foreach(matsubara_freq_mesh const &m, Lambda F) { + for (auto const &w : m) F(w); + } } } diff --git a/triqs/gfs/meshes/matsubara_time.hpp b/triqs/gfs/meshes/matsubara_time.hpp new file mode 100644 index 00000000..7809f027 --- /dev/null +++ b/triqs/gfs/meshes/matsubara_time.hpp @@ -0,0 +1,56 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2013 by O. Parcollet + * + * TRIQS is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * TRIQS. If not, see . + * + ******************************************************************************/ +#pragma once +#include "./mesh_tools.hpp" +#include "../domains/matsubara.hpp" +#include "./linear.hpp" + +namespace triqs { +namespace gfs { + + struct matsubara_time_mesh : linear_mesh { + private: + using B = linear_mesh; + + public: + matsubara_time_mesh() = default; + matsubara_time_mesh(matsubara_time_mesh const &x) = default; + matsubara_time_mesh(matsubara_time_domain d, int n_time_slices, mesh_kind mk = half_bins) : B(d, 0, d.beta, n_time_slices, mk) {} + matsubara_time_mesh(double beta, statistic_enum S, int n_time_slices, mesh_kind mk = half_bins) + : matsubara_time_mesh({beta, S}, n_time_slices, mk) {} + }; + + + //------------------------------------------------------- + + /** \brief foreach for this mesh + * + * @param m : a mesh + * @param F : a function of synopsis auto F (matsubara_time_mesh::mesh_point_t) + * + * Calls F on each point of the mesh, in arbitrary order. + **/ + template void foreach(matsubara_time_mesh const &m, Lambda F) { + for (auto const &w : m) F(w); + } +} +} + diff --git a/triqs/gfs/meshes/product.hpp b/triqs/gfs/meshes/product.hpp index 12b4107e..ac8f4a42 100644 --- a/triqs/gfs/meshes/product.hpp +++ b/triqs/gfs/meshes/product.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,170 +18,254 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_MESH_PRODUCT_H -#define TRIQS_GF_MESH_PRODUCT_H +#pragma once #include "./mesh_tools.hpp" #include "../domains/product.hpp" #include #include #include -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - template struct mesh_product : tag::composite { - typedef domain_product domain_t; - typedef std::c14::tuple index_t; - typedef std::tuple m_tuple_t; - typedef std::tuple m_pt_tuple_t; - typedef typename domain_t::point_t domain_pt_t; + /** Cartesian product of meshes + */ + template struct mesh_product : tag::composite { + using domain_t = domain_product; + using index_t = std::c14::tuple; + using m_tuple_t = std::tuple; + using m_pt_tuple_t = std::tuple; + using domain_pt_t = typename domain_t::point_t; - static constexpr int dim = sizeof...(Meshes); + static constexpr int dim = sizeof...(Meshes); - mesh_product () {} - mesh_product (Meshes const & ... meshes) : m_tuple(meshes...), _dom(meshes.domain()...) {} + mesh_product() {} + mesh_product(Meshes const &... meshes) : m_tuple(meshes...), _dom(meshes.domain()...) {} - domain_t const & domain() const { return _dom;} - m_tuple_t const & components() const { return m_tuple;} - m_tuple_t & components() { return m_tuple;} - - /// size of the mesh is the product of size - struct _aux0 { template size_t operator()(M const & m, size_t R) { return R*m.size();}}; - size_t size() const { return triqs::tuple::fold(_aux0(), m_tuple, 1);} - - /// Conversions point <-> index <-> linear_index - struct _aux1 { template void operator()(P & p, M const & m, I const& i) {p = m.index_to_point(i);}}; - typename domain_t::point_t index_to_point(index_t const & ind) const { domain_pt_t res; triqs::tuple::apply_on_zip(_aux1(), res,m_tuple,ind); return res;} - - // index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....)) - struct _aux2 { template size_t operator()(M const & m, I const & i,size_t R) {return m.index_to_linear(i) + R * m.size();}}; - size_t index_to_linear(index_t const & ii) const { return triqs::tuple::fold_on_zip(_aux2(), reverse(m_tuple), reverse(ii), size_t(0)); } -// size_t index_to_linear(index_t const & ii) const { return triqs::tuple::fold_on_zip([](auto const &m, auto const &i, auto R){return m.index_to_linear(i) + R * m.size();} , m_tuple, ii, size_t(0)); } - - - // Same but a tuple of mesh_point_t - struct _aux3 { template size_t operator()(M const & m, P const & p,size_t R) {return p.linear_index() + R * m.size();}}; - size_t mp_to_linear(m_pt_tuple_t const & mp) const { return triqs::tuple::fold_on_zip(_aux3(), reverse(m_tuple), reverse(mp), size_t(0)); } - - // - struct _aux4 { template< typename M, typename V> V * operator()(M const & m, V * v) {*v = m.size(); return ++v;}}; - utility::mini_vector all_size_as_mini_vector () const { - utility::mini_vector res; - triqs::tuple::fold(_aux4(), m_tuple, &res[0] ); - return res; - } - - // Same but a variadic list of mesh_point_t - template size_t mesh_pt_components_to_linear(MP const & ... mp) const { - static_assert(std::is_same< std::tuple, m_pt_tuple_t>::value, "Call incorrect "); - //static_assert(std::is_same< std::tuple::type>::type...>, m_pt_tuple_t>::value, "Call incorrect "); - return mp_to_linear(std::forward_as_tuple(mp...)); - } // speed test ? or make a variadic fold... - - /// The wrapper for the mesh point - class mesh_point_t : tag::mesh_point{ - const mesh_product * m; - m_pt_tuple_t _c; bool _atend; - struct F2 { template typename M::mesh_point_t operator()(M const & m, typename M::index_t const & i) const { return m[i];}}; - struct F1 { template typename M::mesh_point_t operator()(M const & m) const { return m[typename M::index_t()];}}; - public : - mesh_point_t() = default; - mesh_point_t(mesh_product const & m_, index_t index_ ) : m(&m_), _c (triqs::tuple::apply_on_zip(F2(), m_.m_tuple, index_)), _atend(false) {} - mesh_point_t(mesh_product const & m_) : m(&m_), _c (triqs::tuple::apply(F1(), m_.m_tuple)), _atend(false) {} - m_pt_tuple_t const & components_tuple() const { return _c;} - size_t linear_index() const { return m->mp_to_linear(_c);} - const mesh_product * mesh() const { return m;} - - typedef domain_pt_t cast_t; - operator cast_t() const { return m->index_to_point(index);} - - // index[0] +=1; if index[0]==m.component[0].size() { index[0]=0; index[1] +=1; if ....} and so on until dim - struct _aux1 { template bool operator()(P & p, bool done) - {if (done) return true; p.advance(); if (p.at_end()) {p.reset(); return false;} return true;} - }; - void advance() { triqs::tuple::fold(_aux1(), _c, false);} - - //index_t index() const { return _index;} // not implemented yet - bool at_end() const { return _atend;} - - struct _aux{ template size_t operator()(M & m,size_t ) { m.reset(); return 0;}}; - void reset() { _atend = false; triqs::tuple::fold(_aux(), _c,0);} - };// end mesh_point_t - - /// Accessing a point of the mesh - mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i);} - mesh_point_t operator()(typename Meshes::index_t ... i) const { return (*this)[std::make_tuple(i...)];} - - /// Iterating on all the points... - typedef mesh_pt_generator const_iterator; - const_iterator begin() const { return const_iterator (this);} - const_iterator end() const { return const_iterator (this, true);} - const_iterator cbegin() const { return const_iterator (this);} - const_iterator cend() const { return const_iterator (this, true);} - - /// Mesh comparison - friend bool operator == (mesh_product const & M1, mesh_product const & M2) { return M1.m_tuple==M2.m_tuple; } - - /// Write into HDF5 - struct _auxh5w { - h5::group gr; _auxh5w( h5::group gr_) : gr(gr_) {} //icc has yet another bug on new initialization form with {}... - template size_t operator()(M const & m, size_t N) { std::stringstream fs;fs <<"MeshComponent"<< N; h5_write(gr,fs.str(), m); return N+1; } - }; - friend void h5_write (h5::group fg, std::string subgroup_name, mesh_product const & m) { - h5::group gr = fg.create_group(subgroup_name); - //h5_write(gr,"domain",m.domain()); - triqs::tuple::fold(_auxh5w(gr), m.components(), size_t(0)); - } - - /// Read from HDF5 - struct _auxh5r { - h5::group gr;_auxh5r( h5::group gr_) : gr(gr_) {} - template size_t operator()(M & m, size_t N) { std::stringstream fs;fs <<"MeshComponent"<< N; h5_read(gr,fs.str(), m); return N+1; } - }; - friend void h5_read (h5::group fg, std::string subgroup_name, mesh_product & m){ - h5::group gr = fg.open_group(subgroup_name); - //h5_read(gr,"domain",m._dom); - triqs::tuple::fold(_auxh5r(gr), m.components(), size_t(0)); - } - - // BOOST Serialization - friend class boost::serialization::access; - template struct _aux_ser { - Archive & ar;_aux_ser( Archive & ar_) : ar(ar_) {} - template size_t operator()(M & m, size_t N) { - std::stringstream fs;fs <<"MeshComponent"<< N; - ar & boost::serialization::make_nvp(fs.str().c_str(),m); - return N+1; - } - }; - template - void serialize(Archive & ar, const unsigned int version) { - triqs::tuple::fold(_aux_ser(ar), m_tuple, size_t(0)); - } - - friend std::ostream &operator <<(std::ostream &sout, mesh_product const & m){return sout << "Product Mesh"; } + domain_t const &domain() const { return _dom; } + m_tuple_t const &components() const { return m_tuple; } + m_tuple_t &components() { return m_tuple; } private: - m_tuple_t m_tuple; - domain_t _dom; + struct _aux0 { + template size_t operator()(M const &m, size_t R) { return R * m.size(); } + }; + + public: + /// size of the mesh is the product of size + size_t size() const { return triqs::tuple::fold(_aux0(), m_tuple, 1); } + + private: + struct _aux1 { + template void operator()(P &p, M const &m, I const &i) { p = m.index_to_point(i); } + }; + + public: + /// Conversions point <-> index <-> linear_index + typename domain_t::point_t index_to_point(index_t const &ind) const { + domain_pt_t res; + triqs::tuple::apply_on_zip(_aux1(), res, m_tuple, ind); + return res; + } + + private: + struct _aux2 { + template size_t operator()(M const &m, I const &i, size_t R) { + return m.index_to_linear(i) + R * m.size(); + } + }; + + public: + /// Flattening index to linear : index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....)) + size_t index_to_linear(index_t const &ii) const { + return triqs::tuple::fold_on_zip(_aux2(), reverse(m_tuple), reverse(ii), size_t(0)); + } + // size_t index_to_linear(index_t const & ii) const { return triqs::tuple::fold_on_zip([](auto const &m, auto const &i, auto R) + //{return m.index_to_linear(i) + R * m.size();} , m_tuple, ii, size_t(0)); } + + private: + struct _aux3 { + template size_t operator()(M const &m, P const &p, size_t R) { + return p.linear_index() + R * m.size(); + } + }; + + public: + /// Flattening index to linear : index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....)) + size_t mp_to_linear(m_pt_tuple_t const &mp) const { + return triqs::tuple::fold_on_zip(_aux3(), reverse(m_tuple), reverse(mp), size_t(0)); + } + + // + private: + struct _aux4 { + template V *operator()(M const &m, V *v) { + *v = m.size(); + return ++v; + } + }; + + public: + utility::mini_vector all_size_as_mini_vector() const { + utility::mini_vector res; + triqs::tuple::fold(_aux4(), m_tuple, &res[0]); + return res; + } + + // Same but a variadic list of mesh_point_t + template size_t mesh_pt_components_to_linear(MP const &... mp) const { + static_assert(std::is_same, m_pt_tuple_t>::value, "Call incorrect "); + // static_assert(std::is_same< std::tuple::type>::type...>, + // m_pt_tuple_t>::value, "Call incorrect "); + return mp_to_linear(std::forward_as_tuple(mp...)); + } // speed test ? or make a variadic fold... + + /// The wrapper for the mesh point + class mesh_point_t : tag::mesh_point { + const mesh_product *m; + m_pt_tuple_t _c; + bool _atend; + struct F2 { + template typename M::mesh_point_t operator()(M const &m, typename M::index_t const &i) const { return m[i]; } + }; + struct F1 { + template typename M::mesh_point_t operator()(M const &m) const { return m[typename M::index_t()]; } + }; + + public: + mesh_point_t() = default; + mesh_point_t(mesh_product const &m_, index_t index_) + : m(&m_), _c(triqs::tuple::apply_on_zip(F2(), m_.m_tuple, index_)), _atend(false) {} + mesh_point_t(mesh_product const &m_) : m(&m_), _c(triqs::tuple::apply(F1(), m_.m_tuple)), _atend(false) {} + m_pt_tuple_t const &components_tuple() const { return _c; } + size_t linear_index() const { return m->mp_to_linear(_c); } + const mesh_product *mesh() const { return m; } + + using cast_t = domain_pt_t; + operator cast_t() const { return m->index_to_point(index); } + + // index[0] +=1; if index[0]==m.component[0].size() { index[0]=0; index[1] +=1; if ....} and so on until dim + private: + struct _aux1 { + template bool operator()(P &p, bool done) { + if (done) return true; + p.advance(); + if (!p.at_end()) return true; + p.reset(); + return false; + } + }; + + public: + void advance() { triqs::tuple::fold(_aux1(), _c, false); } + + // index_t index() const { return _index;} // not implemented yet + bool at_end() const { return _atend; } + + private: + struct _aux { + template size_t operator()(M &m, size_t) { + m.reset(); + return 0; + } + }; + + public: + void reset() { + _atend = false; + triqs::tuple::fold(_aux(), _c, 0); + } + }; // end mesh_point_t + + /// Accessing a point of the mesh + mesh_point_t operator[](index_t i) const { return mesh_point_t(*this, i); } + mesh_point_t operator()(typename Meshes::index_t... i) const { return (*this)[std::make_tuple(i...)]; } + + /// Iterating on all the points... + using const_iterator = mesh_pt_generator; + const_iterator begin() const { return const_iterator(this); } + const_iterator end() const { return const_iterator(this, true); } + const_iterator cbegin() const { return const_iterator(this); } + const_iterator cend() const { return const_iterator(this, true); } + + /// Mesh comparison + friend bool operator==(mesh_product const &M1, mesh_product const &M2) { return M1.m_tuple == M2.m_tuple; } + + private: + + /// Write into HDF5 + struct _auxh5w { + h5::group gr; + _auxh5w(h5::group gr_) : gr(gr_) {} // icc has yet another bug on new initialization form with {}... + template size_t operator()(M const &m, size_t N) { + std::stringstream fs; + fs << "MeshComponent" << N; + h5_write(gr, fs.str(), m); + return N + 1; + } + }; + + friend void h5_write(h5::group fg, std::string subgroup_name, mesh_product const &m) { + h5::group gr = fg.create_group(subgroup_name); + // h5_write(gr,"domain",m.domain()); + triqs::tuple::fold(_auxh5w(gr), m.components(), size_t(0)); + } + + /// Read from HDF5 + struct _auxh5r { + h5::group gr; + _auxh5r(h5::group gr_) : gr(gr_) {} + template size_t operator()(M &m, size_t N) { + std::stringstream fs; + fs << "MeshComponent" << N; + h5_read(gr, fs.str(), m); + return N + 1; + } + }; + friend void h5_read(h5::group fg, std::string subgroup_name, mesh_product &m) { + h5::group gr = fg.open_group(subgroup_name); + // h5_read(gr,"domain",m._dom); + triqs::tuple::fold(_auxh5r(gr), m.components(), size_t(0)); + } + + // BOOST Serialization + friend class boost::serialization::access; + template struct _aux_ser { + Archive &ar; + _aux_ser(Archive &ar_) : ar(ar_) {} + template size_t operator()(M &m, size_t N) { + std::stringstream fs; + fs << "MeshComponent" << N; + ar &boost::serialization::make_nvp(fs.str().c_str(), m); + return N + 1; + } + }; + template void serialize(Archive &ar, const unsigned int version) { + triqs::tuple::fold(_aux_ser(ar), m_tuple, size_t(0)); + } + + friend std::ostream &operator<<(std::ostream &sout, mesh_product const &m) { return sout << "Product Mesh"; } + + private: + m_tuple_t m_tuple; + domain_t _dom; }; - template - auto get_index(P const & p) DECL_AND_RETURN( std::get(p.components_tuple()).index()); + template auto get_index(P const &p) DECL_AND_RETURN(std::get(p.components_tuple()).index()); + + template + auto get_point(P const &p) + DECL_AND_RETURN(std::get(p.mesh() -> components()).index_to_point(std::get(p.components_tuple()).index())); + + template auto get_component(P const &p) DECL_AND_RETURN(std::get(p.components_tuple())); - template - auto get_point(P const & p) DECL_AND_RETURN( std::get( p.mesh()->components() ).index_to_point( std::get(p.components_tuple()).index() ) ); - - template - auto get_component(P const & p) DECL_AND_RETURN( std::get(p.components_tuple())); - // Given a composite mesh m , and a linear array of storage A // reinterpret_linear_array(m,A) returns a d-dimensionnal view of the array // with indices egal to the indices of the components of the mesh. // Very useful for slicing, currying functions. - template - arrays::array_view - reinterpret_linear_array(mesh_product const & m, arrays::array_view A) { - return { {join (m.all_size_as_mini_vector(), get_shape(A).front_pop())}, A.storage()}; - } - -}} -#endif + template + arrays::array_view + reinterpret_linear_array(mesh_product const &m, arrays::array_view A) { + return {{join(m.all_size_as_mini_vector(), get_shape(A).front_pop())}, A.storage()}; + } +} +} diff --git a/triqs/gfs/meshes/segment.hpp b/triqs/gfs/meshes/segment.hpp new file mode 100644 index 00000000..180bb145 --- /dev/null +++ b/triqs/gfs/meshes/segment.hpp @@ -0,0 +1,49 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2013 by O. Parcollet + * + * TRIQS is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * TRIQS. If not, see . + * + ******************************************************************************/ +#pragma once +#include "./mesh_tools.hpp" +#include "./linear.hpp" + +namespace triqs { +namespace gfs { + + /** A linear mesh on a segment on R */ + struct segment_mesh : linear_mesh { + using B = linear_mesh; + segment_mesh() = default; + segment_mesh(double x_min, double x_max, int n_freq, mesh_kind mk = full_bins) : B(typename B::domain_t(), x_min, x_max, n_freq, mk) {} + }; + + //------------------------------------------------------- + + /** \brief foreach for this mesh + * + * @param m : a mesh + * @param F : a function of synopsis auto F (matsubara_time_mesh::mesh_point_t) + * + * Calls F on each point of the mesh, in arbitrary order. + **/ + template void foreach(segment_mesh const &m, Lambda F) { + for (auto const &w : m) F(w); + } +} +} + diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp index 21f72f4d..0f4f11cd 100644 --- a/triqs/gfs/product.hpp +++ b/triqs/gfs/product.hpp @@ -1,5 +1,5 @@ /******************************************************************************* - * + * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2013 by O. Parcollet @@ -18,86 +18,86 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_PRODUCT_H -#define TRIQS_GF_PRODUCT_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" #include "./meshes/product.hpp" #include "./evaluators.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - template struct cartesian_product{ - typedef std::tuple type; - static constexpr size_t size = sizeof...(Ms); + template struct cartesian_product { + using type = std::tuple; + static constexpr size_t size = sizeof...(Ms); }; - // use alias - template struct cartesian_product > : cartesian_product{}; + // use alias + template struct cartesian_product> : cartesian_product {}; // the mesh is simply a cartesian product - template struct gf_mesh,Opt> : mesh_product< gf_mesh ... > { - typedef mesh_product< gf_mesh ... > B; - typedef std::tuple mesh_name_t; + template struct gf_mesh, Opt> : mesh_product...> { + // using mesh_product< gf_mesh ... >::mesh_product< gf_mesh ... > ; + using B = mesh_product...>; gf_mesh() = default; - gf_mesh (gf_mesh ... ms) : B {std::move(ms)...} {} + gf_mesh(gf_mesh... ms) : B{std::move(ms)...} {} }; - namespace gfs_implementation { + namespace gfs_implementation { /// --------------------------- hdf5 --------------------------------- // h5 name : name1_x_name2_..... - template struct h5_name,matrix_valued,Opt> { - static std::string invoke(){ - return triqs::tuple::fold( - [](std::string a, std::string b) { return a + std::string(b.empty()?"" : "_x_") + b;}, - std::make_tuple(h5_name::invoke()...), - std::string()); + template struct h5_name, matrix_valued, Opt> { + static std::string invoke() { + return triqs::tuple::fold([](std::string a, std::string b) { return a + std::string(b.empty() ? "" : "_x_") + b; }, + std::make_tuple(h5_name::invoke()...), std::string()); } }; - template struct h5_name,tensor_valued,Opt> : h5_name,matrix_valued,Opt> {}; + template + struct h5_name, tensor_valued, Opt> : h5_name, matrix_valued, Opt> {}; // a slight difference with the generic case : reinterpret the data array to avoid flattening the variables - template - struct h5_rw,tensor_valued,Opt> { - typedef gf,tensor_valued,Opt> g_t; + template struct h5_rw, tensor_valued, Opt> { + using g_t = gf, tensor_valued, Opt>; - static void write (h5::group gr, typename g_t::const_view_type g) { - h5_write(gr,"data",reinterpret_linear_array(g.mesh(), g().data())); - h5_write(gr,"singularity",g._singularity); - h5_write(gr,"mesh",g._mesh); - h5_write(gr,"symmetry",g._symmetry); - } + static void write(h5::group gr, typename g_t::const_view_type g) { + h5_write(gr, "data", reinterpret_linear_array(g.mesh(), g().data())); + h5_write(gr, "singularity", g._singularity); + h5_write(gr, "mesh", g._mesh); + h5_write(gr, "symmetry", g._symmetry); + } - template - static void read (h5::group gr, gf_impl,tensor_valued,Opt,IsView,false> & g) { - using G_t= gf_impl,tensor_valued,Opt,IsView,false> ; - h5_read(gr,"mesh",g._mesh); - auto arr = arrays::array{}; - h5_read(gr,"data",arr); - auto sh = arr.shape(); - arrays::mini_vector sh2; - sh2[0] = g._mesh.size(); - for (int u=1; u{sh2, std::move(arr.storage())}; - h5_read(gr,"singularity",g._singularity); - h5_read(gr,"symmetry",g._symmetry); - } - - }; + template + static void read(h5::group gr, gf_impl, tensor_valued, Opt, IsView, false> &g) { + using G_t = gf_impl, tensor_valued, Opt, IsView, false>; + h5_read(gr, "mesh", g._mesh); + auto arr = arrays::array{}; + h5_read(gr, "data", arr); + auto sh = arr.shape(); + arrays::mini_vector sh2; + sh2[0] = g._mesh.size(); + for (int u = 1; u < R + 1; ++u) sh2[u] = sh[sizeof...(Ms) - 1 + u]; + g._data = arrays::array{sh2, std::move(arr.storage())}; + h5_read(gr, "singularity", g._singularity); + h5_read(gr, "symmetry", g._symmetry); + } + }; /// --------------------------- data access --------------------------------- - template struct data_proxy,scalar_valued,Opt> : data_proxy_array,1> {}; - template struct data_proxy,matrix_valued,Opt> : data_proxy_array,3> {}; - template struct data_proxy,tensor_valued,Opt> : data_proxy_array,R+1> {}; + template + struct data_proxy, scalar_valued, Opt> : data_proxy_array, 1> {}; + template + struct data_proxy, matrix_valued, Opt> : data_proxy_array, 3> {}; + template + struct data_proxy, tensor_valued, Opt> : data_proxy_array, R + 1> {}; /// --------------------------- evaluator --------------------------------- - /** + /** * This the multi-dimensional evaluator. - * It combine the evaluator of each components, as long as they are a linear form + * It combine the evaluator of each components, as long as they are a linear form * eval(g, x) = \sum_i w_i g( n_i(x)) , with w some weight and n_i some points on the grid. * Mathematically, it is written as (example of evaluating g(x1,x2,x3,x4)). * Notation : eval(X) : g -> g(X) @@ -110,54 +110,63 @@ namespace triqs { namespace gfs { * p_i are points on the grids, x_i points in the domain. * * Unrolling the formula gives (for 2 variables, with 2 points interpolation) - * eval(xa,xb) (g) = eval (xa) ( binder ( g, (), (xb)) ) = w_1(xa) binder ( g, (), (xb))( n_1(xa)) + w_2(xa) binder ( g, (), (xb))( n_2(xa)) - * = w_1(xa) ( eval(xb)( binder ( g, (n_1(xa) ), ()))) + 1 <-> 2 - * = w_1(xa) ( W_1(xb) * binder ( g, (n_1(xa) ), ())(N_1(xb)) + 1<->2 ) + 1 <-> 2 - * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + 1<->2 ) + 1 <-> 2 - * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + W_2(xb) * g[n_1(xa), N_2(xb)] ) + 1 <-> 2 + * eval(xa,xb) (g) = eval (xa) ( binder ( g, (), (xb)) ) = + * w_1(xa) binder ( g, (), (xb))( n_1(xa)) + w_2(xa) binder ( g, (), (xb))( n_2(xa)) + * = w_1(xa) ( eval(xb)( binder ( g, (n_1(xa) ), ()))) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * binder ( g, (n_1(xa) ), ())(N_1(xb)) + 1<->2 ) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + 1<->2 ) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + W_2(xb) * g[n_1(xa), N_2(xb)] ) + 1 <-> 2 * which is the expected formula */ - // implementation : G = gf, Tn : tuple of n points, Ev : tuple of evaluators (the evals functions), pos = counter from #args-1 =>0 + // implementation : G = gf, Tn : tuple of n points, Ev : tuple of evaluators (the evals functions), + // pos = counter from #args-1 =>0 // NB : the tuple is build in reverse with respect to the previous comment. - template struct binder; + template struct binder; - template - binder make_binder(G const * g, Tn tn, Ev const & ev) { return binder{g, std::move(tn), ev}; } + template binder make_binder(G const *g, Tn tn, Ev const &ev) { + return binder{g, std::move(tn), ev}; + } - template struct binder { - G const * g; Tn tn; Ev const & evals; - auto operator()(size_t p) const DECL_AND_RETURN( std::get(evals) ( make_binder(g, triqs::tuple::push_front(tn,p), evals) )); + template struct binder { + G const *g; + Tn tn; + Ev const &evals; + auto operator()(size_t p) const + DECL_AND_RETURN(std::get(evals)(make_binder(g, triqs::tuple::push_front(tn, p), evals))); }; - template struct binder { - G const * g; Tn tn; Ev const & evals; - auto operator()(size_t p) const DECL_AND_RETURN( triqs::tuple::apply(on_mesh(*g), triqs::tuple::push_front(tn,p))); + template struct binder { + G const *g; + Tn tn; + Ev const &evals; + auto operator()(size_t p) const DECL_AND_RETURN(triqs::tuple::apply(on_mesh(*g), triqs::tuple::push_front(tn, p))); }; // now the multi d evaluator itself. - template - struct evaluator,Target,Opt> { - static constexpr int arity = sizeof...(Ms); - mutable std::tuple< evaluator_fnt_on_mesh ... > evals; + template struct evaluator, Target, Opt> { + static constexpr int arity = sizeof...(Ms); + mutable std::tuple...> evals; - struct _poly_lambda {// replace by a polymorphic lambda in C++14 - template void operator()(A & a, B const & b, C const & c) const { a = A{b,c};} - }; - - template - //std::complex operator() (G const * g, Args && ... args) const { - auto operator() (G const * g, Args && ... args) const - -> decltype (std::get(evals) (make_binder (g, std::make_tuple(), evals) )) - // when do we get C++14 decltype(auto) ...!? - { - static constexpr int R = sizeof...(Args); - // build the evaluators, as a tuple of ( evaluator ( mesh_component, args)) - triqs::tuple::call_on_zip(_poly_lambda(), evals, g->mesh().components(), std::make_tuple(args...)); - return std::get(evals) (make_binder (g, std::make_tuple(), evals) ); - } + struct _poly_lambda { // replace by a polymorphic lambda in C++14 + template void operator()(A &a, B const &b, C const &c) const { + a = A{b, c}; + } }; - } // gf_implementation - }} -#endif + template + // std::complex operator() (G const * g, Args && ... args) const { + auto operator()(G const *g, Args &&... args) + const -> decltype(std::get(evals)(make_binder(g, std::make_tuple(), evals))) + // when do we get C++14 decltype(auto) ...!? + { + static constexpr int R = sizeof...(Args); + // build the evaluators, as a tuple of ( evaluator ( mesh_component, args)) + triqs::tuple::call_on_zip(_poly_lambda(), evals, g->mesh().components(), std::make_tuple(args...)); + return std::get(evals)(make_binder(g, std::make_tuple(), evals)); + } + }; + + } // gf_implementation +} +} diff --git a/triqs/gfs/refreq.hpp b/triqs/gfs/refreq.hpp index 0b1f9bbb..43c66387 100644 --- a/triqs/gfs/refreq.hpp +++ b/triqs/gfs/refreq.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,47 +18,52 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_FREQ_H -#define TRIQS_GF_FREQ_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" #include "./domains/R.hpp" -#include "./meshes/linear.hpp" +#include "./meshes/segment.hpp" #include "./evaluators.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { struct refreq {}; - template struct gf_mesh : linear_mesh { - typedef linear_mesh B; - gf_mesh() = default; - gf_mesh (double wmin, double wmax, int n_freq, mesh_kind mk=full_bins) : - B(typename B::domain_t(), wmin, wmax, n_freq, mk){} + template struct gf_mesh : segment_mesh { + template gf_mesh(T &&... x) : segment_mesh(std::forward(x)...) {} + //using segment_mesh::segment_mesh; }; - namespace gfs_implementation { + namespace gfs_implementation { - // singularity - template struct singularity { typedef local::tail type;}; - template struct singularity { typedef local::tail type;}; + // singularity + template struct singularity { + using type = local::tail; + }; + template struct singularity { + using type = local::tail; + }; // h5 name - template struct h5_name { static std::string invoke(){ return "ReFreq";}}; + template struct h5_name { + static std::string invoke() { return "ReFreq"; } + }; /// --------------------------- evaluator --------------------------------- - template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation); + template <> + struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, + evaluator_grid_linear_interpolation); - template struct evaluator : evaluator_one_var{}; + template struct evaluator : evaluator_one_var {}; /// --------------------------- data access --------------------------------- - template struct data_proxy : data_proxy_array,3> {}; - template struct data_proxy : data_proxy_array,1> {}; - + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; } -}} -#endif +} +} diff --git a/triqs/gfs/retime.hpp b/triqs/gfs/retime.hpp index 612f5f7a..e15d690e 100644 --- a/triqs/gfs/retime.hpp +++ b/triqs/gfs/retime.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2013 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,44 +18,50 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_ONE_REAL_TIME_H -#define TRIQS_GF_ONE_REAL_TIME_H +#pragma once #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" #include "./domains/R.hpp" -#include "./meshes/linear.hpp" +#include "./meshes/segment.hpp" #include "./evaluators.hpp" -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { struct retime {}; - template struct gf_mesh : linear_mesh { - typedef linear_mesh B; - gf_mesh() = default; - gf_mesh(double tmin, double tmax, int n_points, mesh_kind mk=full_bins) : B (typename B::domain_t(), tmin, tmax, n_points, mk){} + template struct gf_mesh : segment_mesh { + template gf_mesh(T &&... x) : segment_mesh(std::forward(x)...) {} + //using segment_mesh::segment_mesh; }; - namespace gfs_implementation { + namespace gfs_implementation { - // singularity - template struct singularity { typedef local::tail type;}; - template struct singularity { typedef local::tail type;}; + // singularity + template struct singularity { + using type = local::tail; + }; + template struct singularity { + using type = local::tail; + }; // h5 name - template struct h5_name { static std::string invoke(){ return "ReTime";}}; + template struct h5_name { + static std::string invoke() { return "ReTime"; } + }; /// --------------------------- evaluator --------------------------------- - template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation); + template <> + struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, + evaluator_grid_linear_interpolation); - template struct evaluator : evaluator_one_var{}; + template struct evaluator : evaluator_one_var {}; /// --------------------------- data access --------------------------------- - template struct data_proxy : data_proxy_array,3> {}; - template struct data_proxy : data_proxy_array,1> {}; + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; } // gfs_implementation -}} -#endif - +} +} diff --git a/triqs/gfs/tools.hpp b/triqs/gfs/tools.hpp index 8f3fe431..3b7294ac 100644 --- a/triqs/gfs/tools.hpp +++ b/triqs/gfs/tools.hpp @@ -18,90 +18,89 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_TOOLS_H -#define TRIQS_GF_TOOLS_H -#include -#include -#include -#include -#include -#include -#include +#pragma once +#include #include "triqs/utility/complex_ops.hpp" #include #include -#include +#include +#include -namespace triqs { namespace gfs { - namespace tqa= triqs::arrays; - namespace mpl=boost::mpl; +namespace triqs { +namespace gfs { + namespace mpl = boost::mpl; - namespace tag { struct composite{}; struct mesh_point{};} - - struct scalar_valued {}; - template struct tensor_valued {static_assert( R>0, "tensor_valued only for rank >0");}; - struct matrix_valued{}; - - //------------------------------------------------------ - - typedef std::complex dcomplex; - - enum statistic_enum {Boson,Fermion}; - - struct freq_infty{}; // the point at infinity - - //------------------------------------------------------ - - inline std::vector split(const std::string &s, char delim){ - std::vector elems; - std::stringstream ss(s); - std::string item; - while(std::getline(ss, item, delim)) { elems.push_back(item); } - return elems; + namespace tag { + struct composite {}; + struct mesh_point {}; } - //------------------------------------------------------ - template struct closest_pt_wrap; - - template struct closest_pt_wrap : tag::mesh_point { T value; template explicit closest_pt_wrap(U&&x): value(std::forward(x)){} }; + // scalar_valued, matrix_valued, tensor_valued + struct scalar_valued {}; - template struct closest_pt_wrap : tag::mesh_point { - std::tuple value_tuple; - template explicit closest_pt_wrap(U&& ... x): value_tuple (std::forward(x) ... ){} + template struct tensor_valued { + static_assert(R > 0, "tensor_valued only for rank >0"); }; - - template closest_pt_wrap closest_mesh_pt(T && ... x) { return closest_pt_wrap (std::forward(x)...);} + + struct matrix_valued {}; //------------------------------------------------------ + using dcomplex = std::complex; + + /** The statistics : Boson or Fermion + */ + enum statistic_enum { + Boson, + Fermion + }; + + struct freq_infty {}; // the point at infinity + + //------------------------------------------------------ + + template struct closest_pt_wrap; + + template struct closest_pt_wrap : tag::mesh_point { + T value; + template explicit closest_pt_wrap(U &&x) : value(std::forward(x)) {} + }; + + template struct closest_pt_wrap : tag::mesh_point { + std::tuple value_tuple; + template explicit closest_pt_wrap(U &&... x) : value_tuple(std::forward(x)...) {} + }; + + template closest_pt_wrap closest_mesh_pt(T &&... x) { + return closest_pt_wrap{std::forward(x)...}; + } + + //------------------------------------------------------ + + // A simple replacement of tail when there is none to maintain generic code simple... struct nothing { - template explicit nothing(Args&&...) {} // takes anything, do nothing.. + template explicit nothing(Args &&...) {} // takes anything, do nothing.. nothing() {} - typedef nothing view_type; - typedef nothing regular_type; - void rebind (nothing){} - template< typename RHS> void operator=(RHS && ) {} - friend void h5_write (h5::group, std::string subgroup_name, nothing ) {} - friend void h5_read (h5::group, std::string subgroup_name, nothing ) {} - // BOOST Serialization + using view_type = nothing; + using regular_type = nothing; + void rebind(nothing) {} + template void operator=(RHS &&) {} + friend void h5_write(h5::group, std::string subgroup_name, nothing) {} + friend void h5_read(h5::group, std::string subgroup_name, nothing) {} friend class boost::serialization::access; - template - void serialize(Archive & ar, const unsigned int version) { - } - friend nothing operator +( nothing, nothing) { return nothing();} - template friend void assign_from_expression(nothing & ,RHS) {} + template void serialize(Archive &ar, const unsigned int version) {} + friend nothing operator+(nothing, nothing) { return nothing(); } + template friend void assign_from_expression(nothing &, RHS) {} + }; -}; - -template nothing slice_target(nothing, T...) { return nothing(); } -template nothing operator+(nothing, T const &) { return nothing(); } -template nothing operator-(nothing, T const &) { return nothing(); } -template nothing operator*(nothing, T const &) { return nothing(); } -template nothing operator/(nothing, T const &) { return nothing(); } -template TYPE_DISABLE_IF(nothing, std::is_same) operator+(T const &, nothing) { return nothing(); } -template TYPE_DISABLE_IF(nothing, std::is_same) operator-(T const &, nothing) { return nothing(); } -template TYPE_DISABLE_IF(nothing, std::is_same) operator*(T const &, nothing) { return nothing(); } -template TYPE_DISABLE_IF(nothing, std::is_same) operator/(T const &, nothing) { return nothing(); } - -}} -#endif + template nothing slice_target(nothing, T...) { return nothing(); } + template nothing operator+(nothing, T const &) { return nothing(); } + template nothing operator-(nothing, T const &) { return nothing(); } + template nothing operator*(nothing, T const &) { return nothing(); } + template nothing operator/(nothing, T const &) { return nothing(); } + template TYPE_DISABLE_IF(nothing, std::is_same) operator+(T const &, nothing) { return nothing(); } + template TYPE_DISABLE_IF(nothing, std::is_same) operator-(T const &, nothing) { return nothing(); } + template TYPE_DISABLE_IF(nothing, std::is_same) operator*(T const &, nothing) { return nothing(); } + template TYPE_DISABLE_IF(nothing, std::is_same) operator/(T const &, nothing) { return nothing(); } +} +}