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 <iostream>
#include <sysexits.h>
#include <cmath>
 
static inline bool Initialize(int const argc, char const*const* argv, int const Nargs_common) { return true; }
static inline bool Finalize(int argc, char** argv)   { return true; }
 
static inline Integer_t SectorDimension(TransSector const& Sector) {
  return 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);
 
  mat.at(0,0) = -1.1217*ComplexOne<>;
  mat.at(1,0) = 0.0*ComplexOne<>; mat.at(1,1) = -0.4034*ComplexOne<>;
  mat.at(2,0) = 0.0*ComplexOne<>; mat.at(2,1) = 0.2730*ComplexOne<>; mat.at(2,2) = -0.4034*ComplexOne<>;
  mat.at(3,0) = -0.1222*ComplexOne<>; mat.at(3,1) = 0.0*ComplexOne<>; mat.at(3,2) = 0.0*ComplexOne<>; mat.at(3,3) = -1.1217*ComplexOne<>;
  for(int i = 0;i < 4; ++i) for(int j = i+1;j < 4; ++j) mat.at(i,j) = conj(mat.at(j,i));
  mat.print(d_lmat,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);
  RandomMatrix_GUE(mat, d_lmat);
};
 
static inline void MatrixElementsInSector(matrix<Complex_t>& matTot, matrix<Complex_t> const& matLoc, Integer_t const numLoc, TransSector const& 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){
  MatrixElementsInSector(global, local, numLoc, Sector);
  #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
  return Sector.dim();
};
 
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 constructGlobal_h(global, local, numLoc, Sector);
};
 
static inline void MatrixElementsInSector(matrix<Complex_t>& matTot, matrix<Complex_t> const& matLoc, Integer_t const numLoc, TransSector const& Sector) {
  std::fill(matTot.begin(), matTot.end(), ComplexZero<>);
 
  double theta;
  Integer_t Id1, Id2, locId1, locId2;
  Integer_t const dmatLoc = (Integer_t)pow(Sector.dloc(), numLoc);
 
  #pragma omp parallel for schedule(dynamic, 1) private(Id1, Id2, locId1, locId2, theta)
  for(int idx = 0;idx < Sector.dim(); ++idx) for(int idy = 0;idy <= idx; ++idy) {
    for(int trans1 = 0;trans1 < Sector.periodicity(idx); ++trans1) {
          Id1 = transSpin(Sector.representative(idx),trans1,Sector.dloc(),Sector.N());
       locId1 = Id1%dmatLoc;
      for(int trans2 = 0;trans2 < Sector.periodicity(idy); ++trans2) {
          Id2 = transSpin(Sector.representative(idy),trans2,Sector.dloc(),Sector.N());
       locId2 = Id2%dmatLoc;
        if(Id1/dmatLoc != Id2/dmatLoc) continue;
        theta = 2.0*M_PI*Sector.momentum()*(trans2-trans1)/(Real_t)Sector.N();
        matTot.at(idx,idy) += matLoc.at(locId1,locId2) *std::exp(I<>*theta);
      }
    }
    matTot.at(idx,idy) /= sqrt( (double)Sector.periodicity(idx)*Sector.periodicity(idy) );
    matTot.at(idy,idx)  = conj(matTot.at(idx,idy));
  }
}
 
#endif