HilbertSpace.hpp

Go to the documentation of this file.

#pragma once
 
#include "Headers/mytypes.hpp"
#include "Headers/mydebug.hpp"
#include "Headers/Algorithm/powi.hpp"
#include "Headers/Algorithm/BaseNnumber.hpp"
#include "Headers/Algorithm/IntegerComposition.hpp"
#include "Headers/MatrixUtils/MatrixUtils.hpp"
#include "Eigen/SparseCore"
#include <iostream>
#include <vector>
#include <complex>
#include <cerrno>
 
#ifdef GPU
    #include "Headers/MatrixUtils/MatrixUtils.cuh"
    #include <thrust/host_vector.h>
    #include <thrust/device_vector.h>
    #include <thrust/execution_policy.h>
#endif
 
template<class T>
struct HilbertSpaceTraits;
 
template<class Derived = int>
class HilbertSpace;
template<class Derived>
class HilbertSpace {
    public:
        __host__ __device__ int dim() const {
            return static_cast<Derived const*>(this)->dim_impl();
        };
 
        /* @{ */
        __host__ __device__ bool operator==(HilbertSpace const& other) const {
            return this->dim() == other.dim();
        }
        /* @} */
};
 
template<>
class HilbertSpace<int> {
    private:
        int m_dim = 0;
 
    public:
        __host__ __device__ HilbertSpace(int dim) : m_dim{dim} { debug_constructor_printf(1); }
        __host__ __device__ HilbertSpace() = default;
        __host__ __device__ HilbertSpace(HilbertSpace const& other) = default;
        __host__ __device__ HilbertSpace(HilbertSpace&& other) = default;
        __host__ __device__ ~HilbertSpace() = default;
 
        /* @{ */
        __host__ __device__ HilbertSpace& operator=(HilbertSpace const& other) = default;
        __host__ __device__ HilbertSpace& operator=(HilbertSpace&& other) = default;
        __host__ __device__ bool operator==(HilbertSpace const& other) const {
            return this->dim() == other.dim();
        }
        /* @} */
 
        __host__ __device__ int dim() const { return m_dim; };
};
 
#ifdef GPU
// Forward declarelation for class ManyBodySpaceBase
template<typename Derived>
class ManyBodySpaceBase;
template<typename Derived>
__global__ void setTransEqDim_kernel(ManyBodySpaceBase<Derived>* obj, int transEqDim);
template<typename Derived>
__global__ void copyTransEqClass_kernel(ManyBodySpaceBase<Derived>* obj, int* srcTransEqClassRep,
                                        int* srcTransPeriod);
#endif
template<class Derived>
class ManyBodySpaceBase : public HilbertSpace<Derived> {
    private:
        using LocSpace_t = typename HilbertSpaceTraits<Derived>::LocalSpace;
#ifdef GPU
        mutable struct cudaFuncAttributes attr;
        mutable bool                      gotKernelAttribute = false;
#endif
 
    protected:
        using Vector_t = Eigen::VectorX<int>;
 
        int              m_systemSize = 0;  // Number of local Hilbert spaces
        LocSpace_t       m_locSpace;
        mutable int      m_transEqDim = -1;
        mutable Vector_t m_transEqClassRep;
        mutable Vector_t m_transPeriod;
        mutable Vector_t m_stateToTransEqClass;
        mutable Vector_t m_parityPair;
 
    public:
        __host__ __device__ ManyBodySpaceBase(int systemSize, LocSpace_t const& locSpace)
            : m_systemSize{systemSize}, m_locSpace{locSpace} {
            debug_constructor_printf(1);
            debug_printf("\tlocSpace.dim() = %d\n\n", locSpace.dim());
        }
        __host__ __device__ ManyBodySpaceBase(int systemSize, LocSpace_t&& locSpace)
            : m_systemSize{systemSize}, m_locSpace{std::move(locSpace)} {
            debug_constructor_printf(2);
            debug_printf("\tlocSpace.dim() = %d\n\n", locSpace.dim());
        }
 
        template<typename... Args>
        __host__ __device__ ManyBodySpaceBase(int systemSize, Args... args)
            : ManyBodySpaceBase(systemSize, LocSpace_t(args...)) {
            debug_constructor_printf(3);
        }
        __host__ __device__ ManyBodySpaceBase() = default;
        __host__ __device__ ManyBodySpaceBase(ManyBodySpaceBase const& other)
            : HilbertSpace<Derived>(other),
              m_systemSize(other.m_systemSize),
              m_locSpace(other.m_locSpace),
              m_transEqDim(other.m_transEqDim),
              m_transEqClassRep{other.m_transEqClassRep},
              m_transPeriod{other.m_transPeriod},
              m_stateToTransEqClass{other.m_stateToTransEqClass} {
            debug_constructor_printf((Copy));
        }
        __host__ __device__ ManyBodySpaceBase(ManyBodySpaceBase&& other)
            : HilbertSpace<Derived>(other),
              m_systemSize(other.m_systemSize),
              m_locSpace(other.m_locSpace),
              m_transEqDim(other.m_transEqDim),
              m_transEqClassRep{std::move(other.m_transEqClassRep)},
              m_transPeriod{std::move(other.m_transPeriod)},
              m_stateToTransEqClass{std::move(other.m_stateToTransEqClass)} {
            debug_constructor_printf((Move));
        }
        __host__ __device__ ~ManyBodySpaceBase() = default;
 
        /* @{ */
        __host__ __device__ ManyBodySpaceBase& operator=(ManyBodySpaceBase const& other) = delete;
        __host__ __device__ Derived& operator=(ManyBodySpaceBase&& other) {
            debug_printf("%s:\n\tMove assignment operator.\n\n", __PRETTY_FUNCTION__);
            m_systemSize          = other.m_systemSize;
            m_locSpace            = other.m_locSpace;
            m_transEqDim          = other.m_transEqDim;
            m_transEqClassRep     = std::move(other.m_transEqClassRep);
            m_transPeriod         = std::move(other.m_transPeriod);
            m_stateToTransEqClass = std::move(other.m_stateToTransEqClass);
            return *static_cast<Derived*>(this);
        }
        __host__ __device__ bool operator==(ManyBodySpaceBase const& other) const {
            return m_locSpace == other.locSpace() && m_systemSize == other.sysSize();
        }
        /* @} */
 
        __host__ __device__ LocSpace_t const& locSpace() const { return m_locSpace; }
        __host__ __device__ int               sysSize() const { return m_systemSize; }
        __host__ __device__ int               dimLoc() const { return m_locSpace.dim(); }
        __host__ __device__ int               locState(int state, int pos) const {
                          return static_cast<Derived const*>(this)->locState_impl(state, pos);
        }
 
        __host__ Derived const& printInString(int basisNum) const {
            Eigen::RowVectorXi config = this->ordinalToConfig(basisNum);
            for(int l = 0; l != this->sysSize(); ++l) std::cout << " " << config(l);
            std::cout << std::endl;
            return static_cast<Derived const&>(*this);
        }
        __host__ __device__ Eigen::RowVectorXi ordinalToConfig(int basisNum) const {
            return static_cast<Derived const*>(this)->ordinalToConfig_impl(basisNum);
        }
        template<class EigenDerived>
        __host__ __device__ int configToOrdinal(Eigen::DenseBase<EigenDerived>& config) const {
            return static_cast<Derived const*>(this)->configToOrdinal_impl(config);
        }
 
        /* @{ */
        __host__ __device__ int translate(int state, int trans) const {
            return static_cast<Derived const*>(this)->translate_impl(state, trans);
        }
        template<class EigenDerived>
        __host__ __device__ int translate(int state, int trans,
                                          Eigen::DenseBase<EigenDerived>& work) const {
            return static_cast<Derived const*>(this)->translate_impl(state, trans, work);
        }
 
        __host__ __device__ int transEqDim() const { return m_transEqDim; }
 
        __host__ __device__ Vector_t const& transEqClassRep() const { return m_transEqClassRep; }
        __host__ __device__ int transEqClassRep(int id) const { return m_transEqClassRep[id]; }
        __host__ __device__ int transEqClassRep(int id, int trans) const {
            return this->translate(m_transEqClassRep[id], trans);
        }
 
        __host__ __device__ Vector_t const& transPeriod() const { return m_transPeriod; }
        __host__ __device__ int             transPeriod(int id) const { return m_transPeriod[id]; }
 
        __host__ __device__ Vector_t const& stateToTransEqClass() const {
            return m_stateToTransEqClass;
        }
        __host__ __device__ int stateToTransEqClass(int state) const {
            return m_stateToTransEqClass[state];
        };
        __host__ __device__ int stateToTransPeriod(int state) const {
            auto eqClass = this->stateToTransEqClass(state);
            return this->transPeriod(eqClass);
        };
 
        __host__ __device__ void computeTransEqClass() const;
#ifdef GPU
        friend void setTransEqDim_kernel<Derived>(ManyBodySpaceBase*, int);
        friend void copyTransEqClass_kernel<Derived>(ManyBodySpaceBase*, int*, int*);
#endif
        /* @} */
 
        /* @{ */
        __host__ __device__ int reverse(int state) const {
            return static_cast<Derived const*>(this)->reverse_impl(state);
        }
        __host__ __device__ Vector_t const& parityPair() const { return m_parityPair; }
        __host__ __device__ int  parityPair(int state) const { return m_parityPair[state]; }
        __host__ __device__ void computeParityPair() const;
        /* @} */
};
 
template<class Derived>
__host__ void ManyBodySpaceBase<Derived>::computeTransEqClass() const {
    debug_printf("%s:\n\tComputing trans equivalent class... (this->dim()=%d, m_transEqDim=%d)\n",
                 __PRETTY_FUNCTION__, this->dim(), m_transEqDim);
    if(m_transEqDim != -1) {
        debug_printf("%s:\tTrans equivalent class is already computed. m_transEqDim=%d\n", __func__,
                     m_transEqDim);
        return;
    }
    if(this->dim() <= 0) {
        debug_printf("%s:\tDimension of the space is ZERO (this->dim() = %d)\n", __func__,
                     this->dim());
        return;
    }
    int                 translated;
    int                 period     = 0;
    Eigen::ArrayX<bool> calculated = Eigen::ArrayX<bool>::Zero(this->dim());
    m_transEqDim                   = 0;
    for(int basis = 0; basis < this->dim(); ++basis) {
        if(calculated(basis)) continue;
        calculated(basis) = true;
        for(period = 1; period < this->sysSize(); ++period) {
            translated = this->translate(basis, period);
            if(basis == translated) break;
            calculated(translated) = true;
        }
        m_transEqDim += 1;
    }
 
    m_transEqClassRep.resize(m_transEqDim);
    m_transPeriod.resize(m_transEqDim);
    m_stateToTransEqClass.resize(this->dim());
    int num    = 0;
    calculated = Eigen::ArrayX<bool>::Zero(this->dim());
    for(int basis = 0; basis < this->dim(); ++basis) {
        if(calculated(basis)) continue;
        calculated(basis)            = true;
        m_stateToTransEqClass(basis) = num;
 
        for(period = 1; period < this->sysSize(); ++period) {
            translated = this->translate(basis, period);
            if(basis == translated) break;
            calculated(translated)            = true;
            m_stateToTransEqClass(translated) = num;
        }
        m_transEqClassRep(num) = basis;
        m_transPeriod(num)     = period;
        num += 1;
    }
    debug_printf("\t(END) Computed trans equivalent classes (m_transEqDim=%d)\n\n", m_transEqDim);
}
 
template<class Derived>
__host__ __device__ void ManyBodySpaceBase<Derived>::computeParityPair() const {
    m_parityPair.resize(this->dim());
    for(auto j = 0; j != this->dim(); ++j) {
        auto reversed          = this->reverse(j);
        m_parityPair(j)        = reversed;
        m_parityPair(reversed) = j;
    }
}
 
class ManyBodySpinSpace;
template<>
struct HilbertSpaceTraits<ManyBodySpinSpace> {
        using LocalSpace = HilbertSpace<int>;
};
class ManyBodySpinSpace : public ManyBodySpaceBase<ManyBodySpinSpace> {
    private:
        BaseConverter<int> m_bConv;
 
    public:
        using Base       = ManyBodySpaceBase<ManyBodySpinSpace>;
        using LocalSpace = typename HilbertSpaceTraits<ManyBodySpinSpace>::LocalSpace;
        __host__ __device__ ManyBodySpinSpace(int systemSize, LocalSpace const& locSpace)
            : Base(systemSize, locSpace), m_bConv(locSpace.dim(), systemSize) {
            debug_constructor_printf(1);
        }
        __host__ __device__ ManyBodySpinSpace(int systemSize, LocalSpace&& locSpace)
            : Base(systemSize, std::move(locSpace)), m_bConv(locSpace.dim(), systemSize) {
            debug_constructor_printf(2);
        }
        __host__ __device__ ManyBodySpinSpace(int systemSize = 0, int dimLoc = 0)
            : Base(systemSize, LocalSpace(dimLoc)), m_bConv(dimLoc, systemSize) {
            debug_constructor_printf((Default));
        }
 
    private:
        friend HilbertSpace<ManyBodySpinSpace>;
        __host__ __device__ int dim_impl() const { return m_bConv.maxVal(); }
 
        friend ManyBodySpaceBase<ManyBodySpinSpace>;
        __host__ __device__ int locState_impl(int state, int pos) const {
            return m_bConv.digit(state, pos);
        }
 
        __host__ __device__ Eigen::RowVectorXi ordinalToConfig_impl(int basisNum) const {
            Eigen::RowVectorXi res;
            m_bConv.printInDigits(res, basisNum);
            return res;
        }
        template<class EigenDerived>
        __host__ __device__ int configToOrdinal_impl(Eigen::DenseBase<EigenDerived>& config) const {
            return m_bConv.digitsToNumber(config);
        }
 
        __host__ __device__ int translate_impl(int state, int trans) const {
            return m_bConv.shiftDigits(state, trans);
        }
        template<class EigenDerived>
        __host__ __device__ int translate_impl(int state, int trans,
                                               Eigen::DenseBase<EigenDerived>&) const {
            return this->translate_impl(state, trans);
        }
 
        __host__ __device__ int reverse_impl(int state) const {
            int res = 0, base = 1;
            for(int pos = 0; pos != this->sysSize(); ++pos, base *= this->dimLoc()) {
                res += base * this->locState(state, this->sysSize() - 1 - pos);
            }
            return res;
        }
};
 
class ManyBodyBosonSpace;
template<>
struct HilbertSpaceTraits<ManyBodyBosonSpace> {
        using LocalSpace = HilbertSpace<int>;
};
class ManyBodyBosonSpace : public ManyBodySpaceBase<ManyBodyBosonSpace> {
    private:
        IntegerComposition m_iComp;
 
    public:
        using Base       = ManyBodySpaceBase<ManyBodyBosonSpace>;
        using LocalSpace = typename HilbertSpaceTraits<ManyBodyBosonSpace>::LocalSpace;
        __host__ __device__ ManyBodyBosonSpace(int systemSize = 0, int nBosons = 0)
            : Base(systemSize, LocalSpace(nBosons + 1)), m_iComp(nBosons, systemSize, nBosons) {
            debug_constructor_printf((Default));
        }
        __host__ __device__ ManyBodyBosonSpace(int systemSize, int nBosons, int max)
            : Base(systemSize, LocalSpace(max + 1)), m_iComp(nBosons, systemSize, max) {
            debug_constructor_printf(1);
        }
 
    private:
        friend HilbertSpace<ManyBodyBosonSpace>;
        __host__ __device__ int dim_impl() const { return m_iComp.dim(); }
 
        friend ManyBodySpaceBase<ManyBodyBosonSpace>;
        __host__ __device__ int locState_impl(int state, int pos) const {
            return m_iComp.locNumber(state, pos);
        }
 
        __host__ __device__ Eigen::RowVectorXi ordinalToConfig_impl(int basisNum) const {
            Eigen::RowVectorXi res;
            m_iComp.ordinalToPart(res, basisNum);
            return res;
        }
        template<class EigenDerived>
        __host__ __device__ int configToOrdinal_impl(Eigen::DenseBase<EigenDerived>& config) const {
            return m_iComp.partToOrdinal(config);
        }
 
        __host__ __device__ int translate_impl(int state, int trans) const {
            return m_iComp.translate(state, trans);
        }
        template<class EigenDerived>
        __host__ __device__ int translate_impl(int state, int trans,
                                               Eigen::DenseBase<EigenDerived>&) const {
            return this->translate_impl(state, trans);
        }
 
        __host__ __device__ int reverse_impl(int state) const {
            Eigen::ArrayXi config1(this->sysSize());
            m_iComp.ordinalToPart(config1, state);
            Eigen::ArrayXi config2(this->sysSize());
            for(int pos = 0; pos != this->sysSize(); ++pos)
                config2[pos] = config1[this->sysSize() - 1 - pos];
            return m_iComp.partToOrdinal(config2);
        }
};
 
template<class TotalSpace_, typename ScalarType_ = Complex_t<Real_t>>
class SubSpace : public HilbertSpace< SubSpace<TotalSpace_, ScalarType_> > {
    public:
        using TotalSpace = TotalSpace_;
        using Scalar     = ScalarType_;
        using Real       = typename Eigen::NumTraits<ScalarType_>::Real;
 
    protected:
#ifdef __CUDA_ARCH__
        using Matrix_t = SparseCompressed<Scalar>;
#else
        using Matrix_t = Eigen::SparseMatrix<Scalar>;  // default is column major
#endif
        TotalSpace m_totalSpace;
        Matrix_t   m_basisStates;
 
    public:
        __host__ __device__ SubSpace(TotalSpace const& totalSpace)
            : m_totalSpace{totalSpace}, m_basisStates(totalSpace.dim(), 0) {
            debug_constructor_printf(1);
        }
        __host__ __device__ SubSpace(TotalSpace&& totalSpace)
            : m_totalSpace{std::move(totalSpace)}, m_basisStates(totalSpace.dim(), 0) {
            debug_constructor_printf(2);
        }
        __host__ __device__ SubSpace() = default;
        __host__ __device__ SubSpace(SubSpace const& other) = default;
        template<typename Scalar_>
        __host__ __device__ SubSpace(SubSpace<TotalSpace_, Scalar_> const& other)
            : m_totalSpace{other.totalSpace()},
              m_basisStates{other.basis().template cast<Scalar>()} {};
        __host__ __device__ SubSpace(SubSpace&& other) = default;
        __host__ __device__ ~SubSpace() = default;
 
        /* @{ */
        __host__ __device__ SubSpace& operator=(SubSpace const& other) {
            debug_printf("%s:\n\tCopy assignment operator.\n\n", __PRETTY_FUNCTION__);
            m_totalSpace  = other.m_totalSpace;
            m_basisStates = other.m_basisStates;
            return *this;
        }
        __host__ __device__ SubSpace& operator=(SubSpace&& other) {
            debug_printf("%s:\n\tMove assignment operator.\n\n", __PRETTY_FUNCTION__);
            m_totalSpace  = std::move(other.m_totalSpace);
            m_basisStates = std::move(other.m_basisStates);
            return *this;
        }
        __host__ __device__ bool operator==(SubSpace const& other) const {
            return m_totalSpace == other.totalSpace() && m_basisStates == other.basis();
        }
        /* @} */
 
        __host__ __device__ Matrix_t&       basis() { return m_basisStates; }
        __host__ __device__ Matrix_t const& basis() const { return m_basisStates; }
 
        __host__ __device__ TotalSpace const& totalSpace() const { return m_totalSpace; }
        __host__ __device__ int               dimTot() const { return m_totalSpace.dim(); }
 
    private:
        friend HilbertSpace< SubSpace<TotalSpace, Scalar> >;
        __host__ __device__ int dim_impl() const { return m_basisStates.cols(); }
};