ParityTransSector.hpp

Go to the documentation of this file.

#pragma once
 
#include "Headers/mytypes.hpp"
#include "HilbertSpace.hpp"
#include "Eigen/SparseCore"
#ifdef GPU
    #include <thrust/sort.h>
    #include <thrust/execution_policy.h>
#endif
 
template<typename RealType = Real_t>
class ParityTransSector : public SubSpace< ManyBodySpinSpace, Complex_t<RealType> > {
    private:
        int m_parity;
        int m_numParityPairs = 0;
 
    public:
        using TotalSpace = ManyBodySpinSpace;
        __host__ __device__ ParityTransSector(int parity, ManyBodySpinSpace const& mbHSpace);
        __host__ __device__ ParityTransSector(int parity = 0, int systemSize = 0, int dimLoc = 0)
            : ParityTransSector(parity, ManyBodySpinSpace(systemSize, dimLoc)) {
            debug_constructor_printf((Default));
        }
 
        __host__ __device__ constexpr int momentum() const { return 0; };
 
        __host__ __device__ constexpr int parity() const { return m_parity; };
 
        __host__ __device__ int rep(int j, int trans = 0) const {
            int innerId = this->basis().outerIndexPtr()[j] + (trans % this->period(j));
            assert(innerId < this->basis().nonZeros() && "innerId < this->basis().nonZeros()");
            return this->basis().innerIndexPtr()[innerId];
        };
 
        __host__ __device__ int period(int j) const {
            assert(j < this->dim() && "j < this->dim()");
            return this->basis().outerIndexPtr()[j + 1] - this->basis().outerIndexPtr()[j];
        };
};
 
template<typename RealType>
__host__ __device__
ParityTransSector<RealType>::ParityTransSector(int const parity, ManyBodySpinSpace const& mbHSpace)
    : SubSpace< ManyBodySpinSpace, Complex_t<RealType> >{mbHSpace},
      m_parity((parity >= 0 ? +1 : -1)) {
    debug_constructor_printf(1);
    constexpr Complex_t<RealType> I = Complex_t<RealType>(0, 1);
    if(this->totalSpace().sysSize() == 0) { return; }
 
    this->totalSpace().computeTransEqClass();
    debug_printf("%s:\n\t this->totalSpace().dim()=%d, this->totalSpace().transEqDim()=%d, "
                 "this->totalSpace().transPeriod().sum()=%d\n",
                 __PRETTY_FUNCTION__, this->totalSpace().dim(), this->totalSpace().transEqDim(),
                 this->totalSpace().transPeriod().sum());
    assert(this->totalSpace().transPeriod().sum() == this->totalSpace().dim());
 
    int L            = this->totalSpace().sysSize();
    int dim          = 0;
    m_numParityPairs = 0;
    Eigen::ArrayX<bool> calculated;
    calculated = Eigen::ArrayX<bool>::Constant(this->totalSpace().transEqDim(), false);
    for(int eqClass = 0; eqClass < this->totalSpace().transEqDim(); ++eqClass) {
        if(calculated(eqClass)) continue;
        calculated(eqClass)  = true;
        auto repState        = this->totalSpace().transEqClassRep(eqClass);
        auto revState        = this->totalSpace().reverse(repState);
        auto reversed        = this->totalSpace().stateToTransEqClass(revState);
        calculated(reversed) = true;
        // std::cout << "reversed = " << reversed << std::endl;
        // std::cout << this->totalSpace().ordinalToConfig(repState) << std::endl;
        // std::cout << this->totalSpace().ordinalToConfig(revState) << "\n" << std::endl;
        if(eqClass != reversed) { m_numParityPairs += 1; }
        dim += 1;
    }
    if(m_parity == -1) dim = m_numParityPairs;
    std::cout << "m_parity = " << m_parity << ", dim = " << dim
              << ", m_numParityPairs = " << m_numParityPairs << std::endl;
 
#ifdef __CUDA_ARCH__
    this->basis().resize(this->totalSpace().dim(), dim, this->totalSpace().dim());
 
    int nonZeroId = 0, classIdx = 0;
    if(m_parity != -1) {
        calculated = Eigen::ArrayX<bool>::Constant(this->totalSpace().transEqDim(), false);
        for(int eqClass = 0; eqClass != this->totalSpace().transEqDim(); ++eqClass) {
            if(calculated(eqClass)) continue;
            calculated(eqClass)  = true;
            auto repState        = this->totalSpace().transEqClassRep(eqClass);
            auto revState        = this->totalSpace().reverse(repState);
            auto reversed        = this->totalSpace().stateToTransEqClass(revState);
            calculated(reversed) = true;
            if(eqClass != reversed) continue;
 
            this->basis().outerIndexPtr()[classIdx] = nonZeroId;
            for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
                this->basis().innerIndexPtr()[nonZeroId]
                    = this->totalSpace().transEqClassRep(eqClass, trans);
                this->basis().valuePtr()[nonZeroId]
                    = 1.0 / RealType(sqrt(RealType(this->totalSpace().transPeriod(eqClass))));
                ++nonZeroId;
            }
            ++classIdx;
        }
    }
 
    calculated = Eigen::ArrayX<bool>::Constant(this->totalSpace().transEqDim(), false);
    for(int eqClass = 0; eqClass != this->totalSpace().transEqDim(); ++eqClass) {
        if(calculated(eqClass)) continue;
        calculated(eqClass)  = true;
        auto repState        = this->totalSpace().transEqClassRep(eqClass);
        auto revState        = this->totalSpace().reverse(repState);
        auto reversed        = this->totalSpace().stateToTransEqClass(revState);
        calculated(reversed) = true;
        if(eqClass == reversed) continue;
 
        this->basis().outerIndexPtr()[classIdx] = nonZeroId;
        for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
            this->basis().innerIndexPtr()[nonZeroId]
                = this->totalSpace().transEqClassRep(eqClass, trans);
            this->basis().valuePtr()[nonZeroId]
                = 1.0 / RealType(sqrt(RealType(2 * this->totalSpace().transPeriod(eqClass))));
            ++nonZeroId;
        }
        for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
            this->basis().innerIndexPtr()[nonZeroId]
                = this->totalSpace().transEqClassRep(reversed, trans);
            this->basis().valuePtr()[nonZeroId]
                = m_parity / RealType(sqrt(RealType(2 * this->totalSpace().transPeriod(reversed))));
            ++nonZeroId;
        }
        ++classIdx;
    }
    assert(dim == classIdx && "dim == classIdx");
    this->basis().outerIndexPtr()[dim] = nonZeroId;
    this->basis().setNonZeros(this->basis().outerIndexPtr()[dim]);
    for(int eqClass = 0; eqClass != dim; ++eqClass) {
        thrust::sort_by_key(
            thrust::seq, this->basis().innerIndexPtr() + this->basis().outerIndexPtr()[eqClass],
            this->basis().innerIndexPtr() + this->basis().outerIndexPtr()[eqClass + 1],
            this->basis().valuePtr() + this->basis().outerIndexPtr()[eqClass]);
    }
#else   // #ifdef __CUDA_ARCH__
    this->basis().resize(this->totalSpace().dim(), dim);
    this->basis().reserve(Eigen::VectorXi::Constant(dim, 2 * L));
 
    int classIdx = 0;
    if(m_parity != -1) {
        calculated = Eigen::ArrayX<bool>::Constant(this->totalSpace().transEqDim(), false);
        for(int eqClass = 0; eqClass != this->totalSpace().transEqDim(); ++eqClass) {
            if(calculated(eqClass)) continue;
            calculated(eqClass)  = true;
            auto repState        = this->totalSpace().transEqClassRep(eqClass);
            auto revState        = this->totalSpace().reverse(repState);
            auto reversed        = this->totalSpace().stateToTransEqClass(revState);
            calculated(reversed) = true;
            if(eqClass != reversed) continue;
 
            for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
                auto state = this->totalSpace().transEqClassRep(eqClass, trans);
                this->basis().insert(state, classIdx)
                    = 1.0 / RealType(sqrt(RealType(this->totalSpace().transPeriod(eqClass))));
            }
            ++classIdx;
        }
    }
 
    calculated = Eigen::ArrayX<bool>::Constant(this->totalSpace().transEqDim(), false);
    for(int eqClass = 0; eqClass != this->totalSpace().transEqDim(); ++eqClass) {
        if(calculated(eqClass)) continue;
        calculated(eqClass)  = true;
        auto repState        = this->totalSpace().transEqClassRep(eqClass);
        auto revState        = this->totalSpace().reverse(repState);
        auto reversed        = this->totalSpace().stateToTransEqClass(revState);
        calculated(reversed) = true;
        if(eqClass == reversed) continue;
 
        for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
            auto state = this->totalSpace().transEqClassRep(eqClass, trans);
            this->basis().insert(state, classIdx)
                = 1.0 / RealType(sqrt(RealType(2 * this->totalSpace().transPeriod(eqClass))));
        }
        for(int trans = 0; trans != this->totalSpace().transPeriod(reversed); ++trans) {
            auto state = this->totalSpace().transEqClassRep(reversed, trans);
            this->basis().insert(state, classIdx)
                = m_parity / RealType(sqrt(RealType(2 * this->totalSpace().transPeriod(reversed))));
        }
        ++classIdx;
    }
 
    this->basis().makeCompressed();
#endif  // #ifdef __CUDA_ARCH__
}
 
#ifdef GPU
    #include "ParityTransSector.cuh"
#endif