Ensemble copy.hpp

Go to the documentation of this file.

#include "Headers/mytypes.hpp"
#include "Headers/MatrixUtils/MatrixUtils.hpp"
#include "Headers/RandomMatrices.hpp"
#include "Headers/HilbertSpace/HilbertSpace.hpp"
#include "unsupported/Eigen/KroneckerProduct"
#include <iostream>
 
#ifdef GPU
    #include "Headers/ObjectOnGPU.cuh"
    #include <cuda_runtime.h>
 
template<class Derived1, class Derived2, class TotalSpace, typename Scalar>
__global__ void constructHamiltonian_kernel(Eigen::MatrixBase<Derived1>*        dTotMat,
                                            Eigen::MatrixBase<Derived2> const*  dLocMat,
                                            SubSpace<TotalSpace, Scalar> const* subSpace) {
    int const idx = blockIdx.x * blockDim.x + threadIdx.x;
    int const idy = blockIdx.y * blockDim.y + threadIdx.y;
    if(idx >= (*dTotMat).rows() || idy >= (*dTotMat).cols()) return;
    if(idx < idy) return;
 
    int const dimLoc = (*dLocMat).cols();
 
    (*dTotMat)(idx, idy) = {0, 0};
    magma_index_t innerX, innerY, locIdX, locIdY;
    for(auto valuePosX = (*subSpace).basis().outerIndexPtr()[idx];
        valuePosX != (*subSpace).basis().outerIndexPtr()[idx + 1]; ++valuePosX) {
        innerX          = (*subSpace).basis().innerIndexPtr()[valuePosX];
        auto conjCoeffX = conj((*subSpace).basis().valuePtr()[valuePosX]);
 
        for(auto valuePosY = (*subSpace).basis().outerIndexPtr()[idy];
            valuePosY != (*subSpace).basis().outerIndexPtr()[idy + 1]; ++valuePosY) {
            innerY      = (*subSpace).basis().innerIndexPtr()[valuePosY];
            auto CoeffY = (*subSpace).basis().valuePtr()[valuePosY];
 
            // if((innerX - innerY) / dimLoc != 0) continue; // WRONG FORMULA used until 2020/10/07
            if(innerX / dimLoc != innerY / dimLoc) continue;
            locIdX = innerX % dimLoc;
            locIdY = innerY % dimLoc;
            (*dTotMat)(idx, idy) += conjCoeffX * (*dLocMat)(locIdX, locIdY) * CoeffY;
        }
    }
 
    (*dTotMat)(idy, idx) = conj((*dTotMat)(idx, idy));
    // printf("LDT=%d, id=(%d,%d), val=%+f%+f*i\t %s is yet to be tested.\n", LDT, idx, idy,
    //        real(dTotMat[idx + LDT * idy]), imag(dTotMat[idx + LDT * idy]), __func__);
    // printf("\t id=(%d,%d)\n", idx, idy);
    return;
}
#endif
 
template<typename Scalar_t = Complex_t<Real_t>>
class RandomMatrixEnsemble {
    private:
        HilbertSpace<int> const                     m_locHSpace;
        GUEgenerator<double>                        m_GUE;
        mutable Eigen::MatrixX< Complex_t<double> > m_locMat;
#ifndef GPU
        mutable Eigen::MatrixX<Scalar_t> m_dLocMat;
#else
        ObjectOnGPU< Eigen::MatrixX<Scalar_t> > m_dLocMat;
#endif
        mutable unsigned long m_count;
 
    public:
        RandomMatrixEnsemble(HilbertSpace<int> const& locHSpace)
            : m_locHSpace{locHSpace},
              m_GUE{},
              m_locMat(locHSpace.dim() * locHSpace.dim(), locHSpace.dim() * locHSpace.dim()),
              m_dLocMat(locHSpace.dim() * locHSpace.dim(), locHSpace.dim() * locHSpace.dim()),
              m_count{0} {
            debug_constructor_printf(1);
            debug_print("\tlocHSpace.dim()=" << locHSpace.dim()
                                             << ", m_dLocMat.ptr()=" << m_dLocMat.ptr());
        }
 
        void discardSamples(int num) {
            if(num <= 0) return;
            debug_print(__PRETTY_FUNCTION__);
            m_GUE.discard(m_locMat, num);
            m_count += num;
        };
        void getNewSample(void) {
            debug_print(__PRETTY_FUNCTION__);
            m_GUE(m_locMat);
            debug_print("Matrix on CPU:\n" << m_locMat);
#ifdef GPU
            m_dLocMat.set(m_locMat);
            debug_print("Matrix on GPU:\n" << m_dLocMat);
#endif
            ++m_count;
            return;
        };
 
            //      template<typename EigenMatrix, class TotalSpace,
            //               typename std::enable_if_t<!on_GPU_v<EigenMatrix>>* = nullptr>
            //      void constructHamiltonian(EigenMatrix&                          mat,
            //                                SubSpace<TotalSpace, Scalar_t> const& subSpace) const {
            //          static_assert(!on_GPU_v<EigenMatrix>,
            //                        "constructHamiltonian-CPU: Input matrix is on GPU.");
            //          debug_print(__func__ << ": Input matrix is NOT on GPU. subSpace.basis():\n"
            //                               << subSpace.basis());
            //          Eigen::SparseMatrix<Scalar_t> Identity(subSpace.totalSpace().dim() / m_locMat.rows(),
            //                                                 subSpace.totalSpace().dim() / m_locMat.rows());
            //          Identity.setIdentity();
            //          mat.resize(subSpace.dim(), subSpace.dim());
            //          debug_print("subSpace.basis().size = (" << subSpace.basis().rows() << ","
            //                                                  << subSpace.basis().cols() << ")");
            //          mat = subSpace.basis().adjoint() * Eigen::kroneckerProduct(m_locMat, Identity)
            //                * subSpace.basis();
            //      };
 
#ifdef GPU
        template<typename EigenMatrix, class TotalSpace, typename Scalar_,
                 typename std::enable_if_t<on_GPU_v<EigenMatrix>>* = nullptr>
        void constructHamiltonian(EigenMatrix&                         mat,
                                  SubSpace<TotalSpace, Scalar_> const& subSpace) const {
            debug_print(__func__ << ": Input matrix is on GPU. (Not tested)\n" << mat);
            debug_print(__func__ << ":\n\tsubSpace.dim()=" << subSpace.dim()
                                 << ", subSpace.totalSpace().dim()=" << subSpace.totalSpace().dim()
                                 << ", \n"
                                 << mat);
 
            mat.resize(subSpace.dim(), subSpace.dim());
 
            constexpr void (*kernel)(
                Eigen::MatrixBase< std::remove_reference_t<decltype(*mat.ptr())> >*,
                Eigen::MatrixBase< std::remove_reference_t<decltype(*m_dLocMat.ptr())> > const*,
                SubSpace<TotalSpace, Scalar_t> const*)
                = &constructHamiltonian_kernel;
            struct cudaFuncAttributes attr;
            cuCHECK(cudaFuncGetAttributes(&attr, kernel));
            int const nThread = min(subSpace.dim(), (int)sqrt(attr.maxThreadsPerBlock));
            int const nBlock  = (int)(subSpace.dim() % nThread == 0 ? subSpace.dim() / nThread
                                                                    : subSpace.dim() / nThread + 1);
            debug_print(__func__ << ":\t subSpace.dim()=" << subSpace.dim() << ", nBlock=" << nBlock
                                 << ", nThread=" << nThread);
 
            ObjectOnGPU< SubSpace<TotalSpace, Scalar_t> > dSubSpace(subSpace);
 
            constructHamiltonian_kernel<<<dim3(nBlock, nBlock, 1), dim3(nThread, nThread, 1)>>>(
                mat.ptr(), m_dLocMat.ptr(), dSubSpace.ptr());
            cuCHECK(cudaPeekAtLastError());
            cuCHECK(cudaDeviceSynchronize());
            debug_print(__func__ << ": mat.data() = " << mat.data());
            debug_print(__func__ << ": Result\n" << mat);
        };
#endif
 
        //      // Specific member functions for debugging
        //      Eigen::MatrixXd const& localHamiltonian() const { return m_locMat; }
 
        //      void status() const {
        //          std::cout << "Count:\t" << m_count << std::endl;
        //          m_GUE.status();
        //      }
};