MicroCanonicalAverage.hpp

Go to the documentation of this file.

 
#pragma once
 
#ifndef __NVCC__
    #define __host__
    #define __device__
#endif
#include "Eigen/Dense"
 
class MicroCanonicalAverage {
    private:
        template<class T>
        using Vector_t = Eigen::VectorX<T>;
 
        int              m_dim;
        double           m_shellWidth;
        Vector_t<double> m_eigVal;
        Vector_t<int>    m_idMin, m_idMax;
        Vector_t<int>    m_dimShell;
 
    public:
        template<typename InVector_t>
        __host__ __device__ MicroCanonicalAverage(InVector_t const& eigVal, double const shellWidth)
            : m_dim(eigVal.size()), m_shellWidth(shellWidth), m_eigVal{eigVal} {
            m_idMin.resize(m_dim);
            m_idMax.resize(m_dim);
            m_dimShell.resize(m_dim);
 
            // clang-format off
            #pragma omp parallel for
            // clang-format on
            for(auto j = 0; j < m_dim; ++j) {
                int idMin, idMax;
                for(idMin = j; idMin >= 0 && eigVal[j] - eigVal[idMin] <= shellWidth; --idMin) {};
                m_idMin[j] = (++idMin);
                for(idMax = j; idMax < m_dim && eigVal[idMax] - eigVal[j] <= shellWidth; ++idMax) {};
                m_idMax[j]    = (--idMax);
                m_dimShell[j] = idMax - idMin + 1;
            }
        }
 
        __host__ void status(void) const {
            std::cout << "MicroCanonicalAverage.status()"
                      << "\n"
                      << "\tDim:   " << this->dim() << "\n"
                      << "\tShellWidth: " << this->shellWidth() << std::endl;
        }
 
        __host__ __device__ int    dim() const { return m_dim; }
        __host__ __device__ double shellWidth() const { return m_shellWidth; }
 
        __host__ __device__ int dimShell(int index) const { return m_dimShell[index]; }
        __host__ __device__ Vector_t<int>& dimShell() { return m_dimShell; }
        __host__ __device__ Vector_t<int> const& dimShell() const { return m_dimShell; }
 
        __host__ __device__ Vector_t<double>& eigVal() { return m_eigVal; }
        __host__ __device__ Vector_t<double> const& eigVal() const { return m_eigVal; }
 
        template<class VectorType>
        __host__ __device__ void operator()(VectorType& res, VectorType const& expVal) const {
            res.resize(m_dim);
            // clang-format off
            #pragma omp parallel for
            // clang-format on
            for(int j = 0; j < m_dim; ++j) {
                res[j] = 0;
                for(auto id = m_idMin[j]; id <= m_idMax[j]; ++id) res[j] += expVal[id];
                res[j] /= double(m_dimShell[j]);
            }
        }
};