TransSector.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 TransSector : public SubSpace< ManyBodySpinSpace, Complex_t<RealType> > {
    private:
        int m_momentum;
 
    public:
        using TotalSpace = ManyBodySpinSpace;
        __host__ __device__ TransSector(int k, ManyBodySpinSpace const& mbHSpace);
        __host__ __device__ TransSector(int k = 0, int systemSize = 0, int dimLoc = 0)
            : TransSector(k, ManyBodySpinSpace(systemSize, dimLoc)) {
            debug_constructor_printf((Default));
        }
 
        __host__ __device__ int momentum() const { return m_momentum; };
 
        __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__ TransSector<RealType>::TransSector(int const                   k,
                                                       ManyBodySpinSpace const& mbHSpace)
    : SubSpace< ManyBodySpinSpace, Complex_t<RealType> >{mbHSpace}, m_momentum(k) {
    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;
    for(int eqClass = 0; eqClass < this->totalSpace().transEqDim(); ++eqClass) {
        if((this->totalSpace().transPeriod(eqClass) * k) % L != 0) continue;
        dim += 1;
    }
 
#ifdef __CUDA_ARCH__
    this->basis().resize(this->totalSpace().dim(), dim, this->totalSpace().dim());
#else
    this->basis().resize(this->totalSpace().dim(), dim);
    this->basis().reserve(Eigen::VectorXi::Constant(dim, L));
#endif  // #ifdef __CUDA_ARCH__
 
    int numCompatible = 0, nonZeroId = 0;
    for(int eqClass = 0; eqClass != this->totalSpace().transEqDim(); ++eqClass) {
        if((this->totalSpace().transPeriod(eqClass) * k) % L != 0) continue;
#ifdef __CUDA_ARCH__
        assert(nonZeroId < this->basis().reserved());
        this->basis().outerIndexPtr()[numCompatible] = nonZeroId;
        for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
            this->basis().innerIndexPtr()[nonZeroId]
                = this->totalSpace().transEqClassRep(eqClass, trans);
            this->basis().valuePtr()[nonZeroId]
                = exp(-I * RealType(M_PI * 2 * k * trans / RealType(L)))
                  / RealType(sqrt(RealType(this->totalSpace().transPeriod(eqClass))));
            ++nonZeroId;
        }
#else
        for(int trans = 0; trans != this->totalSpace().transPeriod(eqClass); ++trans) {
            this->basis().insert(this->totalSpace().transEqClassRep(eqClass, trans), numCompatible)
                = exp(-I * RealType(M_PI * 2 * k * trans / RealType(L)))
                  / RealType(sqrt(RealType(this->totalSpace().transPeriod(eqClass))));
        }
#endif  // #ifdef __CUDA_ARCH__
        numCompatible += 1;
        assert(numCompatible <= dim);
    }
    assert(dim == numCompatible && "dim == numCompatible");
    this->basis().outerIndexPtr()[dim] = nonZeroId;
 
#ifdef __CUDA_ARCH__
    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().makeCompressed();
#endif  // #ifdef __CUDA_ARCH__
}
 
#ifdef GPU
    #include "TransSector.cuh"
#endif