ETHmeasure.hpp File Reference

Go to the source code of this file.

Classes
class	FuncETHmeasure< SubSpace< TotalSpace_, Scalar_ > >

Functions
__device__ unsigned	dynamic_smem_size ()
__device__ unsigned	get_smid (void)
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)

Function Documentation

dynamic_smem_size()

__device__ unsigned dynamic_smem_size

(

)

                                        {
    unsigned ret;
    asm volatile("mov.u32 %0, %dynamic_smem_size;" : "=r"(ret));
    return ret;
}

ETHmeasure_kernel()

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

get_smid()

__device__ unsigned get_smid

(

void

)

                                   {
    unsigned ret;
    asm("mov.u32 %0, %smid;" : "=r"(ret));
    return ret;
}