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/translation.hpp"
#include "Headers/mersenne_twister.hpp"
#include "Headers/mydebug.hpp"
#include <iostream>
#include <cassert>
#include <sysexits.h>
#include <cmath>
 
#define RANDTYPE argv[6] << argv[7] << "_"
 
constexpr char USAGE_sp[] = "";
constexpr int  Nargs_sp = 0;
 
matrix<double> prefactor;
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 void generateLocal_h(matrix<Complex_t<double>>& mat, Integer_t const d_lmat, Integer_t const seed){
  if(seed >= 0) {
    init_genrand(seed);
    std::cout << "# init_genrand(" << seed << ")\n";
  }
  // std::fill(mat.begin(), mat.end(), DoubleComplexZero);
  // mat.at(0,0) = sqrt(2)*Gaussian();
  // mat.at(0,1) = Complex_t<double>(Gaussian(),Gaussian());
  // mat.at(0,2) = mat.at(0,1)
  // mat.at(0,3) = Complex_t<double>(Gaussian(),Gaussian());
  // mat.at(1,1) = sqrt(2)*Gaussian(); // = (2,2)
  // mat.at(1,2) = Complex_t<double>(Gaussian(),Gaussian()); // = (2,1)
  // mat.at(1,3) = Complex_t<double>(Gaussian(),Gaussian()); // = (2.3)
  // mat.at(2,2) = mat.at(1,1);
  // mat.at(2,3) = mat.at(1,3);
  // mat.at(3,3) = sqrt(2)*Gaussian();
 
  if(d_lmat != 4) {
    std::cerr << "Error: Currently, only two-body operators on spin-1/2 chain is implemented." << std::endl;
    std::exit(1);
  }
  RandomMatrix_GUE(mat, d_lmat);
  // mat.at(0,0) = mat.at(0,0);
  mat.at(0,1) = (mat.at(0,1)+mat.at(0,2))/2.0;
  mat.at(0,2) = mat.at(0,1);
  // mat.at(0,3) = mat.at(0,3);
  mat.at(1,1) = (mat.at(1,1)+mat.at(2,2))/2.0;
  mat.at(1,2) = (mat.at(1,2)+mat.at(2,1))/2.0;
  mat.at(1,3) = (mat.at(1,3)+mat.at(2,3))/2.0;
  mat.at(2,2) = mat.at(1,1);
  mat.at(2,3) = mat.at(1,3);
  // mat.at(3,3) = mat.at(3,3);
 
  for(int i = 0;i < 4; ++i) for(int j = 0;j < i; ++j) {
    mat.at(i,j) = conj(mat.at(j,i));
  }
 
};
 
static inline void generateLocal_op(matrix<Complex_t<double>>& mat, Integer_t const d_lmat, Integer_t const seed){
  generateLocal_h(mat, d_lmat, seed);
};
 
static inline Integer_t constructGlobal_h(matrix<Complex_t<double>>& global, const matrix<Complex_t<double>>& local, double const s, TransSector const Sector) {
  Integer_t const imax = (Integer_t)(2*s);
  double    const n = Sector.N()/2.0;
  double m;
 
  std::fill(global.begin(), global.end(), DoubleComplexZero);
 
  // #pragma omp parallel for private(m,n,s)
  for(Integer_t i = 0;i <= imax; ++i) {
    m = i-s;
    if(i+2 <= imax) {
      global.at(i+2,i)  = sqrt((s-m)*(s-m-1)*(s+m+1)*(s+m+2)) *local.at(0,3); // 2nd degree
      global.at(i,i+2)  = conj(global.at(i+2,i));
    }
    if(i+1 <= imax) {
      global.at(i+1,i)  = 2*sqrt((s-m)*(s+m+1)) *( (n+m)*local.at(0,1) +(n-m-1)*local.at(1,3) ); // 2nd degree
      global.at(i,i+1)  = conj(global.at(i+1,i));
    }
    global.at(i,i)  = 2.0*(s*(s+1) -m*m -n) *local.at(1,2) +2.0*(n*n-m*m)*local.at(1,1); // 2nd degree
    global.at(i,i) += (n+m)*(n+m-1) *local.at(0,0) +(n-m)*(n-m-1) *local.at(3,3); // 2nd degree
  }
 
  #ifndef NDEBUG
    if(imax <= 6) {
      std::cout << "imax=" << imax << "\n";
      local.print(4,4);
      global.print(imax+1,imax+1);
      std::cout << "\n" << std::endl;
    }
  #endif
 
  return imax+1;
};
 
static inline Integer_t constructGlobal_op(matrix<Complex_t<double>>& global, const matrix<Complex_t<double>>& local, double const s, TransSector const Sector) {
  return constructGlobal_h(global, local, s, Sector);
};
 
#endif