Overlap_vsOp.cpp File Reference

Functions
Integer_t	Bin (double x)
int	main (int argc, char **argv)

Variables
namespace	filesystem = std::experimental::filesystem

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"
 
  Integer_t dim_sub; //全ヒルベルト空間の次元
  Integer_t failed=0; // カウンタ
  double EnergyRange1;
  #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 OutMultiFractal, OutSample, OutParam;
  //******************** (END) Check for the directory structure ********************
 
  //***************************************************************************
  //******************** Allocation & Initialization **************************
  //***************************************************************************
#ifdef GPU
  magma_init();
  magma_queue_t queue = NULL;
  magma_int_t     dev = 0;
  magma_getdevice( &dev );
  magma_queue_create( dev, &queue );
#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>         Op1EigenValue(dim_max);
  std::vector<double>         Op2EigenValue(dim_max);
  matrix<Complex_t>  Op1(dloc_op, dloc_op);
  matrix<Complex_t>  Op2(dloc_op, dloc_op);
  matrix<Complex_t>  Overlap(dim_max, dim_max);
  std::vector<double> aveOverlap(NDIV);
#ifndef GPU
  matrix<Complex_t>  Op1_tot(dim_max, dim_max);
  matrix<Complex_t>  Op2_tot(dim_max, dim_max);
  #define dOp1 Op1
  #define dOp2 Op2
  #define dOp1_tot Op1_tot
  #define dOp2_tot Op2_tot
  #define dOverlap Overlap
#endif
  //********** (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>     dOp1(dloc_op, dloc_op);
  matrix_gpu<Complex_t>     dOp2(dloc_op, dloc_op);
  matrix_gpu<Complex_t> dOp1_tot(LDT, dim_max);
  matrix_gpu<Complex_t> dOp2_tot(LDT, dim_max);
  matrix_gpu<Complex_t> dOverlap(LDT, dim_max);
  //********** (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_op(Op1, dloc_op, -1);
    generateLocal_op(Op2, dloc_op, -1);
  }
  end = getETtime();
  std::cout << "(init_genrand): time=" << std::fixed << (end-start) << std::endl;
  //***************************************************************************
  //******************** (END) Allocation & Initialization ********************
  //***************************************************************************
// #ifndef NDEBUG
//   std::string   InFilename("/Users/Shoki/GitHub/temp/matrix.txt");
//   std::ifstream InFile(InFilename); checkIsFileOpen(InFile, InFilename);
//   debug_print("Input: " << InFilename);
// #endif
 
  start = getETtime();
  end   = start;
  for(Integer_t repetition = repMin;repetition <= repMax; ++repetition) {
    // 局所ハミルトニアンと局所物理量をランダムにとる ******************************//
    generateLocal_op(Op1, dloc_op, -1);
    generateLocal_op(Op2, dloc_op, -1);
 
    debug_print("# rep=" << repetition);
 
  #ifdef GPU
    magma_setmatrix(dloc_op, dloc_op, sizeof(Complex_t), &*Op1.begin(), dloc_op, dOp1.ptr(), dloc_op, queue);
    magma_setmatrix(dloc_op, dloc_op, sizeof(Complex_t), &*Op2.begin(), dloc_op, dOp2.ptr(), dloc_op, queue);
  #endif
 
    debug_print("# rep=" << repetition);
 
    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);
    }
 
    debug_print("# rep=" << repetition);
    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 = constructGlobal_op(dOp1_tot, dOp1, num_op, Sector.at(n), GPUconf);
                  constructGlobal_h( dOp2_tot, dOp2, num_op, Sector.at(n), GPUconf);
        #ifdef GPU
          magma_queue_sync(queue);
        #endif
      }
      T_pre += getETtime() -temp_t;
 
      debug_print("# Calculating overlaps of eigenvectors.");
      temp_t = getETtime();
      {
        if(Diagonalize(dim_sub, dOp1_tot, Op1EigenValue) != 0) {
          std::cerr << "# Warning: Failed in diagonalizing Op1 for (rep,n) = (" << repetition << "," << n << ")" << std::endl;
          failed += 1;
          continue;
        }
        if(Diagonalize(dim_sub, dOp2_tot, Op2EigenValue) != 0) {
          std::cerr << "# Warning: Failed in diagonalizing Op2 for (rep,n) = (" << repetition << "," << n << ")" << std::endl;
          failed += 1;
          continue;
        }
        matrixProduct_gemm(transC, transN, dim_sub, dim_sub, dim_sub, dOp1_tot, dOp2_tot, dOverlap, GPUconf);
      #ifdef GPU
        magma_getmatrix(dim_sub, dim_sub, sizeof(Complex_t), dOverlap.ptr(), dOverlap.LD(), &*Overlap.begin(), Overlap.LD(), queue);
      #endif
        EnergyRange1 = Op1EigenValue[dim_sub-1] - Op1EigenValue[0];
      }
      T_diag += getETtime() -temp_t;
 
      debug_print("# Writing results to a file.");
      temp_t  = getETtime();
        { // Calculate Multti-fractal dimensions
          std::stringstream buff("");
          buff << "/MultiFractal" << PRECISION << "_vsOp_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, Op2EigenValue, Overlap, dim_sub);
          OutMultiFractal.close();
        }
        {
          #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, Op1EigenValue, Op2EigenValue, Overlap);
          auto [mean, sigma, DOS] = OverlapSqVsEnergy(NDIV, dim_sub, Op1EigenValue, Op2EigenValue, Overlap);
 
          {
            std::stringstream buff(""); buff << "/OverlapSq" << PRECISION << "_vsOp_N" << n << ".txt";
            std::string filename(outDirName+buff.str());
            OutSample.open(filename); checkIsFileOpen(OutSample,filename);
            OutSample << "# 1.(Energy2) 2.(Energy1) 3.(Overlap^2) 4.(Num. data) 5.(Dos of Op2)\n" << std::endl
                      << std::scientific;
            for(size_t bin2 = 0;bin2 < NDIV; ++bin2) for(size_t bin1 = 0;bin1 < NDIV; ++bin1) {
              if( isnan(real(BinAverage.at(bin1,bin2))) ) continue;
              OutSample << std::setw(13) << (bin2+0.5)/(double)NDIV  << " "
                        << std::setw(13) << (bin1+0.5)/(double)NDIV  << " "
                        << std::setw(13) << real(BinAverage.at(bin1,bin2)) << " "
                        << std::setw(5)  << BinDataNum.at(bin1,bin2) << " "
                        << std::setw(13) << DOS.at(bin2)
                        << "\n" << std::noshowpos;
            }
            OutSample.close();
          }
          {
            std::stringstream buff(""); buff << "/Parameter" << PRECISION << "_vsOp_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 bin2 = 0;bin2 < NDIV; ++bin2) {
              if(DOS[bin2] == 0) continue;
              OutParam << std::setw(13) << (bin2+0.5)/(double)NDIV << " "
                       << std::setw(13) << mean[bin2]  << " "
                       << std::setw(13) << sigma[bin2] << " "
                       << std::setw(13) << EnergyRange1*sqrt(sigma[bin2]/DOS[bin2]/M_PI) << " "
                       << std::setw(13) << DOS[bin2]
                       << "\n" << std::noshowpos;
            }
            OutParam.close();
          }
        }
      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(3) << 100*T_pre /(end-start) << "%)"
                << ", T_diag="      << std::setprecision(6)  << std::setw(8)  << T_diag << "(" << std::setprecision(3) << 100*T_diag/(end-start) << "%)"
                << ", T_process="   << std::setprecision(6)  << std::setw(8)  << T_post << "(" << std::setprecision(3) << 100*T_post/(end-start) << "%)"
                << ", failed=" << failed
                << std::endl;
    }
  }
 
  Finalize(argc, argv);
  #ifdef GPU
    magma_finalize();
  #endif
  return 0;
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem