OffDiagonalElements.cpp File Reference

Functions
int	main (int argc, char **argv)

Variables
namespace	filesystem = std::experimental::filesystem

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
  #define MAXDATA  100
  #define SEED 0
  #include "PBC_TI/Framework/Fragments/setVariablesForEnsemble.cpp"
 
  Integer_t dim_sub; //全ヒルベルト空間の次元
  Integer_t info=0, failed=0; // カウンタ
  double OpRange;
 
  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 OutMatElems, OutEnergy, OutOpEigV;
  //******************** (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<Integer_t>   NdataInShell(dim_max);
  std::vector<double>         eigenEnergy(dim_max);
  std::vector<double>               MCAverage(dim_max);
  std::vector<double>        OpEigenValue(dim_max);
  matrix<Complex_t>        h(dloc_h , dloc_h );
  matrix<Complex_t>      loc(dloc_op, dloc_op);
  matrix<Complex_t>  loc_tot(dim_max, dim_max);
#ifndef GPU
  matrix<Complex_t>    h_tot(dim_max, dim_max);
  #define dh   h
  #define dloc loc
  #define dh_tot   h_tot
  // #define dloc_tot loc_tot
  matrix<Complex_t> dloc_tot(dim_max, dim_max);
#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>       dh(dloc_h , dloc_h );
  matrix_gpu<Complex_t>     dloc(dloc_op, dloc_op);
  matrix_gpu<Complex_t>   dh_tot(LDT, dim_max);
  matrix_gpu<Complex_t> dloc_tot(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_h( h,   dloc_h,  -1);
    generateLocal_op(loc, 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_h( h,   dloc_h,  -1);
    generateLocal_op(loc, dloc_op, -1);
 
  #ifdef GPU
    magma_setmatrix(dloc_h,  dloc_h,  sizeof(Complex_t), &*h.begin(),   dloc_h,  dh.ptr(),   dloc_h,  GPUconf.queue());
    magma_setmatrix(dloc_op, dloc_op, sizeof(Complex_t), &*loc.begin(), dloc_op, dloc.ptr(), dloc_op, GPUconf.queue());
  #endif
    std::string outDirName(baseDirName);
    {
      std::stringstream buff;
      buff << "/RawData/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 = constructGlobal_h( dh_tot,   dh,   num_h,  Sector.at(n), GPUconf);
                  constructGlobal_op(dloc_tot, dloc, num_op, Sector.at(n), GPUconf);
        #ifdef GPU
          magma_queue_sync(queue);
        #endif
      }
      T_pre += getETtime() -temp_t;
 
      debug_print("# Calculating eigenstate matrix elements.");
      temp_t = getETtime();
      {
          info = EigenMatrixElements(eigenEnergy, dim_sub, dh_tot, dloc_tot, GPUconf);
            if(info != 0) {
              #include "PBC_TI/Framework/Fragments/writeFailedMatrixToFile.cpp"
              continue;
            }
        #ifdef GPU
          magma_getmatrix(dim_sub, dim_sub, sizeof(Complex_t), dloc_tot.ptr(), dloc_tot.LD(), &*loc_tot.begin(), loc_tot.LD(), GPUconf.queue());
        #else
          loc_tot.copy(dloc_tot);
        #endif
      }
      T_diag += getETtime() -temp_t;
 
#ifndef NDEBUG
      if(dim_sub < 200) {
          debug_print("# h");
          h.print(std::cout, 4, 4, "All", GPUconf);
          debug_print("# loc");
          loc.print(std::cout, 4, 4, "All", GPUconf);
          debug_print("# loc_tot");
          loc_tot.print(dim_sub, dim_sub);
          debug_print("# dloc_tot");
          dloc_tot.print(std::cout, dim_sub, dim_sub, "All", GPUconf);
      }
#endif
 
      debug_print("# Writing results to a file.");
      temp_t  = getETtime();
      double gE = eigenEnergy.at(0);
      double EnergyRange = eigenEnergy.at(dim_sub-1) - gE;
    #pragma omp parallel sections
    {
      #pragma omp section
      {
        std::stringstream buff("");
        buff << "/OffDiagonalElements" << PRECISION << "_N" << n << ".txt";
        std::string filename(outDirName); filename += buff.str();
        OutMatElems.open(filename); checkIsFileOpen(OutMatElems, filename);
        OutMatElems << std::right << std::scientific << std::setprecision(6);
        OutMatElems << "                              \n";
        // if(argc > Nargs) loc_tot.printAbs(OutMatElems, dim_sub, dim_sub, "Lower");
        // else loc_tot.print(OutMatElems, dim_sub, dim_sub, "Lower");
        for(Integer_t j = 0; j < dim_sub; ++j) for(Integer_t k = 0; k < j; ++k) {
          OutMatElems << (eigenEnergy.at(j) - gE) / EnergyRange << " "
                      << (eigenEnergy.at(k) - gE) / EnergyRange << " "
                      << eigenEnergy.at(j) << " "
                      << eigenEnergy.at(k) << " "
                      << real(fabs(loc_tot.at(j,k))) << "\n";
        }
      }
      #pragma omp section
      {
        info = EigenValues(OpEigenValue, dim_sub, dloc_tot);
            if(info != 0) {
            #include "PBC_TI/Framework/Fragments/writeFailedMatrixToFile.cpp"
          }
        OpRange = OpEigenValue[dim_sub-1] - OpEigenValue[0];
        std::stringstream buff("");
        buff << "/OperatorEigenValues" << PRECISION << "_N" << n << ".txt";
        std::string filename(outDirName); filename += buff.str();
        OutOpEigV.open(filename); checkIsFileOpen(OutOpEigV, filename);
        OutOpEigV << std::right << std::scientific << std::showpos << std::setprecision(6);;
        OutOpEigV << "# dim = " << dim_sub << "\n";
        for(int i = 0;i < OpEigenValue.size(); ++i) OutOpEigV << OpEigenValue[i] << "\n";
        OutOpEigV.close();
      }
      #pragma omp section
      {
        std::stringstream buff("");
        buff << "/EigenEnergies" << PRECISION << "_N" << n << ".txt";
        std::string filename(outDirName); filename += buff.str();
        OutEnergy.open(filename); checkIsFileOpen(OutEnergy, filename);
        OutEnergy << std::right << std::scientific << std::showpos << std::setprecision(6);
        OutEnergy << "# dim = " << dim_sub << "\n";
        for(int i = 0;i < eigenEnergy.size(); ++i) OutEnergy << eigenEnergy[i] << "\n";
        OutEnergy.close();
      }
    }
        OutMatElems.seekp(0, std::ios_base::beg);
        OutMatElems << "# OpRange= " << OpRange << std::endl;
        OutMatElems.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(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
    magma_finalize();
  #endif
  return 0;
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem