OperatorSpace.hpp

Go to the documentation of this file.

 
#pragma once
 
#define _USE_MATH_DEFINES
#include "Headers/mytypes.hpp"
#include "Headers/mydebug.hpp"
#include "Headers/HilbertSpace/HilbertSpace.hpp"
#include "Headers/Algorithm/IntegerComposition.hpp"
#include <cmath>
 
template<class Derived>
class OperatorSpaceBase {
    private:
        using BaseSpace  = typename HilbertSpaceTraits<Derived>::BaseSpace;
        using ScalarType = typename HilbertSpaceTraits<Derived>::Scalar;
        BaseSpace m_baseSpace;
 
    public:
        typedef typename Eigen::NumTraits<ScalarType>::Real Real;
 
        __host__ __device__ OperatorSpaceBase(BaseSpace const& baseSpace) : m_baseSpace{baseSpace} {
            debug_constructor_printf(1(Custom Copy));
        }
        __host__ __device__ OperatorSpaceBase(BaseSpace&& baseSpace)
            : m_baseSpace{std::move(baseSpace)} {
            debug_constructor_printf(2(Custom Move));
        }
        __host__ __device__ OperatorSpaceBase() = default;
        __host__ __device__ OperatorSpaceBase(OperatorSpaceBase const& other) = default;
        __host__ __device__ OperatorSpaceBase(OperatorSpaceBase&& other) = default;
        __host__ __device__ ~OperatorSpaceBase() = default;
 
        /* @{ */
        __host__ __device__ OperatorSpaceBase& operator=(OperatorSpaceBase const& other) = default;
        __host__ __device__ OperatorSpaceBase& operator=(OperatorSpaceBase&& other) = default;
        /* @} */
 
        __host__ __device__ int dim() const {
            return static_cast<Derived const*>(this)->dim_impl();
        }
 
        __host__ __device__ BaseSpace const& baseSpace() const { return m_baseSpace; }
        __host__ __device__ int              baseDim() const { return m_baseSpace.dim(); }
 
        template<typename RealType>
        __host__ __device__ void action(int& resBasisNum, Complex_t<RealType>& coeff, int opNum,
                                        int basisNum) const {
            static_cast<Derived const*>(this)->action_impl(resBasisNum, coeff, opNum, basisNum);
        };
 
        template<typename RealType, class EigenDerived>
        __host__ __device__ void action(int& resBasisNum, Complex_t<RealType>& coeff, int opNum,
                                        int                             basisNum,
                                        Eigen::DenseBase<EigenDerived>& config) const {
            static_cast<Derived const*>(this)->action_impl(resBasisNum, coeff, opNum, basisNum,
                                                           config);
        };
 
        __host__ __device__ std::pair<int, Complex_t<Real>> action(int opNum, int basisNum) const {
            // debug_print(__PRETTY_FUNCTION__);
            int             resBasisNum;
            Complex_t<Real> coeff;
            this->action(resBasisNum, coeff, opNum, basisNum);
            return std::make_pair(resBasisNum, coeff);
        }
 
        __host__ void basisOp(Eigen::SparseMatrix<ScalarType>& res, int opNum) const {
            res.resize(this->baseDim(), this->baseDim());
#pragma omp parallel for
            for(int basisNum = 0; basisNum < this->baseDim(); basisNum++) {
                auto [resBasisNum, coeff]           = this->action(opNum, basisNum);
                res.coeffRef(resBasisNum, basisNum) = coeff;
            };
        }
 
        __host__ Eigen::SparseMatrix<ScalarType> basisOp(int opNum) const {
            Eigen::SparseMatrix<ScalarType> res(this->baseDim(), this->baseDim());
            for(int basisNum = 0; basisNum < this->baseDim(); basisNum++) {
                // debug_print("\t basisNum=" << basisNum);
                auto [resBasisNum, coeff]         = this->action(opNum, basisNum);
                res.insert(resBasisNum, basisNum) = coeff;
            }
            return res;
        }
};
 
template<typename ScalarType>
class OperatorSpace;
template<typename ScalarType>
struct HilbertSpaceTraits< OperatorSpace<ScalarType> > {
        using BaseSpace = HilbertSpace<int>;
        using Scalar    = ScalarType;
};
template<typename ScalarType = Complex_t<Real_t>>
class OperatorSpace : public OperatorSpaceBase< OperatorSpace<ScalarType> > {
    public:
        using OperatorSpaceBase< OperatorSpace >::OperatorSpaceBase;
        using BaseSpace = typename HilbertSpaceTraits<OperatorSpace>::BaseSpace;
        using Scalar    = ScalarType;
 
    private:
        friend HilbertSpace< OperatorSpace >;
        __host__ __device__ int dim_impl() const {
            return this->baseSpace().dim() * this->baseSpace().dim();
        }
 
        friend OperatorSpaceBase< OperatorSpace >;
        template<typename RealType>
        __host__ __device__ void action_impl(int& resBasisNum, Complex_t<RealType>& coeff,
                                             int opNum, int basisNum) const {
            // debug_printf("%s\n", __PRETTY_FUNCTION__);
            resBasisNum = basisNum;
            coeff       = {0.0, 0.0};
            int Digit1, Digit2;
 
            if(opNum < this->baseSpace().dim()) {
                resBasisNum = basisNum;
                if(opNum == 0) { coeff = 1.0 / RealType(sqrt(RealType(this->baseSpace().dim()))); }
                else if(opNum == 1 && this->baseSpace().dim() % 2 == 0) {
                    coeff = (basisNum % 2 == 0 ? 1.0 : -1.0)
                            / RealType(sqrt(RealType(this->baseSpace().dim())));
                }
                else {
                    coeff = sin(M_PI
                                * (2 * (opNum / 2) * basisNum / RealType(this->baseSpace().dim())
                                   + (opNum % 2) / 2.0))
                            / RealType(sqrt(RealType(this->baseSpace().dim()) / 2.0));
                    ;
                }
            }
            else if(opNum < this->baseSpace().dim() * (this->baseSpace().dim() + 1) / 2) {
                opNum -= this->baseSpace().dim();
                Digit1 = int(sqrt(RealType(2 * opNum) + 0.25) + 0.5);
                Digit2 = opNum - Digit1 * (Digit1 - 1) / 2;
                if(Digit1 == basisNum) {
                    resBasisNum = Digit2;
                    coeff       = Complex_t<RealType>(1.0 / sqrt(RealType(2)), 0.0);
                }
                else if(Digit2 == basisNum) {
                    resBasisNum = Digit1;
                    coeff       = Complex_t<RealType>(1.0 / sqrt(RealType(2)), 0.0);
                }
            }
            else if(opNum < this->baseSpace().dim() * this->baseSpace().dim()) {
                opNum -= this->baseSpace().dim() * (this->baseSpace().dim() + 1) / 2;
                Digit1 = int(sqrt(RealType(2 * opNum) + 0.25) + 0.5);
                Digit2 = opNum - Digit1 * (Digit1 - 1) / 2;
                if(Digit1 == basisNum) {
                    resBasisNum = Digit2;
                    coeff       = Complex_t<RealType>(0.0, -1.0 / sqrt(RealType(2)));
                }
                else if(Digit2 == basisNum) {
                    resBasisNum = Digit1;
                    coeff       = Complex_t<RealType>(0.0, +1.0 / sqrt(RealType(2)));
                }
            }
            else {
                // clang-format off
                #ifndef __CUDA_ARCH__
                std::cerr << "Error: " << __PRETTY_FUNCTION__ << "\n\t"
                          << " opNum(" << opNum << ") must be less than "
                          << this->baseSpace().dim() * this->baseSpace().dim() << std::endl;
                std::exit(EXIT_FAILURE);
                #endif
                // clang-format on
            }
            return;
        };
};
 
template<class Derived>
class ManyBodyOperatorSpaceBase : public OperatorSpaceBase<Derived>,
                                  public ManyBodySpaceBase<Derived> {
    public:
        using ManyBodySpaceBase<Derived>::dim;
        using BaseSpace  = typename HilbertSpaceTraits< Derived >::BaseSpace;
        using LocalSpace = typename HilbertSpaceTraits< Derived >::LocalSpace;
 
        __host__ __device__ ManyBodyOperatorSpaceBase(BaseSpace const& baseSpace, int systemSize,
                                                      LocalSpace const& locSpace)
            : OperatorSpaceBase<Derived>(baseSpace),
              ManyBodySpaceBase<Derived>(systemSize, locSpace) {
            debug_constructor_printf(1(Custom Copy));
        }
        __host__ __device__ ManyBodyOperatorSpaceBase(BaseSpace&& baseSpace, int systemSize,
                                                      LocalSpace&& locSpace)
            : OperatorSpaceBase<Derived>(std::move(baseSpace)),
              ManyBodySpaceBase<Derived>(systemSize, std::move(locSpace)) {
            debug_constructor_printf(2(Custom Move));
        }
 
        __host__ __device__ ManyBodyOperatorSpaceBase() = default;
        __host__ __device__ ManyBodyOperatorSpaceBase(ManyBodyOperatorSpaceBase const& other) = default;
        __host__ __device__ ManyBodyOperatorSpaceBase(ManyBodyOperatorSpaceBase&& other) = default;
        __host__ __device__ ~ManyBodyOperatorSpaceBase() = default;
 
        /* @{ */
        __host__ __device__ ManyBodyOperatorSpaceBase& operator=(ManyBodyOperatorSpaceBase const& other) = default;
        __host__ __device__ ManyBodyOperatorSpaceBase& operator=(ManyBodyOperatorSpaceBase&& other) = default;
        /* @} */
};
 
template<typename ScalarType = Complex_t<Real_t>>
class ManyBodyOperatorSpace;
template<typename ScalarType>
struct HilbertSpaceTraits< ManyBodyOperatorSpace<ScalarType> > {
        using BaseSpace  = ManyBodySpinSpace;
        using LocalSpace = OperatorSpace<ScalarType>;
        using Scalar     = ScalarType;
};
template<typename ScalarType>
class ManyBodyOperatorSpace
    : public ManyBodyOperatorSpaceBase< ManyBodyOperatorSpace<ScalarType> > {
    private:
        BaseConverter<int> m_bConv;
 
    public:
        using Scalar = ScalarType;
        using Base   = ManyBodyOperatorSpaceBase< ManyBodyOperatorSpace<ScalarType> >;
        using BaseSpace =
            typename HilbertSpaceTraits< ManyBodyOperatorSpace<ScalarType> >::BaseSpace;
        using LocalSpace =
            typename HilbertSpaceTraits< ManyBodyOperatorSpace<ScalarType> >::LocalSpace;
 
        __host__ __device__ ManyBodyOperatorSpace(BaseSpace const& baseSpace)
            : Base(baseSpace, baseSpace.sysSize(), OperatorSpace<ScalarType>(baseSpace.locSpace())),
              m_bConv(baseSpace.dimLoc() * baseSpace.dimLoc(), baseSpace.sysSize()) {
            debug_constructor_printf(1(Custom Copy));
        }
        __host__ __device__ ManyBodyOperatorSpace(BaseSpace&& baseSpace)
            : Base(std::move(baseSpace), baseSpace.sysSize(),
                   OperatorSpace<ScalarType>(baseSpace.locSpace())),
              m_bConv(baseSpace.dimLoc() * baseSpace.dimLoc(), baseSpace.sysSize()) {
            debug_constructor_printf(2(Custom Move));
        }
        __host__ __device__ ManyBodyOperatorSpace() = default;
        __host__ __device__ ManyBodyOperatorSpace(ManyBodyOperatorSpace const& other) = default;
        __host__ __device__ ManyBodyOperatorSpace(ManyBodyOperatorSpace&& other) = default;
        __host__ __device__ ~ManyBodyOperatorSpace() = default;
 
        /* @{ */
        __host__ __device__ ManyBodyOperatorSpace& operator=(ManyBodyOperatorSpace const& other) = default;
        __host__ __device__ ManyBodyOperatorSpace& operator=(ManyBodyOperatorSpace&& other) = default;
        /* @} */
 
    private:
        friend HilbertSpace< ManyBodyOperatorSpace >;
        __host__ __device__ int dim_impl() const {
            return this->baseSpace().dim() * this->baseSpace().dim();
        }
 
        /* @{ */
        friend OperatorSpaceBase< ManyBodyOperatorSpace >;
        template<typename RealType>
        __host__ __device__ void action_impl(int& resBasisNum, Complex_t<RealType>& coeff,
                                             int opNum, int basisNum) const {
            // debug_print(__func__ << " at Class:ManyBodyOperatorSpace");
            resBasisNum = basisNum;
            coeff       = {1, 0};
            int locOpNum, locBasisNum;
            int base = 1;
            for(int pos = 0; pos != this->sysSize(); ++pos, base *= this->baseSpace().dimLoc()) {
                locOpNum                           = this->locState(opNum, pos);
                locBasisNum                        = this->baseSpace().locState(basisNum, pos);
                auto [resLocBasisNum, resLocCoeff] = this->locSpace().action(locOpNum, locBasisNum);
                if(abs(resLocCoeff) == 0) {
                    coeff       = 0;
                    resBasisNum = basisNum;
                    break;
                }
                coeff *= resLocCoeff;
                resBasisNum += (resLocBasisNum - locBasisNum) * base;
            }
        }
        /* @} */
 
        /* @{ */
        friend ManyBodySpaceBase< ManyBodyOperatorSpace >;
        __host__ __device__ int locState_impl(int state, int pos) const {
            return this->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);
        }
 
        template<class EigenDerived>
        __host__ __device__ int translate_impl(int state, int trans) const {
            return this->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);
        }
        /* @} */
};
 
template<typename ScalarType = Complex_t<Real_t>>
class mBodyOperatorSpace;
template<typename ScalarType>
struct HilbertSpaceTraits< mBodyOperatorSpace<ScalarType> > {
        using BaseSpace  = ManyBodySpinSpace;
        using LocalSpace = OperatorSpace<ScalarType>;
        using Scalar     = ScalarType;
};
template<typename ScalarType>
class mBodyOperatorSpace : public ManyBodyOperatorSpaceBase< mBodyOperatorSpace<ScalarType> > {
    private:
        template<class T>
        using Vector_t = Eigen::VectorX<T>;
        template<class T>
        using Array_t = Eigen::ArrayX<T>;
 
        int                m_mBody     = 0;
        int                m_opConfDim = 0;
        IntegerComposition m_ActingSites;
        BaseConverter<int> m_opConfig;
#ifndef __CUDA_ARCH__
        mutable Array_t< Array_t<int> > m_config;
#endif
 
    public:
        using Scalar = ScalarType;
        using Base      = ManyBodyOperatorSpaceBase< mBodyOperatorSpace<ScalarType> >;
        using BaseSpace = typename HilbertSpaceTraits< mBodyOperatorSpace<ScalarType> >::BaseSpace;
        using LocalSpace =
            typename HilbertSpaceTraits< mBodyOperatorSpace<ScalarType> >::LocalSpace;
        using ManyBodySpaceBase< mBodyOperatorSpace<ScalarType> >::dim;
 
        __host__ __device__ mBodyOperatorSpace(int m, BaseSpace const& baseSpace)
            : Base(baseSpace, baseSpace.sysSize(), OperatorSpace<ScalarType>(baseSpace.locSpace())),
              m_mBody{m},
              m_opConfDim(powi(baseSpace.dimLoc() * baseSpace.dimLoc() - 1, m)),
              m_ActingSites(m, baseSpace.sysSize(), 1),
              m_opConfig(baseSpace.dimLoc() * baseSpace.dimLoc() - 1, m) {
            debug_constructor_printf(1);
#ifndef __CUDA_ARCH__
            if(m > baseSpace.sysSize()) {
                std::cerr << "ERROR: " << __PRETTY_FUNCTION__
                          << "\nMessage:\t m(" << m << ") > "
                          << "baseSpace.sysSize()(" << baseSpace.sysSize() << ")" << std::endl;
                std::exit(EXIT_FAILURE);
            }
            m_config.resize(get_max_threads());
            for(int j = 0; j != get_max_threads(); ++j) m_config[j].resize(this->sysSize());
#endif
        };
        __host__ __device__ mBodyOperatorSpace(int m, int sysSize, int dimLoc)
            : mBodyOperatorSpace(m, BaseSpace(sysSize, dimLoc)) {
            debug_constructor_printf(2);
        };
        __host__ __device__ mBodyOperatorSpace() = default;
        __host__ __device__ mBodyOperatorSpace(mBodyOperatorSpace const& other) = default;
        __host__ __device__ mBodyOperatorSpace(mBodyOperatorSpace&& other) = default;
        __host__ __device__ ~mBodyOperatorSpace() = default;
 
        /* @{ */
        __host__ __device__ mBodyOperatorSpace& operator=(mBodyOperatorSpace const& other) = default;
        __host__ __device__ mBodyOperatorSpace& operator=(mBodyOperatorSpace&& other) = default;
        /* @} */
 
 
 
        __host__ __device__ int m() const { return m_mBody; }
 
    private:
        friend HilbertSpace< mBodyOperatorSpace >;
        __host__ __device__ int dim_impl() const { return m_ActingSites.dim() * m_opConfDim; }
 
        friend OperatorSpaceBase< mBodyOperatorSpace >;
        template<typename RealType>
        __host__ __device__ void action_impl(int& resBasisNum, Complex_t<RealType>& coeff,
                                             int opOrdinal, int basisNum) const {
#ifndef __CUDA_ARCH__
            action_impl(resBasisNum, coeff, opOrdinal, basisNum, m_config[get_thread_num()]);
#else
            Eigen::ArrayXi config(this->sysSize());
            action_impl(resBasisNum, coeff, opOrdinal, basisNum, config);
#endif
        };
 
        template<typename RealType, class Derived>
        __host__ __device__ void action_impl(int& resBasisNum, Complex_t<RealType>& coeff,
                                             int opOrdinal, int basisNum,
                                             Eigen::DenseBase<Derived>& config) const;
 
        friend ManyBodySpaceBase< mBodyOperatorSpace >;
        __host__ __device__ int locState_impl(int opOrdinal, int pos) const {
            int posConfNum = opOrdinal / m_opConfDim;
            if(this->m_ActingSites.locNumber(posConfNum, pos) == 0)
                return 0;
            else {
                int index = -1;
                for(int l = 0; l <= (pos % this->m_systemSize); ++l) {
                    index += this->m_ActingSites.locNumber(posConfNum, l);
                }
                int opConfNum = opOrdinal % m_opConfDim;
                return this->m_opConfig.digit(opConfNum, index) + 1;
            }
        }
 
        __host__ __device__ int translate_impl(int opOrdinal, int trans) const {
            Eigen::ArrayXi work(this->sysSize());
            return this->translate(opOrdinal, trans, work);
        };
        template<class EigenDerived>
        __host__ __device__ int translate_impl(int opOrdinal, int trans,
                                               Eigen::DenseBase<EigenDerived>& workSpace) const;
 
        __host__ __device__ Eigen::RowVectorXi ordinalToConfig_impl(int basisNum) const;
        template<class EigenDerived>
        __host__ __device__ int configToOrdinal_impl(Eigen::DenseBase<EigenDerived>& config) const;
};
 
template<typename ScalarType>
template<typename RealType, class Derived>
__host__ __device__ void mBodyOperatorSpace<ScalarType>::action_impl(
    int& resBasisNum, Complex_t<RealType>& coeff, int opOrdinal, int basisNum,
    Eigen::DenseBase<Derived>& workSpace) const {
// debug_print(__PRETTY_FUNCTION__);
#ifndef NDEBUG
    if(workSpace.size() < this->sysSize())
        debug_printf("%s\n\tworkSpace.size()=%d\n", __PRETTY_FUNCTION__, (int)workSpace.size());
#endif
    assert(workSpace.size() >= this->sysSize());
    resBasisNum = basisNum;
    coeff       = {1.0, 0};
    int                 resLocBasisNum;
    Complex_t<RealType> resLocCoeff;
 
    int posConfNum = opOrdinal / m_opConfDim;
    int opConfNum  = opOrdinal % m_opConfDim;
    this->m_ActingSites.ordinalToPart(workSpace, posConfNum);
    int locOpNum, locBasisNum;
 
    int index = 0;
    int base  = 1;
    for(int pos = 0; pos != this->sysSize(); ++pos, base *= this->baseSpace().dimLoc()) {
        locBasisNum = this->baseSpace().locState(basisNum, pos);
        if(workSpace[pos] == 0) {
            this->locSpace().action(resLocBasisNum, resLocCoeff, 0, locBasisNum);
            coeff *= resLocCoeff;
            continue;
        }
        locOpNum = 1 + this->m_opConfig.digit(opConfNum, index++);
#ifndef NDEBUG
        if(locOpNum >= 4) {
            printf("Error: mBodyOperatorSpace::%s Out of bounds locOpNum = %d\n\n", __func__,
                   locOpNum);
    #ifndef __CUDA_ARCH__
            std::cout << "locOpNum = " << locOpNum << std::endl;
            this->m_opConfig.status();
            std::exit(EXIT_FAILURE);
    #endif
        }
#endif
 
        this->locSpace().action(resLocBasisNum, resLocCoeff, locOpNum, locBasisNum);
        if(abs(resLocCoeff) < 1.0e-8) {
            coeff = {0, 0};
            return;
        }
        coeff *= resLocCoeff;
        resBasisNum += (resLocBasisNum - locBasisNum) * base;
    }
    return;
}
 
template<typename ScalarType>
template<class EigenDerived>
__host__ __device__ int mBodyOperatorSpace<ScalarType>::translate_impl(
    int opOrdinal, int trans, Eigen::DenseBase<EigenDerived>& workSpace) const {
    int posConfNum = opOrdinal / m_opConfDim;
    int opConfNum  = opOrdinal % m_opConfDim;
#ifndef __CUDA_ARCH__
    workSpace.resize(this->sysSize());
#else
    assert(workSpace.size() >= this->sysSize());
#endif
#ifndef NDEBUG
    if(workSpace.size() < this->sysSize())
        debug_printf("%s\n\tworkSpace.size()=%d\n", __PRETTY_FUNCTION__, (int)workSpace.size());
#endif
    m_ActingSites.ordinalToPart(workSpace, posConfNum);
    int opTrans = 0;
    for(int l = 0; l != trans; ++l) opTrans += workSpace[l];
 
    posConfNum = m_ActingSites.translate(posConfNum, trans, workSpace);
    opConfNum  = this->m_opConfig.shiftDigits(opConfNum, opTrans);
    return opConfNum + posConfNum * m_opConfDim;
}
 
template<typename ScalarType>
__host__ __device__ Eigen::RowVectorXi mBodyOperatorSpace<ScalarType>::ordinalToConfig_impl(
    int basisNum) const {
    int posConfNum = basisNum / m_opConfDim;
    int opConfNum  = basisNum % m_opConfDim;
 
    Eigen::RowVectorXi config(this->sysSize());
    this->m_ActingSites.ordinalToPart(config, posConfNum);
 
    int index = 0;
    for(int pos = 0; pos != this->sysSize(); ++pos) {
        if(config[pos] == 0) continue;
        config[pos] = 1 + this->m_opConfig.digit(opConfNum, index++);
    }
    return config;
}
 
#ifdef GPU
    #include "ManyBodyOperatorSpaceBase.cuh"
#endif