generateRM.cuh File Reference

Functions
__global__ void	MatrixElementsInSector (Complex_t *dmatTot, Integer_t const LDT, Complex_t const *__restrict__ dmatLoc, Integer_t const LDL, Integer_t const numLoc, Integer_t const n, Integer_t const dloc, Integer_t const momentum, Integer_t const SecDim, Integer_t const numPairs, Integer_t const *representatives, Integer_t const *periodicity, Real_t const *__restrict__ dfactor)
__device__ Complex_t< double >	MAGMA_CEXP (Real_t theta)
__device__ Complex_t	transEigenMatrixElements (Integer_t const idx, Integer_t const idy, Complex_t const *__restrict__ dmatLoc, Integer_t const LDL, Integer_t const n, Integer_t const dloc, Integer_t const momentum, Integer_t const *representatives, Integer_t const *periodicity, Real_t const *__restrict__ dfactor)

Variables
std::vector< Real_t * >	prefactorH
std::vector< Real_t * >	prefactorOp

Function Documentation

MAGMA_CEXP()

__device__ Complex_t< float > MAGMA_CEXP

(

Real_t

theta

)

                                                      {
  return MAGMA_Z_MAKE(cos(theta),sin(theta));
}

MatrixElementsInSector()

__global__ void MatrixElementsInSector	(	Complex_t *	dmatTot,
		Integer_t const	LDT,
		Complex_t const *__restrict__	dmatLoc,
		Integer_t const	LDL,
		Integer_t const	numLoc,
		Integer_t const	n,
		Integer_t const	dloc,
		Integer_t const	momentum,
		Integer_t const	SecDim,
		Integer_t const	numPairs,
		Integer_t const *	representatives,
		Integer_t const *	periodicity,
		Real_t const *__restrict__	dfactor
	)

                                                                                                                                                                                                                                                                                                                                                                                             {
  Integer_t const dim = (momentum==0) ? (SecDim - numPairs) : SecDim;
  Integer_t const numParityEigens = (momentum==0) ? (SecDim - 2*numPairs) : SecDim;
 
  int const idx = blockIdx.x*blockDim.x +threadIdx.x;
  int const idy = blockIdx.y*blockDim.y +threadIdx.y;
  if(idx>=dim || idy>=dim) return;
  if(idx < idy) {
    // dmatTot[idx+LDT*idy] = {NAN,NAN};
    return;
  }
 
  if(idy < numParityEigens) {
    if(idx < numParityEigens) {
      dmatTot[idx+LDT*idy] = transEigenMatrixElements(idx, idy, dmatLoc, LDL, n, dloc, momentum, representatives, periodicity, dfactor);
    } else {
      Complex_t tempMatElem1 = transEigenMatrixElements(idx, idy, dmatLoc, LDL, n, dloc, momentum, representatives, periodicity, dfactor);
      Complex_t tempMatElem2 = transEigenMatrixElements(idx+numPairs, idy, dmatLoc, LDL, n, dloc, momentum, representatives, periodicity, dfactor);
      dmatTot[idx+LDT*idy] = (tempMatElem1 + tempMatElem2)/sqrt(2.0);
    }
  } else {
    Complex_t tempMatElem1 = transEigenMatrixElements(idx, idy, dmatLoc, LDL, n, dloc, momentum, representatives, periodicity, dfactor);
    Complex_t tempMatElem2 = transEigenMatrixElements(idx+numPairs, idy, dmatLoc, LDL, n, dloc, momentum, representatives, periodicity, dfactor);
    dmatTot[idx+LDT*idy] = tempMatElem1 + tempMatElem2;
  }
 
  dmatTot[idy+LDT*idx] = conj(dmatTot[idx+LDT*idy]);
  if( isnan(real(conj(dmatTot[idy+LDT*idx])*dmatTot[idy+LDT*idx])) ) printf("(idx,idy) = (%d,%d)\n", idx, idy);
  if(idx==idy) dmatTot[idx+LDT*idx] = dComplexOne<>*real(dmatTot[idx+LDT*idx]);
 
  return;
}

transEigenMatrixElements()

__device__ Complex_t transEigenMatrixElements	(	Integer_t const	idx,
		Integer_t const	idy,
		Complex_t const *__restrict__	dmatLoc,
		Integer_t const	LDL,
		Integer_t const	n,
		Integer_t const	dloc,
		Integer_t const	momentum,
		Integer_t const *	representatives,
		Integer_t const *	periodicity,
		Real_t const *__restrict__	dfactor
	)

                                                                                                                                                                                                                                                                                                                           {
  Integer_t trans1, trans2;
  Integer_t Id1, Id2, locId1, locId2;
  Integer_t Ndiff, diff, pos;
  double theta;
 
  Complex_t res = dComplexZero<>;
  for(trans1 = 0;trans1 < periodicity[idx]; ++trans1) {
        Id1 = transSpin(representatives[idx],trans1,dloc,n);
    for(trans2 = 0;trans2 < periodicity[idy]; ++trans2) {
        Id2 = transSpin(representatives[idy],trans2,dloc,n);
      theta = 2.0*M_PI*momentum*(trans2-trans1)/(Real_t)n;
      Ndiff = 0;
      for(pos = 1;pos<n && Ndiff<2; ++pos) {
        if( Bit(Id1,pos,dloc,n)!=Bit(Id2,pos,dloc,n) ) { Ndiff+=1; diff=pos; }
      }
      switch(Ndiff){
        case 0:
          for(diff = 1;diff < n; ++diff) {
            locId1 = Bit(Id1,0,dloc,n) + dloc*Bit(Id1,diff,dloc,n);
            locId2 = Bit(Id2,0,dloc,n) + dloc*Bit(Id2,diff,dloc,n);
            if(locId1+LDL*locId2 >= 16) printf("# Case0: locId1+LDL*locId2 = %lld+%lld*%lld = %lld\n", locId1, LDL, locId2, locId1+LDL*locId2);
            res += dfactor[diff]*dmatLoc[locId1+LDL*locId2]*MAGMA_CEXP(theta);
          }
          break;
        case 1:
            locId1 = Bit(Id1,0,dloc,n) + dloc*Bit(Id1,diff,dloc,n);
            locId2 = Bit(Id2,0,dloc,n) + dloc*Bit(Id2,diff,dloc,n);
            if(locId1+LDL*locId2 >= 16) printf("# Case1: locId1+LDL*locId2 = %lld+%lld*%lld = %lld\n", locId1, LDL, locId2, locId1+LDL*locId2);
            res += dfactor[diff]*dmatLoc[locId1+LDL*locId2]*MAGMA_CEXP(theta);
          break;
        default:
          continue;
      }
    }
  }
  res /= sqrt( (Real_t)periodicity[idx]*periodicity[idy] );
 
  // if( isnan(real(conj(res)*res)) ) printf("(idx,idy) = (%lld,%lld)\n", idx, idy);
  return res;
}

Variable Documentation

prefactorH

std::vector<Real_t*> prefactorH

prefactorOp

std::vector<Real_t*> prefactorOp