FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > > Class Template Reference

#include <ETHmeasure.hpp>

Collaboration diagram for FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >:

Public Member Functions
	FuncETHmeasure ()
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)
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)

Private Types
using	TotalSpace = TotalSpace_
using	Scalar = Scalar_
using	Real = typename SubSpace< TotalSpace_, Scalar_ >::Real
using	Vector = Eigen::VectorXd

Private Attributes
std::vector< Vector >	m_expVal
std::vector< Vector >	m_mcAverage

Member Typedef Documentation

Real

template<class TotalSpace_ , typename Scalar_ >

using FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::Real = typename SubSpace<TotalSpace_, Scalar_>::Real

private

Scalar

template<class TotalSpace_ , typename Scalar_ >

using FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::Scalar = Scalar_

private

TotalSpace

template<class TotalSpace_ , typename Scalar_ >

using FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::TotalSpace = TotalSpace_

private

Vector

template<class TotalSpace_ , typename Scalar_ >

using FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::Vector = Eigen::VectorXd

private

Constructor & Destructor Documentation

FuncETHmeasure()

template<class TotalSpace_ , typename Scalar_ >

FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::FuncETHmeasure ( )

inline

                         : m_expVal(omp_get_max_threads()), m_mcAverage(omp_get_max_threads()) {
            debug_constructor_printf(1);
        }

Member Function Documentation

operator()() [1/2]

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...");
}

operator()() [2/2]

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();
}

Member Data Documentation

m_expVal

template<class TotalSpace_ , typename Scalar_ >

std::vector<Vector> FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::m_expVal

private

m_mcAverage

template<class TotalSpace_ , typename Scalar_ >

std::vector<Vector> FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >::m_mcAverage

private

---

The documentation for this class was generated from the following file:

• /Users/shoki/GitHub/Locality/Headers/StatMech/ETHmeasure.hpp