Overlap_vsH0.cpp File Reference

Classes
struct	statm_t

Functions
void	read_off_memory_status (statm_t &result)
Integer_t	Bin (double x)
int	main (int argc, char **argv)

Variables
namespace	filesystem = std::experimental::filesystem
statm_t	Usage

Function Documentation

Bin()

Integer_t Bin

(

double

x

)

                       {
  Integer_t res = (Integer_t)(x*NDIV);
  return (res==NDIV ? res-1: res);
};

main()

int main	(	int	argc,
		char **	argv
	)

                                {
  #ifndef MAX
  #  define MAX  100
  #endif
  #ifndef SEED
  #  define SEED 0
  #endif
  #include "PBC_TI/Framework/Fragments/setVariablesForEnsemble.cpp"
  std::string const Type(argv[Nargs_common+0]);
  double      const param = atof(argv[Nargs_common+1]);
 
  Integer_t dim_sub; //全ヒルベルト空間の次元
  double EnergyRange0, EnergyRange1;
  double ggE, gE1;
  #ifdef GPU
    constexpr magma_trans_t transC = MagmaConjTrans;
    constexpr magma_trans_t transN = MagmaNoTrans;
  #else
    char transC = 'C', transN = 'N';
  #endif
 
  if( !Initialize(argc, argv, Nargs_common) ) {
    std::cerr << "Error: Initialization failed." << std::endl;
    std::exit(EX_USAGE);
  }
  debug_print("# Successfully initialized.");
 
  //******************** Check for the directory structure ********************
  debug_print("# Checking for the directory structure.");
  std::ofstream OutSample, OutParam, OutMultiFractal;
  //******************** (END) Check for the directory structure ********************
 
  //***************************************************************************
  //******************** Allocation & Initialization **************************
  //***************************************************************************
#ifdef GPU
  magma_init();
  magma_queue_t queue = NULL, queue2 = NULL;
  magma_queue_t& queue1 = queue;
  magma_int_t     dev = 0;
  magma_getdevice( &dev );
  magma_queue_create( dev, &queue1 );
  magma_queue_create( dev, &queue2 );
#else
  void* GPUconf = nullptr;
#endif
 
  //******************** Translation invariance ********************
  debug_print("# Calculating translation-invariant sectors.");
  #include "PBC_TI/Framework/Fragments/TranslationInvariantSectors.cpp"
  //******************** (END)Translation invariance ********************
 
  //********** Allocate CPU memories **********//
  debug_print("# Allocating CPU memories.");
  Integer_t const dim_max = SectorDimension(Sector[n_max]);
  std::vector<double> eigenEnergy1(dim_max);
  std::vector<double> eigenEnergy0(dim_max);
  matrix<Complex_t>       h(dloc_h , dloc_h );
  matrix<Complex_t>  h1_tot(dim_max, dim_max);
  matrix<Complex_t>  h0_tot(dim_max, dim_max);
  matrix<Complex_t> Overlap(dim_max, dim_max);
  std::vector<double> aveOverlap(NDIV);
#ifndef GPU
  #define  dh h
  #define  dh1_tot h1_tot
  #define  dh0_tot h0_tot
  #define dOverlap Overlap
  #define dfactorUnity factorUnity
#endif
 
  constexpr double epsilon = 1.0e-4;
  std::vector<std::vector<Real_t>> factorUnity(n_max+1);
  for(Integer_t n = n_min;n <= n_max; ++n) {
    factorUnity[n].resize(n);
    factorUnity[n][0] = 0;
    for(Integer_t j = 1;j <= n/2; ++j) {
      factorUnity[n][j]      = (Real_t)std::pow(j,-epsilon);
      factorUnity[n].at(n-j) = factorUnity[n][j];
    }
    print(factorUnity[n],n);
  }
  //********** (END) Allocate CPU memories **********//
 
#ifdef GPU
  //********** Allocate GPU memories **********//
  debug_print("# Allocating GPU memories.");
  constexpr Integer_t GPU_UNIT = 32;
  Integer_t const LDT = magma_roundup(dim_max, GPU_UNIT);
  matrix_gpu<Complex_t>       dh(dloc_h, dloc_h);
  matrix_gpu<Complex_t>  dh1_tot(LDT, dim_max);
  matrix_gpu<Complex_t>  dh0_tot(LDT, dim_max);
  matrix_gpu<Complex_t> dOverlap(LDT, dim_max);
 
  std::vector<Real_t*> dfactorUnity(n_max+1);
  for(int n = n_min;n <= n_max; ++n) {
    CHECK( magma_malloc((void**)&dfactorUnity[n], n*sizeof(Real_t)) );
    CHECK( magma_setvector(n, sizeof(Real_t), &*factorUnity[n].begin(), 1, dfactorUnity[n], 1, queue) );
  }
  //********** (END) Allocate GPU memories **********//
 
  //********** Determine GPU configuration **********//
  #include "PBC_TI/Framework/Fragments/getAttributesOfMatrixElementsInSector.cpp"
  //********** (END) Determine GPU configuration **********//
#endif // #ifdef GPU
 
  double start, t_int, end, temp_t;
  double T_diag=0, T_post=0, T_pre=0;
  init_genrand(SEED);
  start = getETtime();
  for(Integer_t repetition = 0;repetition < repMin; ++repetition) {
    generateLocal_h(h, dloc_h, -1);
  }
  end = getETtime();
  std::cout << "(init_genrand): time=" << std::fixed << (end-start) << std::endl;
  //***************************************************************************
  //******************** (END) Allocation & Initialization ********************
  //***************************************************************************
 
  start = getETtime();
  end   = start;
 
  // for(Integer_t repetition = repMin;repetition <= repMin+9; ++repetition) {
  for(Integer_t repetition = repMin;repetition <= repMax; ++repetition) {
  #if !defined(NDEBUG) && defined(__linux__)
    read_off_memory_status(Usage);
        std::cerr << "# (GPU, rep= " << std::setw(6) << repetition
        << ") size="       << Usage.size
        << ", resident="   << Usage.resident
        << ", share="      << Usage.share
        << ", text="       << Usage.text
        << ", lib="        << Usage.lib
        << ", data="       << Usage.data
        << ", dt="         << Usage.dt
        << std::endl;
  #endif
  // 局所ハミルトニアンと局所物理量をランダムにとる ******************************//
    generateLocal_h(h, dloc_h, -1);
  #ifdef GPU
    magma_setmatrix(dloc_h, dloc_h, sizeof(Complex_t), &*h.begin(), dloc_h, dh.ptr(), dloc_h, queue);
  #endif
 
    std::string outDirName(baseDirName);
    {
      std::stringstream buff("");
      buff << "/ProcessedData/Sample_No" << repetition;
      outDirName += buff.str();
      outDirName = std::regex_replace(outDirName, std::regex("//"), "/");
      filesystem::create_directories(outDirName);
    }
 
    for(Integer_t n = n_max;n >= n_min; --n) {
      debug_print("# (rep,n)=(" << repetition << "," << n << ")");
      debug_print("# Constructing global matrices in the sector.");
      temp_t  = getETtime();
      {
        dim_sub = constructGlobalOperator(dh0_tot, dh, num_h, Sector.at(n), dfactorUnity[n], GPUconf);
        debug_print("# dim_sub=" << dim_sub);
        eigenEnergy0.resize(dim_sub);
        eigenEnergy1.resize(dim_sub);
 
        if(Type == "CRMT") {
          matrix<Complex_t>& Pert = Overlap;
          RandomMatrix_GUE(Pert, dim_sub);
          #ifdef GPU
            {
              matrix_gpu<Complex_t>& dPert = dOverlap;
              magma_setmatrix_async(dim_sub, dim_sub, sizeof(Complex_t), &*Pert.begin(), Pert.LD(), dPert.ptr(), dPert.LD(), queue);
              magma_copymatrix_async(dim_sub,dim_sub,sizeof(Complex_t),dh0_tot.ptr(),dh0_tot.LD(),dh1_tot.ptr(),dh1_tot.LD(),queue);
              magma_queue_sync(queue);
              Complex_t alpha = ComplexOne<>*epsilon;
              magma_axpy_wrapper(LDT*dim_max, alpha, dPert.ptr(), 1, dh0_tot.ptr(), 1, queue);
                if(Diagonalize(dim_sub, dh0_tot, eigenEnergy0) != 0) continue;
 
              EnergyRange0 = eigenEnergy0[dim_sub-1] - eigenEnergy0[0];
              alpha = ComplexOne<>*EnergyRange0*param/sqrt((double)dim_sub);
              magma_axpy_wrapper(LDT*dim_max, alpha, dPert.ptr(), 1, dh1_tot.ptr(), 1, queue);
            }
          #else
            {   Integer_t n = Pert.size(), one=1;
                Complex_t alpha = DoubleComplexOne*epsilon;
                h1_tot.copy(h0_tot);
                zaxpy_(&n, &alpha, &*Pert.begin(), &one, &*h0_tot.begin(), &one);
                  if(Diagonalize(dim_sub, dh0_tot, eigenEnergy0) != 0) continue;
 
                EnergyRange0 = eigenEnergy0[dim_sub-1] - eigenEnergy0[0];
                alpha = DoubleComplexOne*EnergyRange0*param/sqrt((double)dim_sub);
                zaxpy_(&n, &alpha, &*Pert.begin(), &one, &*h1_tot.begin(), &one);
            }
          #endif
        } else {
          constructGlobal_h(dh1_tot, dh, num_h, Sector.at(n), GPUconf);
          #ifdef GPU
            magma_queue_sync(queue);
          #endif
          if(Diagonalize(dim_sub, dh0_tot, eigenEnergy0) != 0) continue;
        }
      }
      T_pre += getETtime() -temp_t;
 
      debug_print("# Calculating overlaps of eigenvectors.");
      temp_t = getETtime();
      {
          if(Diagonalize(dim_sub, dh1_tot, eigenEnergy1) != 0) continue;
          matrixProduct_gemm(transC, transN, dim_sub, dim_sub, dim_sub, dh0_tot, dh1_tot, dOverlap, GPUconf);
      }
        debug_print("# EnergyRange0=" << EnergyRange0 << ", EnergyRange1=" << EnergyRange1 << ", eigenEnergy0.size()=" << eigenEnergy0.size() << ", eigenEnergy1.size()=" << eigenEnergy1.size());
        EnergyRange0 = eigenEnergy0.at(dim_sub-1) - eigenEnergy0[0];
        EnergyRange1 = eigenEnergy1.at(dim_sub-1) - eigenEnergy1[0];
        ggE = eigenEnergy0[0];
        gE1 = eigenEnergy1[0];
        #pragma omp parallel for
        for(size_t i = 0;i < dim_sub; ++i) {
          eigenEnergy0.at(i) = (eigenEnergy0.at(i)-ggE)/EnergyRange0;
          eigenEnergy1.at(i) = (eigenEnergy1.at(i)-gE1)/EnergyRange1;
        }
      #ifdef GPU
        magma_getmatrix(dim_sub, dim_sub, sizeof(Complex_t), dOverlap.ptr(), dOverlap.LD(), &*Overlap.begin(), Overlap.LD(), queue);
      #endif
      T_diag += getETtime() -temp_t;
 
      temp_t = getETtime();
      debug_print("# Preprocessing the overlaps.");
      { // Calculate Multi-fractal dimensions
        #ifdef GPU
          magma_getmatrix_async(dim_sub, dim_sub, sizeof(Complex_t), dh1_tot.ptr(), dh1_tot.LD(), &*h1_tot.begin(), h1_tot.LD(), queue1);
          if(Type == "Short") {
            magma_getmatrix_async(dim_sub, dim_sub, sizeof(Complex_t), dh0_tot.ptr(), dh0_tot.LD(), &*h0_tot.begin(), h0_tot.LD(), queue2);
          }
        #endif
        {
          std::stringstream buff("");
          buff << "/MultiFractal" << PRECISION << "_MF_N" << n << ".txt";
          std::string filename(outDirName); filename += buff.str();
          OutMultiFractal.open(filename); checkIsFileOpen(OutMultiFractal,filename);
          OutMultiFractal
                   << "# 1.(Energy) 2.(deg.2) 3.(deg.3) ...\n" << std::endl
                   << std::scientific;
          printMultiFractalDimensions(OutMultiFractal, dim_sub, 8, eigenEnergy1, Overlap, dim_sub);
          OutMultiFractal.close(); OutMultiFractal.clear();
        }
        {
          std::stringstream buff("");
          buff << "/MultiFractal" << PRECISION << "_H1vsFock_N" << n << ".txt";
          std::string filename(outDirName); filename += buff.str();
          OutMultiFractal.open(filename); checkIsFileOpen(OutMultiFractal,filename);
          OutMultiFractal
                   << "# 1.(Energy) 2.(deg.2) 3.(deg.3) ...\n" << std::endl
                   << std::scientific;
          #ifdef GPU
            magma_queue_sync(queue1);
          #endif
          printMultiFractalDimensions(OutMultiFractal, dim_sub, 8, eigenEnergy1, h1_tot, dim_sub);
          OutMultiFractal.close(); OutMultiFractal.clear();
        }
        if(Type == "Short") {
          std::stringstream buff("");
          buff << "/MultiFractal" << PRECISION << "_H0vsFock_N" << n << ".txt";
          std::string filename(outDirName); filename += buff.str();
          OutMultiFractal.open(filename); checkIsFileOpen(OutMultiFractal,filename);
          OutMultiFractal
                   << "# 1.(Energy) 2.(deg.2) 3.(deg.3) ...\n" << std::endl
                   << std::scientific;
          #ifdef GPU
            magma_queue_sync(queue2);
          #endif
          printMultiFractalDimensions(OutMultiFractal, dim_sub, 8, eigenEnergy0, h0_tot, dim_sub);
          OutMultiFractal.close();  OutMultiFractal.clear();
        }
      }
      { debug_print("# Calculating Energy dependence of overlaps.");
        #pragma omp parallel for
        for (size_t i = 0; i < dim_sub; ++i) for (size_t j = 0; j < dim_sub; ++j) {
          Overlap.at(i,j) = conj(Overlap.at(i,j))*Overlap.at(i,j);
        }
        auto [BinAverage, BinDataNum] = coarseGrainingMatrixByEnergy(NDIV, dim_sub, eigenEnergy0, eigenEnergy1, Overlap);
        auto [mean, sigma, DOS] = OverlapSqVsEnergy(NDIV, dim_sub, eigenEnergy0, eigenEnergy1, Overlap);
        {
          debug_print("# Writing results to a file.");
          std::stringstream buff(""); buff << "/Overlap" << PRECISION << "_MF_N" << n << ".txt";
          std::string filename(outDirName+buff.str());
          OutSample.open(filename); checkIsFileOpen(OutSample,filename);
          OutSample << "# 1.(Energy1) 2.(Energy0) 3.(Overlap^2) 4.(Num. data) 5.(Dos of H1)\n" << std::endl
                    << std::scientific;
          for(size_t bin1 = 0;bin1 < NDIV; ++bin1) for(size_t bin0 = 0;bin0 < NDIV; ++bin0) {
            if( isnan(real(BinAverage.at(bin0,bin1))) ) continue;
            OutSample << std::setw(13) << (bin1+0.5)/(double)NDIV  << " "
                      << std::setw(13) << (bin0+0.5)/(double)NDIV  << " "
                      << std::setw(13) << real(BinAverage.at(bin0,bin1)) << " "
                      << std::setw(5)  << BinDataNum.at(bin0,bin1) << " "
                      << std::setw(13) << DOS.at(bin1)
                      << "\n" << std::noshowpos;
          }
          OutSample.close(); OutSample.clear();
        }
        {
          std::stringstream buff(""); buff << "/Parameter" << PRECISION << "_MF_N" << n << ".txt";
          std::string filename(outDirName+buff.str());
          OutParam.open(filename); checkIsFileOpen(OutParam,filename);
          OutParam << "# 1.(Energy1) 2.(mean) 3.(width) 4.(Strength sigma) 5.(density of states)\n" << std::endl
                   << std::scientific;
          for(size_t bin1 = 0;bin1 < NDIV; ++bin1) {
            if(DOS[bin1] == 0) continue;
            OutParam << std::setw(13) << (bin1+0.5)/(double)NDIV << " "
                     << std::setw(13) << mean[bin1]  << " "
                     << std::setw(13) << sigma[bin1] << " "
                     << std::setw(13) << EnergyRange0*sqrt(sigma[bin1]/DOS[bin1]/M_PI) << " "
                     << std::setw(13) << DOS[bin1]
                     << "\n" << std::noshowpos;
          }
          OutParam.close(); OutParam.clear();
        }
      }
      T_post += getETtime() - temp_t;
    }
 
    if(repetition%10 == 9) {
      t_int = end; end = getETtime();
      std::cerr << "(total=" << std::setw(6) << repetition+1
                << "): timeINT="    << std::setprecision(6)  << std::setw(8)  << (end-t_int)
                << ", timeTOT="     << std::setprecision(6)  << std::setw(8)  << (end-start)
                << ", T_construct=" << std::setprecision(6)  << std::setw(10) << T_pre  << "(" << std::setprecision(1) << 100*T_pre /(end-start) << "%)"
                << ", T_diag="      << std::setprecision(6)  << std::setw(8)  << T_diag << "(" << std::setprecision(1) << 100*T_diag/(end-start) << "%)"
                << ", T_process="   << std::setprecision(6)  << std::setw(8)  << T_post << "(" << std::setprecision(1) << 100*T_post/(end-start) << "%)"
                << std::endl;
    }
  }
  Finalize(argc, argv);
  #ifdef GPU
    for(int n = n_min;n <= n_max; ++n) CHECK( magma_free(dfactorUnity[n]) );
    magma_queue_destroy(queue1);
    magma_queue_destroy(queue2);
    magma_finalize();
  #endif
  return 0;
}

read_off_memory_status()

void read_off_memory_status

(

statm_t &

result

)

  {
    const char* statm_path = "/proc/self/statm";
 
    FILE *f = fopen(statm_path,"r");
    if(!f){
      perror(statm_path);
      abort();
    }
    if(7 != fscanf(f,"%ld %ld %ld %ld %ld %ld %ld",
      &result.size,&result.resident,&result.share,&result.text,&result.lib,&result.data,&result.dt))
    {
      perror(statm_path);
      abort();
    }
    fclose(f);
  }

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem

Usage

statm_t Usage