.. _Operator:

########
Operator
########



.. admonition:: Summary

    *   The class can be built from the :ref:`State class<State>` and the :ref:`Basis class<BasisU1>` with number conservation.

    *   The interface is similar to the `ITensor <https://itensor.org/>`_ package and allows a straightforward usage.

    *   It provides different formats to extract the matrix such as :ref:`Eigen<Eigen_setup>` and :ref:`Petsc<Petsc_setup>`.

    *   It handles Hamiltonian and Lindbladian systems

.. _Operator_Type:

=============
Operator Type
=============

The operator can be defined on space of wavefunction **and** on the space of the density matrices.
The latter is required when we deal with dissipative system and the dynamics is described the `Lindblad <https://en.wikipedia.org/wiki/Lindbladian>`_ formalism.

Both can be used **with** and **without** number conservation and all four cases are predefined in the header file :ref:`DanceQ.h<Examples_DanceQ_header>` (**L**, **n** refer to the system size, particle number -- if desired):

    * ``danceq::Lindbladian_U1(L,n)``

    * ``danceq::Lindbladian_NoU1(L)``

    * ``danceq::Hamiltonian_U1(L,n)``

    * ``danceq::Hamiltonian_NoU1(L)``

The class :ref:`Operator class<Operator_Autodocs>` has three template parameters (Basis, ScalarType, and OpertorType).
The OperatorType is a ``enum`` that distinguishes the two different cases:

.. _Operator_OperatorType:

.. doxygenenum:: danceq::OperatorType
    :project: DanceQ


The Lindbladian acts on the space of density matrices :math:`\rho`:

.. math::

    \frac{\text{d}\rho}{\text{d}t} = -i[H,\rho] + \sum_{l} l\rho l^\dagger - \gamma\frac{1}{2}l^\dagger l\rho - \gamma\frac{1}{2}\rho l^\dagger l

.. _Operator_Basis:

=====
Basis
=====

The underlying basis of the operator can be set to the :ref:`State class <State>` with dimension :math:`\mathrm{Q}^L` or to the :ref:`Basis class<Basis_General>` with dimension from :ref:`Eq. (2)<eq:dim>`:

.. code-block:: cpp
    
    /* Operator using the Basis class, system size L, particle number n */
    using MyOperatorU1=danceq::internal::Operator<Basis,double,Hamiltonian> Operator;
    MyOperatorU1 H_U1(L,n);

    /* Operator using the State class, system size L */
    using MyOperatorNoU1=danceq::internal::Operator<State,double,Hamiltonian>;
    MyOperatorNoU1 H_NoU1(L);

The second template parameter defines the scalar, *e.g.*, **complex** and **real**.

.. warning::

    **MaxSites** defined in the :ref:`State class<State>` has to be larger or equal the system size if we deal with Hamiltonian systems and at least double the size if we deal Lindbladian systems.



.. _Operator_Input:

=====
Input
=====

The interface is similar to the `ITensor <https://itensor.org/>`_ and provides an easy way to create arbitrary linear operators from tensor products with the :ref:`add_operator(...)<Operator_Functions_define_operator>` function.
The following example creates a one-dimensional Heisenberg model of size **L** with **n** particles, coupling constant **J**, a local alternating field **h**, and a dissipation rate **gamma**:

.. code-block:: cpp

    /* system size */
    uint64_t L = 10; 

    /* particle number */
    uint64_t n = 5;

    /* field strength */
    double h = 0.1;

    /* coupling strength */
    double J = 1.;

    /* dissipation rate */
    double gamma = 0.01;

    /* defining the Lindbladian with number conservation */
    Lindbladian_U1(L,n);

    /* adding the Hamiltonian terms */
    for(uint64_t i = 0; i < L; i++){
        L.add_operator(J/2., {i,(i+1)%L}, {"S+","S-"});
        L.add_operator(J/2., {i,(i+1)%L}, {"S-","S+"});
        L.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
        L.add_operator(std::pow(-h,i), {i}, {"Sz"});
    }

    /* jump operator acting on site 0 */
    L.add_jump_operator(gamma, {0}, {"Sz"});

 
A list of :ref:`default operators<Operator_Functions_add_default_operators>` is already defined. 
New local operators of size :math:`\textbf{Q}\times\textbf{Q}` acting on a single site can be defined using :ref:`define_operator(...)<Operator_Functions_define_operator>` with a ``std::array< std::array<ScalarType, Q>, Q>``.


.. _Operator_formats:

==============
Matrix formats
==============

The class comes with different formats to represent the operator as a matrix using the provided :ref:`Basis<Operator_Basis>`.
Besides native **cpp** containers, frontends for :ref:`Eigen<Eigen_setup>` and :ref:`Petsc<Petsc_setup>` are available if they are included by :ref:`CMake<Config_and_CMake>`

The dense matrix can be retrieved using a two-dimensional ``std::vector<std::vector<ScalarType>>``.
We further provided an advanced :ref:`SparseMatrix calss<SparseMatrix>` that can be used for simple routines such as **Matrix-Vector** multiplication.
The class works together with the provided :ref:`Vector class<Vector>` that is a primitive extension of ``std::vector``.
It is capable of primitive routines and provides distributed memory if :ref:`MPI<MPI_docs>` is available.
Together with the :ref:`SparseMatrix calss<SparseMatrix>` naive `Krylov space methods <https://en.wikipedia.org/wiki/Krylov_subspace>`_ can be used such as ground state calculation.
No additional third-party libraries like :ref:`Eigen<Eigen_setup>` or :ref:`Petsc and Slepc<Petsc_setup>` are required.
For more complex operators, we recommend working with more advanced linear algebra packages as listed below:

:ref:`Eigen<Eigen_setup>` can be used for dense applications (*e.g.*, calculating the full spectrum) and provides access to `LAPACK <https://de.wikipedia.org/wiki/LAPACK>`_ libraries such as `Intel's MKL <https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html>`_.
For sparse operations (*e.g.*, `Krylov space <https://en.wikipedia.org/wiki/Krylov_subspace>`_ methods), :ref:`Petsc and Slepc<Petsc_setup>` provide a powerful :ref:`MPI-based<MPI_docs>` approach allowing scalable computation across multiple nodes!
Additionally, we have added a **matrix-free** implementation with :ref:`Petsc<Petsc_setup>` to reduce the memory capacity. 
The application does not store the matrix and computes its elements **on-the-fly**. 




.. admonition:: Available matrix formats

    * :ref:`create_DenseMatrix()<Operator_Functions_Frontends_create_DenseMatrix>`: Dense matrix using two-dimensional ``std::vector``

    * :ref:`create_EigenDenseMatrix()<Operator_Functions_Frontends_create_EigenDenseMatrix>`: Dense matrix using `Eigen <https://eigen.tuxfamily.org/index.php?title=Main_Page>`_

    * :ref:`create_SparseMatrix()<Operator_Functions_Frontends_create_SparseMatrix>`: Sparse matrix using the :ref:`SparseMatrix class<SparseMatrix>`

    * :ref:`create_ShellMatrix()<Operator_Functions_Frontends_create_SparseMatrix>`: **Matrix-free** multiplication using the :ref:`ShellMatrix class<SparseMatrix>`

    * :ref:`create_PetscSparseMatrix()<Operator_Functions_Frontends_create_PetscSparseMatrix>`: Sparse matrix using `Petsc <https://petsc.org/release>`_

    * :ref:`create_PetscShellMatrix()<Operator_Functions_Frontends_create_PetscShellMatrix>`: **Matrix-free** multiplication using `Petsc <https://petsc.org/release>`_ with the :ref:`ShellMatrix class<ShellMatrix>`


===============================
Application to a product states
===============================

The core function which is used to create any matrix format is :ref:`apply(...)<Operator_Functions_apply>`.
It returns the elements (indices and scalars) of **i** th column when applied to in the **i** th product state of the underlying basis.

In order to optimize the computation, we exploit a similar divide-and-conquer ansatz acting on non-local subclusters.
When a new operator string is added by :ref:`add_operator(...)<Operator_Functions_add_operator>` it is only stored in :ref:`operator_list<Operator_Members_operator_list>` and is **not** directly used for the computation.
The :ref:`apply(...)<Operator_Functions_apply>` function (and therefore all other more complex functions) is working with the data stored in :ref:`sparse_matrices_for_computation<Operator_Members_operator_sparse_matrices_for_computation>` which is created by the :ref:`merge_terms(...)<Operator_Functions_merge_terms>`.
It bundles all operator strings from :ref:`operator_list<Operator_Members_operator_list>` in :ref:`sparse_matrices_for_computation<Operator_Members_operator_sparse_matrices_for_computation>`. 
**You** do not need to call :ref:`merge_terms(...)<Operator_Functions_merge_terms>` since it is called every time a new operator string was added by :ref:`add_operator(...)<Operator_Functions_add_operator>`.


The :ref:`merge_terms(...)<Operator_Functions_merge_terms>` function identifies clusters of maximal size :ref:`Ntensor_max<Operator_Members_Operator_Information_Ntensor_max>` and integrates all operator strings that are fully supported on such a cluster.
Their additive contribution is stored in a sparse format -- similar to the :ref:`SparseMatrix<SparseMatrix>` --  within :ref:`sparse_matrices_for_computation<Operator_Members_operator_sparse_matrices_for_computation>`.
The default value of :ref:`Ntensor_max<Operator_Members_Operator_Information_Ntensor_max>` is chosen such that the dimension of the underlying sparse matrix does not exceed 2048.
:ref:`Ntensor_max<Operator_Members_Operator_Information_Ntensor_max>` can be set by calling :ref:`set_Ntensor_max(...)<Operator_Functions_set_Ntensor_max>`.

============
Full example
============

.. code-block:: cpp
   :linenos:

    #include "Operator.h"

    /* Defining the state class with a maximal length of 256 sites and a local dimension of 4 */
    using State=danceq::internal::State<256,4>; 

    /* Defining the Basis class from the State class */
    using Basis=danceq::internal::BasisU1<State>; 

    /* Defining the Operator from the Basis class */
    using MyOperator=danceq::internal::Operator<Basis,double,Hamiltonian>;

    /* Defining the sparse matrix class using the scalar from the Operator class */
    using SparseMatrix=danceq::internal::SparseMatrix<Operator::scalartype>;


    int main(int argc, char *argv[]) {

        /* Input parameter */
        uint64_t L = 10;
        uint64_t n = 5;
        double h = 0.1;
        double J = 1.;

        /* Operator class */
        MyOperator H(L,n);

        uint64_t dim = H.get_dim();
        std::cout << "[main] - dim = " << dim << std::endl;

        /* Simply XXX-Heisenberg model with a local field along the z-axis */
        for(uint64_t i = 0; i < L; i++){
            H.add_operator(J/2, {i,(i+1)%L}, {"S+","S-"});
            H.add_operator(J/2, {i,(i+1)%L}, {"S-","S+"});
            H.add_operator(J, {i,(i+1)%L}, {"Sz","Sz"});
            H.add_operator(std::pow(-h,i), {i}, {"Sz"});
        }

        /* Dense matrix using std::vector<std::vector<Operator::scalartype>> */
        auto H_dense = H.create_DenseMatrix();

        /* Sparse matrix using Sprase_Matrix class */
        auto H_sparse = H.create_SparseMatrix();

        /* Dense matrix with Eigen using the type Operator::EigenMatrixType */
        auto H_eigen = H.create_EigenDenseMatrix();

        /* Sparse matrix with Petsc Mat type MATMPIAIJ */
        auto H_petsc = H.create_PetscSparseMatrix();

        /* Sparse matrix with Petsc using a matrix-free shell */
        auto H_petsc_shell = H.create_PetscShellMatrix();

        return 0;
    }


=======
Content
=======

Detailed documentation for the class can be found in :ref:`Data Members<Operator_Members>`, :ref:`Functions<Operator_Functions>`, and :ref:`Autodocs Operator<Operator_Autodocs>`.

.. toctree::
    :maxdepth: 1
    :titlesonly:

    Operator/Members.rst
    Operator/Constructors.rst
    Operator/Functions.rst
    Operator/Autodocs.rst