generateRM.hpp

Go to the documentation of this file.

// Generate a random matrix used by all programs in this directory
#ifndef GENERATE_RM
#define GENERATE_RM
#define _USE_MATH_DEFINES
 
#include "Headers/mytypes.hpp"
#include "Headers/LocalRandomMatrix.hpp"
#include "Headers/matrix_util.hpp"
#include "Headers/spin_int.hpp"
#include "Headers/mersenne_twister.hpp"
#include "Headers/mydebug.hpp"
#include <iostream>
#include <cassert>
#include <sysexits.h>
#include <cmath>
 
#define RANDTYPE argv[Nargs_common+0] << argv[Nargs_common+1] << "H_" << argv[Nargs_common+2] << argv[Nargs_common+3] << "Op"
#define USAGE_sp "sp1.(Decay type of H) sp2.(Decay parameter for H) sp3.(Decay type of Op) sp4.(Decay parameter for Op) "
#define Nargs_sp 4
 
static inline double Prefactor(Integer_t const j, std::string const Type, double const param) {
       if(Type == "Exp")   return exp( -(j-1)*param );
  else if(Type == "Pow")   return pow(j,-param);
  else if(Type == "Short") return (j==1 ? 1: 0);
  else if(Type == "Full")  return 1;
  else if(Type == "CRMT")  return NAN;
  else {
  std::cerr << "Usage: Decay type (argument 7) must be either \"Exp\", \"Pow\", \"Short\", \"Full\", or \"CRMT\"." << std::endl;
  std::exit(EX_USAGE);
  }
}
 
std::vector<std::vector<double>> prefactorH;
std::vector<std::vector<double>> prefactorOp;
static inline bool Initialize(int const argc, char const*const* argv, int const Nargs_common) {
  Integer_t   const n_max = (Integer_t)std::atoll(argv[1]);
  Integer_t   const n_min = (Integer_t)std::atoll(argv[2]);
  std::string const          TypeH(argv[Nargs_common+0]);
  double      const paramH  = atof(argv[Nargs_common+1]);
  std::string const         TypeOp(argv[Nargs_common+2]);
  double      const paramOp = atof(argv[Nargs_common+3]);
 
  prefactorH.resize(n_max+1); prefactorOp.resize(n_max+1);
  for(Integer_t n = n_min;n <= n_max; ++n) {
    prefactorH[n].resize(n);  prefactorH[n][0] = 0;
    prefactorOp[n].resize(n); prefactorOp[n][0] = 0;
    for(Integer_t j = 1;j <= n/2; ++j) {
      prefactorH[n][j]  = Prefactor(j, TypeH,  paramH);
      prefactorOp[n][j] = Prefactor(j, TypeOp, paramOp);
      prefactorH[n].at(n-j)  = prefactorH[n][j];
      prefactorOp[n].at(n-j) = prefactorOp[n][j];
    }
    std::cout << "# Decay factor for Hamiltonians." << std::endl;
    print(prefactorH[n], n);
    std::cout << "# Decay factor for observables." << std::endl;
    print(prefactorOp[n], n);
  }
  return true;
}
static inline bool Finalize(int const argc, char const*const* argv) { return true; }
 
static inline Integer_t SectorDimension(TransSector const& Sector) {
  return (Sector.momentum()==0) ? Sector.dim()-Sector.parityOp().numPairs : Sector.dim();
}
 
static inline void generateLocal_h(matrix<Complex_t>& mat, Integer_t const d_lmat, Integer_t const seed){
  if(seed >= 0) {
    init_genrand(seed);
    std::cout << "# init_genrand(" << seed << ")\n";
  }
  RandomMatrix_GUE(mat, d_lmat);
};
static inline void generateLocal_op(matrix<Complex_t>& mat, Integer_t const d_lmat, Integer_t const seed){
  if(seed >= 0) init_genrand(seed);
  generateLocal_h(mat, d_lmat, seed);
};
 
static inline void MatrixElementsInSector(matrix<Complex_t>& matTot, matrix<Complex_t> const& matLoc, Integer_t const numLoc, TransSector const& Sector, std::vector<double> const& factor);
 
static inline Integer_t constructGlobalOperator(matrix<Complex_t>& global, matrix<Complex_t> const& local, Integer_t const numLoc, TransSector const& Sector, std::vector<double> const& factor, void* GPUconf = nullptr) {
  MatrixElementsInSector(global, local, numLoc, Sector, factor);
  return SectorDimension(Sector);
}
static inline Integer_t constructGlobal_h(matrix<Complex_t>& global, matrix<Complex_t> const& local, Integer_t const numLoc, TransSector const& Sector, void* GPUconf = nullptr){
  return constructGlobalOperator(global, local, numLoc, Sector, prefactorH[Sector.N()]);
 
  #ifndef NDEBUG
    if(Sector.N() <= 4) {
      Integer_t d_loch = pow(Sector.dloc(), numLoc);
      std::cout << "n=" << Sector.N() << "\n";
      local.print(d_loch,d_loch);
      global.print(Sector.dim(),Sector.dim());
      std::cout << "\n" << std::endl;
    }
  #endif
};
 
static inline Integer_t constructGlobal_op(matrix<Complex_t>& global, matrix<Complex_t> const& local, Integer_t const numLoc, TransSector const& Sector, void* GPUconf = nullptr){
  return constructGlobalOperator(global, local, numLoc, Sector, prefactorOp[Sector.N()]);
};
 
static inline Complex_t transEigenMatrixElements(Integer_t const idx, Integer_t const idy, matrix<Complex_t> const& matLoc, TransSector const& Sector, std::vector<double> const& factor) {
  Integer_t trans1, trans2;
  Integer_t Id1, Id2, locId1, locId2;
  Integer_t Ndiff, diff, pos;
  double theta;
 
  Complex_t res = ComplexZero<>;
  for(trans1 = 0;trans1 < Sector.periodicity(idx); ++trans1) {
            Id1 = transSpin(Sector.representative(idx),trans1,Sector.dloc(),Sector.N());
      for(trans2 = 0;trans2 < Sector.periodicity(idy); ++trans2) {
            Id2 = transSpin(Sector.representative(idy),trans2,Sector.dloc(),Sector.N());
          theta = 2.0*M_PI*Sector.momentum()*(trans2-trans1)/(Real_t)Sector.N();
          Ndiff = 0;
          for(pos = 1;pos<Sector.N() && Ndiff<2; ++pos) {
            if( Bit(Id1,pos,Sector.dloc(),Sector.N())!=Bit(Id2,pos,Sector.dloc(),Sector.N()) ) { Ndiff+=1; diff=pos; }
          }
          switch(Ndiff){
            case 0:
              for(diff = 1;diff < Sector.N(); ++diff) {
                locId1 = Bit(Id1,0,Sector.dloc(),Sector.N()) + Sector.dloc()*Bit(Id1,diff,Sector.dloc(),Sector.N());
                locId2 = Bit(Id2,0,Sector.dloc(),Sector.N()) + Sector.dloc()*Bit(Id2,diff,Sector.dloc(),Sector.N());
                res += factor.at(diff)*matLoc.at(locId1,locId2)*std::exp(I<>*theta);
              }
              break;
            case 1:
                locId1 = Bit(Id1,0,Sector.dloc(),Sector.N()) + Sector.dloc()*Bit(Id1,diff,Sector.dloc(),Sector.N());
                locId2 = Bit(Id2,0,Sector.dloc(),Sector.N()) + Sector.dloc()*Bit(Id2,diff,Sector.dloc(),Sector.N());
                res += factor.at(diff)*matLoc.at(locId1,locId2)*std::exp(I<>*theta);
              break;
            default:
              continue;
          }
      }
  }
 
  res /= sqrt( (double)Sector.periodicity(idx)*Sector.periodicity(idy) );
  return res;
}
 
static inline void MatrixElementsInSector(matrix<Complex_t>& matTot, matrix<Complex_t> const& matLoc, Integer_t const numLoc, TransSector const& Sector, std::vector<double> const& factor) {
  Integer_t const dim = (Sector.momentum() == 0) ? (Sector.dim()-Sector.parityOp().numPairs): Sector.dim();
  Integer_t const numParityEigens = (Sector.momentum() == 0) ? (Sector.dim()-2*Sector.parityOp().numPairs) : Sector.dim();
  debug_print("# MatrixElementsInSector: N=" << Sector.N() << ", dim=" << dim);
  std::fill(matTot.begin(), matTot.end(), ComplexZero<>);
 
  Complex_t tempMatElem1, tempMatElem2;
 
  #pragma omp parallel for schedule(dynamic, 1) private(tempMatElem1, tempMatElem2)
  for(int idx = 0;idx < numParityEigens; ++idx) {
    for(int idy = idx;idy < numParityEigens; ++idy) {
      matTot.at(idx,idy) = transEigenMatrixElements(idx, idy, matLoc, Sector, factor);
      matTot.at(idy,idx) = conj(matTot.at(idx,idy));
    }
    for(int idy = numParityEigens;idy < dim; ++idy) {
      tempMatElem1 = transEigenMatrixElements(idx, idy, matLoc, Sector, factor);
      tempMatElem2 = transEigenMatrixElements(idx, Sector.parityOp(idy), matLoc, Sector, factor);
      matTot.at(idx,idy) = (tempMatElem1 + tempMatElem2)/sqrt(2.0);
      matTot.at(idy,idx) = conj(matTot.at(idx,idy));
    }
  }
 
  #pragma omp parallel for schedule(dynamic, 1) private(tempMatElem1, tempMatElem2)
  for(int idx = numParityEigens;idx < dim; ++idx) {
    for(int idy = idx;idy < dim; ++idy) {
      tempMatElem1 = transEigenMatrixElements(idx, idy, matLoc, Sector, factor);
      tempMatElem2 = transEigenMatrixElements(idx, Sector.parityOp(idy), matLoc, Sector, factor);
      matTot.at(idx,idy) = tempMatElem1 + tempMatElem2;
      matTot.at(idy,idx) = conj(matTot.at(idx,idy));
    }
  }
 
  #pragma omp parallel for schedule(dynamic, 1) private(tempMatElem1, tempMatElem2)
  for(int idx = 0;idx < dim; ++idx) matTot.at(idx,idx) = ComplexOne<>*real(matTot.at(idx,idx));
}
 
#endif