ETHmeasure.hpp

Go to the documentation of this file.

#include "Headers/mydebug.hpp"
#include "Headers/MatrixUtils/MatrixUtils.hpp"
#include "Headers/HilbertSpace/HilbertSpace.hpp"
#include "Headers/HilbertSpace/OperatorSpace.hpp"
#include "Headers/StatMech/MicroCanonicalAverage.hpp"
#include <algorithm>
 
#if __has_include(<omp.h>)
    #include <omp.h>
#endif
 
#ifdef GPU
    #include "Headers/ObjectOnGPU.cuh"
    #include "Headers/MatrixUtils/MatrixUtils.cuh"
    #include "Headers/HilbertSpace/ManyBodyOperatorSpaceBase.cuh"
    #include "Headers/HilbertSpace/SubSpace.cuh"
 
__host__ __device__ static inline int round_down(int x, int base) { return (x / base) * base; }
__host__ __device__ static inline int round_up(int x, int base) {
    return (x % base == 0 ? x : (x / base + 1) * base);
}
__host__ __device__ static inline int round_nearest(int x, int base) {
    return (x % base < base / 2 ? (x / base) * base : (x / base + 1) * base);
}
__device__ unsigned dynamic_smem_size() {
    unsigned ret;
    asm volatile("mov.u32 %0, %dynamic_smem_size;" : "=r"(ret));
    return ret;
}
__device__ unsigned get_smid(void) {
    unsigned ret;
    asm("mov.u32 %0, %smid;" : "=r"(ret));
    return ret;
}
template<class Derived1, class Derived2, class Derived3, class Derived4, class Derived5,
         class TotalSpace, typename Scalar>
__global__ void ETHmeasure_kernel(
    Eigen::DenseBase<Derived1> const* __restrict__ resPtr,
    Eigen::DenseBase<Derived2> const* __restrict__ dEigValPtr,
    typename SubSpace<TotalSpace, Scalar>::Real dE,
    Eigen::DenseBase<Derived3> const* __restrict__ dEigVecPtr,
    SubSpace<TotalSpace, Scalar> const* __restrict__ subSpacePtr,
    SparseCompressed<Scalar> const* __restrict__ adjointBasisPtr,
    ManyBodyOperatorSpaceBase<Derived4> const* __restrict__ dmBodyOpSpacePtr, int const transEqDim,
    int const* __restrict__ transEqClassRep, int const* __restrict__ transPeriod,
    Eigen::DenseBase<Derived5>* __restrict__ dWorkPtr) {
    using Real = typename SubSpace<TotalSpace, Scalar>::Real;
 
    Eigen::DenseBase<Derived1>&       res     = *const_cast< Eigen::DenseBase<Derived1>* >(resPtr);
    Eigen::DenseBase<Derived2> const& dEigVal = *dEigValPtr;
    Eigen::DenseBase<Derived3> const& dEigVec = *dEigVecPtr;
    ManyBodyOperatorSpaceBase<Derived4> const& dmBodyOpSpace = *dmBodyOpSpacePtr;
    SparseCompressed<Scalar> const&            adjointBasis  = (*adjointBasisPtr);
    int const                                  dimSpace      = subSpacePtr->dim();
 
    int const opOrdinal = blockIdx.x + gridDim.x * blockIdx.y;
    if(opOrdinal >= transEqDim) return;
    #ifndef NDEBUG
    if(opOrdinal == 0 && threadIdx.x == 0)
        printf("%s\n\t transEqDim=%d, dimSpace=%d, dWorkPtr->rows()=%d\n", __PRETTY_FUNCTION__,
               int(transEqDim), int(dimSpace), int(dWorkPtr->rows()));
    #endif
 
    extern __shared__ int              sharedMem[];
    Real*                              expValPtr = reinterpret_cast<Real*>(sharedMem);
    Eigen::Map< Eigen::VectorX<Real> > expVal(expValPtr, dimSpace);
    Real*                              eigValPtr = expValPtr + expVal.size();
    int                                configStorage[14];
    Eigen::Map< Eigen::VectorX<int> >  config(&configStorage[0], dmBodyOpSpace.sysSize());
 
    for(int j = threadIdx.x; j < dimSpace; j += blockDim.x) { expVal[j] = 0; }
    __syncthreads();
 
    int*    outStatePtr = reinterpret_cast<int*>(expValPtr + expVal.size());
    Scalar* coeffPtr    = reinterpret_cast<Scalar*>(outStatePtr + blockDim.x);
    for(int j = 0; j != round_up(subSpacePtr->totalSpace().dim(), blockDim.x) / blockDim.x; ++j) {
        int const inState = threadIdx.x + j * blockDim.x;
        if(inState < subSpacePtr->totalSpace().dim()) {
            dmBodyOpSpace.action(outStatePtr[threadIdx.x], coeffPtr[threadIdx.x],
                                 transEqClassRep[opOrdinal], inState, config);
            coeffPtr[threadIdx.x] = adjointBasis.valuePtr()[outStatePtr[threadIdx.x]]
                                    * coeffPtr[threadIdx.x]
                                    * conj(adjointBasis.valuePtr()[inState]);
        }
        __syncthreads();
 
        for(int k = 0; k != blockDim.x; ++k) {
            int const inState = k + j * blockDim.x;
            if(inState >= subSpacePtr->totalSpace().dim()) continue;
            int const    id1   = adjointBasis.innerIndexPtr()[inState];
            int const    id2   = adjointBasis.innerIndexPtr()[outStatePtr[k]];
            Scalar const coeff = coeffPtr[k];
            for(int expvalId = threadIdx.x; expvalId < dimSpace; expvalId += blockDim.x) {
                atomicAdd(&expVal[expvalId],
                          real(conj(dEigVec(id2, expvalId)) * coeff * dEigVec(id1, expvalId)));
            }
        }
        __syncthreads();
    }
 
    int const nEigVals
        = (dynamic_smem_size() - sizeof(Real) * dimSpace) / (sizeof(Real) * dimSpace);
    for(int j = threadIdx.x; j < nEigVals * dimSpace; j += blockDim.x) {
        *(eigValPtr + j) = dEigVal(j % dimSpace);
    }
    __syncthreads();
    // Calculate squared differences from Microcanonical average;
    int const                          k = threadIdx.x % nEigVals;
    Eigen::Map< Eigen::VectorX<Real> > eigVal(eigValPtr + k * dimSpace, dimSpace);
    Real const                         opTransPeriod = transPeriod[opOrdinal];
    for(int j = threadIdx.x; j < dimSpace; j += blockDim.x) {
        int idMin = j, idMax = j;
        for(idMin = j; idMin >= 0 && eigVal(j) - eigVal(idMin) <= dE; --idMin) {};
        ++idMin;
        for(idMax = j; idMax < dimSpace && eigVal(idMax) - eigVal(j) <= dE; ++idMax) {};
        --idMax;
        assert(idMin >= 0);
        assert(idMax < dimSpace);
        assert(idMin <= idMax);
 
        Real MCave = 0;
        for(int k = idMin; k != idMax + 1; ++k) MCave += expVal[k];
        MCave /= static_cast<Real>(idMax - idMin + 1);
        Real const diff = opTransPeriod * (expVal[j] - MCave) * (expVal[j] - MCave);
        atomicAdd(&res(j, get_smid()), diff);
    }
}
#endif
 
#pragma omp declare reduction(+ : Eigen::VectorXd : omp_out=omp_out+omp_in) initializer(omp_priv = omp_orig)
 
template<class SubSpace_>
class FuncETHmeasure;
 
template<class TotalSpace_, typename Scalar_>
class FuncETHmeasure< SubSpace<TotalSpace_, Scalar_> > {
    private:
        using TotalSpace = TotalSpace_;
        using Scalar     = Scalar_;
        using Real       = typename SubSpace<TotalSpace_, Scalar_>::Real;
        using Vector     = Eigen::VectorXd;
        std::vector<Vector> m_expVal;
        std::vector<Vector> m_mcAverage;
 
    public:
        FuncETHmeasure() : m_expVal(omp_get_max_threads()), m_mcAverage(omp_get_max_threads()) {
            debug_constructor_printf(1);
        }
 
        template<class Derived>
        void operator()(Eigen::VectorXd& res, Eigen::MatrixBase<Derived> const& eigVector,
                        SubSpace<TotalSpace, Scalar> const& subSpace,
                        mBodyOperatorSpace<Scalar> const&   mBodyOpSpace,
                        MicroCanonicalAverage const&        MCaverage);
 
#ifdef GPU
    public:
        template<typename Matrix_t, class Derived>
        void operator()(Eigen::VectorXd& res, ObjectOnGPU<Matrix_t> const& dEigVector,
                        SubSpace<TotalSpace, Scalar> const&       subSpace,
                        ManyBodyOperatorSpaceBase<Derived> const& mBodyOpSpace,
                        MicroCanonicalAverage const&              MCaverage);
#endif
};
 
template<class TotalSpace_, typename Scalar_>
template<class Derived>
void FuncETHmeasure< SubSpace<TotalSpace_, Scalar_> >::operator()(
    Eigen::VectorXd& res, Eigen::MatrixBase<Derived> const& eigVector,
    SubSpace<TotalSpace, Scalar> const& subSpace, mBodyOperatorSpace<Scalar> const& mBodyOpSpace,
    MicroCanonicalAverage const& MCaverage) {
    debug_print("FuncETHmeasure " << __func__
                                  << ": eigVector is NOT on GPU. Using CPU algorithm...");
 
    mBodyOpSpace.computeTransEqClass();
    std::cout << "FuncETHmeasure():\tm = " << mBodyOpSpace.m()
              << ", \tmBodyOpSpace.dim()        = " << mBodyOpSpace.dim()
              << ", \tmBodyOpSpace.transEqDim() = " << mBodyOpSpace.transEqDim()
              << ", \tsubSpace.dim()            = " << subSpace.dim() << std::endl;
    res = Eigen::VectorXd::Zero(subSpace.dim());
    if(res.norm() > 1.0e-4) {
        std::cerr << "Error(" << __func__
                  << ") : failed to initialize res: res.norm() = " << res.norm() << " is too large."
                  << std::endl;
        std::exit(EXIT_FAILURE);
    }
    std::for_each(m_expVal.begin(), m_expVal.end(),
                  [&eigVector](auto& x) { x.resize(eigVector.cols()); });
 
    debug_print(eigVector);
    debug_print(subSpace.basis());
    debug_print(mBodyOpSpace.basisOp(0));
 
    omp_set_max_active_levels(1);
// #pragma omp parallel for reduction(+ : res)
    for(int opEqClass = 0; opEqClass < mBodyOpSpace.transEqDim(); ++opEqClass) {
        int opNum  = mBodyOpSpace.transEqClassRep(opEqClass);
        int thread = omp_get_thread_num();
 
        m_expVal[thread]
            = (eigVector.adjoint()
               * (subSpace.basis().adjoint() * mBodyOpSpace.basisOp(opNum) * subSpace.basis())
                     .pruned()
                     .eval()
               * eigVector)
                  .diagonal()
                  .real();
        debug_print("\tBefore MCaverage: opNum=" << opNum << ", thread=" << thread);
        MCaverage(m_mcAverage[thread], m_expVal[thread]);
 
        res += mBodyOpSpace.transPeriod(opEqClass)
               * (m_expVal[thread] - m_mcAverage[thread]).cwiseAbs2();
    }
    debug_print("FuncETHmeasure " << __func__
                                  << ": eigVector is NOT on GPU. Using CPU algorithm...");
}
 
#ifdef GPU
template<class TotalSpace_, typename Scalar_>
template<typename Matrix_t, class Derived>
void FuncETHmeasure< SubSpace<TotalSpace_, Scalar_> >::operator()(
    Eigen::VectorXd& res, ObjectOnGPU<Matrix_t> const& dEigVector,
    SubSpace<TotalSpace, Scalar> const&       subSpace,
    ManyBodyOperatorSpaceBase<Derived> const& mBodyOpSpace,
    MicroCanonicalAverage const&              MCaverage) {
 
        // dEigVector should be stored in row-major.
 
    debug_print("FuncETHmeasure " << __func__
                                  << ": dEigVector is on GPU. (Algorithm is NOT implemented)");
    int nGPUs;
    cuCHECK(cudaGetDeviceCount(&nGPUs));
    mBodyOpSpace.computeTransEqClass();
    std::cout << "FuncETHmeasure(): nGPUs = " << nGPUs
              << ", \tmBodyOpSpace.dim()        = " << mBodyOpSpace.dim()
              << ", \tmBodyOpSpace.transEqDim() = " << mBodyOpSpace.transEqDim()
              << ", \tsubSpace.dim()            = " << subSpace.dim() << std::endl;
    res = Eigen::VectorXd::Zero(subSpace.dim());
    if(res.norm() > 1.0e-4) {
        std::cerr << "Error(" << __func__
                  << ") : failed to initialize res: res.norm() = " << res.norm() << " is too large."
                  << std::endl;
        std::exit(EXIT_FAILURE);
    }
    size_t const expValMemSize  = sizeof(Real) * subSpace.dim();
    size_t const eigValMemSize  = sizeof(Real) * subSpace.dim();
    size_t const requiredSmSize = expValMemSize + eigValMemSize;
 
    // GPU-side preparation
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp, 0);
 
    ObjectOnGPU< Eigen::MatrixX<Real> > dRes(
        Eigen::MatrixX<Real>::Zero(subSpace.dim(), deviceProp.multiProcessorCount).eval());
    ObjectOnGPU< SubSpace<TotalSpace, Scalar> > dSubSpace(subSpace);
    ObjectOnGPU< SparseCompressed<Scalar> >     dAdjointBasis(subSpace.basis().adjoint());
    ObjectOnGPU< Eigen::VectorX<Real> >         dEigVal(MCaverage.eigVal());
    ObjectOnGPU< Eigen::MatrixX<int> >          dWork;
    ObjectOnGPU<Derived> dmBodyOpSpace(static_cast<Derived const&>(mBodyOpSpace));
    int*                 transEqClassRep = nullptr;
    int*                 transPeriod     = nullptr;
    cuCHECK(cudaMalloc(&transEqClassRep, mBodyOpSpace.transEqDim() * sizeof(int)));
    cuCHECK(cudaMalloc(&transPeriod, mBodyOpSpace.transEqDim() * sizeof(int)));
    cuCHECK(cudaMemcpyAsync(transEqClassRep, mBodyOpSpace.transEqClassRep().data(),
                            mBodyOpSpace.transEqDim() * sizeof(int), cudaMemcpyHostToDevice));
    cuCHECK(cudaMemcpyAsync(transPeriod, mBodyOpSpace.transPeriod().data(),
                            mBodyOpSpace.transEqDim() * sizeof(int), cudaMemcpyHostToDevice));
 
    void (*m_kernel)(
        Eigen::DenseBase< std::remove_reference_t<decltype(*dRes.ptr())> > const*,
        Eigen::DenseBase< std::remove_reference_t<decltype(*dEigVal.ptr())> > const*, Real,
        Eigen::DenseBase< std::remove_reference_t<decltype(*dEigVector.ptr())> > const*,
        SubSpace<TotalSpace, Scalar> const*, SparseCompressed< Scalar > const*,
        ManyBodyOperatorSpaceBase<Derived> const*, int const, int const*, int const*,
        Eigen::DenseBase< std::remove_reference_t<decltype(*dWork.ptr())> >*)
        = &ETHmeasure_kernel;
 
    // determine the configuration of shared memory
    int shared_memory_size = deviceProp.sharedMemPerMultiprocessor - 1024;
    int nEigVals           = (shared_memory_size - expValMemSize) / eigValMemSize;
    int smSize             = expValMemSize + nEigVals * eigValMemSize;
 
    cuCHECK(cudaFuncSetAttribute(m_kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, smSize));
    struct cudaFuncAttributes m_attr;
    cuCHECK(cudaFuncGetAttributes(&m_attr, m_kernel));
    shared_memory_size = m_attr.maxDynamicSharedSizeBytes;
 
    int constexpr warpSize = 32;
    int const nThread = min(round_up(subSpace.dim(), warpSize), m_attr.maxThreadsPerBlock);
    int const nBlock = static_cast<int>(sqrt(static_cast<double>(mBodyOpSpace.transEqDim()))) + 1;
    nEigVals = 2;
 
    smSize
        = expValMemSize + max(nEigVals * eigValMemSize, (sizeof(int) + sizeof(Scalar)) * nThread);
 
    std::cout << "\tnThread = " << nThread << ", nBlock = " << nBlock
              << ", m_attr.maxThreadsPerBlock = " << m_attr.maxThreadsPerBlock
              << ", requiredSmSize = " << requiredSmSize << ", smSize = " << smSize
              << ", shared_memory_size = " << shared_memory_size << ", nEigVals = " << nEigVals
              << ", deviceProp.sharedMemPerMultiprocessor = "
              << deviceProp.sharedMemPerMultiprocessor << std::endl;
    assert(nThread >= 1);
    assert(nBlock >= 1);
    assert(smSize <= shared_memory_size);
 
    m_kernel<<<dim3(nBlock, nBlock), dim3(nThread, 1), smSize>>>(
        dRes.ptr(), dEigVal.ptr(), static_cast<Real>(MCaverage.shellWidth()), dEigVector.ptr(),
        dSubSpace.ptr(), dAdjointBasis.ptr(), dmBodyOpSpace.ptr(), mBodyOpSpace.transEqDim(),
        transEqClassRep, transPeriod, dWork.ptr());
    cuCHECK(cudaGetLastError());
    cuCHECK(cudaFree(transEqClassRep));
    cuCHECK(cudaFree(transPeriod));
 
    cuCHECK(cudaDeviceSynchronize());
 
    res = dRes.get().template cast<double>().rowwise().sum();
}
#endif  // #ifdef GPU
 
// without the hack to avoid bank conflicts and SM-wise summation (SMid,opOrdinal), 80.8295 sec for L=11 on AI-g11.
// with the hack to avoid bank conflicts and without SM-wise summation (SMid,opOrdinal), 84.4658 sec for L=11 on AI-g11.
// with the hack to avoid bank conflicts and SM-wise summation (SMid,opOrdinal), 84.0022 sec for L=11 on AI-g11.
// without the hack to avoid bank conflicts and with SM-wise summation (SMid,opOrdinal), 81.1254 sec for L=11 on AI-g11.
// without the hack to avoid bank conflicts and with SM-wise summation (opOrdinal,SMid), 81.4284 sec for L=11 on AI-g11.
// with the hack to avoid bank conflicts and SM-wise summation (opOrdinal,SMid), 84.9689 sec for L=11 on AI-g11.