Fluc_randHOp.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 USAGE_opt "(Opt1).dE (Opt2).repMin (Opt3).repMax "
  #define Nargs_opt 3
  #include "PBC_TI/Framework/Fragments/setVariablesForEnsemble.cpp"
  #include "PBC_TI/Framework/Fragments/setVariablesForMCAverage.cpp"
 
  Integer_t dim_sub; //全ヒルベルト空間の次元
  Integer_t ndata, info, failed=0; // カウンタ
  double gE, EnergyRange, OpRange, OpEigMin, OpEigMax, s;
  double Fluc, MaxFluc, MaxPos;
 
  if( num_h!=2 || num_op !=2 ) {
    std::cerr << "Error: Currently, num_h and num_op must be two." << std::endl;
    std::exit(EX_USAGE);
  }
  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.");
  #include "PBC_TI/Framework/Fragments/MeasureOfETH_CreateOutputFiles.cpp"
  std::ofstream OutSample;
  //******************** (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"
  for(Integer_t n = n_max;n >= n_min; --n) {
    OutMeasure[n] << "# (n=" << n << ") dim=" << Sector[n].dim() << "\n"
                  << "# 1.(Sample No.) 2.(Fluc_2) 3.(Fluc_infty) 4.(Pos_max_flux) 5.(ndata)\n"
                  << std::endl << std::scientific;
  }
  //******************** (END)Translation invariance ********************
 
  //******************** SU(2) structure ********************
  debug_print("# Calculating multiplicity of irreducible representations of SU(2) algebra.");
  std::vector<std::vector<Integer_t>> multiplicity(n_max+1);
  {
    Integer_t numUpSpins;
    std::vector<Integer_t> numBasisWithFixedM, spin(n_max);
    for (Integer_t n = n_max;n >= n_min; --n) {
      if(Sector[n].momentum() == 0) dim_sub = Sector[n].dim() - Sector[n].parityOp().numPairs;
      else dim_sub = Sector[n].dim();
 
    spin.resize(n);
    numBasisWithFixedM.resize(n+1);
    std::fill(numBasisWithFixedM.begin(), numBasisWithFixedM.end(), 0);
    for(Integer_t j = 0;j < dim_sub; ++j) {
      i_to_spin(Sector[n].representative(j), spin, 2, n);
      numUpSpins = std::accumulate(spin.begin(), spin.end(), 0);
      numBasisWithFixedM.at(numUpSpins) += 1;
    }
 
    multiplicity[n].resize(n/2+1);
    multiplicity[n][n/2] = numBasisWithFixedM[n];
           for(Integer_t j = n-1;j >= (n+1)/2; --j) {
             multiplicity[n][ j-(n+1)/2 ] = numBasisWithFixedM[j] - numBasisWithFixedM[j+1];
           }
 
    print(numBasisWithFixedM, numBasisWithFixedM.size());
    print(multiplicity[n], multiplicity[n].size());
  }
  }
  //******************** (END) SU(2) structure ********************
 
  //********** Allocate CPU memories **********//
  debug_print("# Allocating CPU memories.");
  Integer_t const dim_max = n_max+1;
  std::vector<Integer_t>   NdataInShell(dim_max);
  std::vector<double>         eigenEnergy(dim_max);
  std::vector<double>            EXPvalue(dim_max);
  std::vector<double>               MCAverage(dim_max);
  std::vector<std::tuple<double,double,Integer_t>> Data;
  if(n_max >= 3) Data.reserve(n_max*n_max); else Data.reserve(4);
  matrix<Complex_t>        h(dloc_h , dloc_h );
  matrix<Complex_t>      loc(dloc_op, dloc_op);
  matrix<Complex_t>    h_tot(dim_max, dim_max);
  matrix<Complex_t>  loc_tot(dim_max, dim_max);
#ifndef GPU
  #define   dh_tot   h_tot
  #define dloc_tot loc_tot
#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 **********//
 
  #inlude "PBC_TI/Framework/Fragments/getAttributesOfMatrixElementsInSector.cpp"
#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 ********************
  //***************************************************************************
 
  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,  queue);
    magma_setmatrix(dloc_op, dloc_op, sizeof(Complex_t), &*loc.begin(), dloc_op, dloc.ptr(), dloc_op, queue);
  #endif
 
    if(repetition < 10000) {
      #include "PBC_TI/Framework/Fragments/createFilesForSampleData.cpp"
    }
    for(Integer_t n = n_max;n >= n_min; --n) {
      debug_print("# (rep,n)=(" << repetition << "," << n << ")");
      debug_print("# Constructing global matrices in the sector.");
      Data.resize(0); OpEigMin = NAN; OpEigMax = NAN;
      for(Integer_t Sint = n/2;Sint >= 0; --Sint) {
        s = n/2.0 - Sint;
        temp_t  = getETtime();
      #ifdef GPU
            // constructGlobal_h(   dh_tot, dh,   num_h,  Sector.at(n), conf);
            // constructGlobal_op(dloc_tot, dloc, num_op, Sector.at(n), conf);
            // magma_queue_sync(queue);
            dim_sub = constructGlobal_h(   h_tot, h,   s, Sector.at(n));
            constructGlobal_h( loc_tot, loc, s, Sector.at(n));
            magma_setmatrix(dim_sub, dim_sub, sizeof(Complex_t), &*h_tot.begin(), h_tot.LD(), dh_tot.ptr(), dh_tot.LD(), queue);
            magma_setmatrix(dim_sub, dim_sub, sizeof(Complex_t), &*loc_tot.begin(), loc_tot.LD(), dloc_tot.ptr(), dh_tot.LD(), queue);
      #else
            dim_sub = constructGlobal_h(   h_tot,   h, s, Sector.at(n));
            constructGlobal_h( loc_tot, loc, s, Sector.at(n));
      #endif
        T_pre += getETtime() -temp_t;
 
 
        debug_print("# Calculating eigenstate expectation values." << dim_sub << " " << dh_tot.LD());
        temp_t = getETtime();
            info = EigenExpValue(EXPvalue, eigenEnergy, dim_sub, dh_tot, dloc_tot, conf.queue());
              if(info != 0) {
                #include "PBC_TI/Framework/Fragments/writeFailedMatrixToFile.cpp"
                continue;
              }
            for(Integer_t j = 0;j < dim_sub; ++j) {
              Data.push_back( std::tie(eigenEnergy[j], EXPvalue[j], multiplicity[n][Sint]) );
            }
            debug_print("# Calculating the spectral range of the observable.");
            Diagonalize(dim_sub, dloc_tot, eigenEnergy);
            OpEigMin = std::fmin(OpEigMin, eigenEnergy[0]);
            OpEigMax = std::fmin(OpEigMax, eigenEnergy[dim_sub-1]);
        T_diag += getETtime() -temp_t;
      }
 
      if(Data.size() != (size_t)( (n/2+1)*( (n+1)/2+1 ) )) {
        std::cerr << "Error:: Data.size(" << Data.size() << ") != (n/2+1)*( (n+1)/2+1 ) = " << (n/2+1)*( (n+1)/2+1 ) << std::endl;
        std::exit(1);
      }
 
      std::sort(Data.begin(), Data.end());
      OpRange = OpEigMax - OpEigMin;
 
      debug_print("# Calculating measures of the ETH.");
      temp_t  = getETtime();
        EnergyRange = std::get<0>(Data[Data.size()-1]) - std::get<0>(Data[0]);
        gE = std::get<0>(Data[0]);
        for(Integer_t i = 0;i < Data.size(); ++i) {
          std::get<0>(Data[i]) = (std::get<0>(Data[i]) -gE)/EnergyRange;
          std::get<1>(Data[i]) = std::get<1>(Data[i])/OpRange;
        }
        MicroCanonicalAverage(MCAverage, NdataInShell, dE, Data);
        ndata = MeasureOfETH(Fluc, MaxFluc, MaxPos, Emin, Emax, MCAverage, Data);
      T_post += getETtime() -temp_t;
 
      OutMeasure[n] << std::setw(6)  << repetition << " " << std::showpos
                    << std::setw(13) << Fluc       << " "
                    << std::setw(13) << MaxFluc    << " "
                    << std::setw(13) << MaxPos     << " " << std::noshowpos << ndata << "\n";
      if(repetition < 10000) {
        for(Integer_t i = 0;i < dim_sub; ++i) {
          OutSample << std::setw(2)  << n << " " << std::showpos
                    << std::setw(13) << std::get<0>(Data[i])*EnergyRange+gE << " "
                    << std::setw(13) << std::get<0>(Data[i]) << " "
                    << std::setw(13) << std::get<1>(Data[i]) << "\n" << std::noshowpos;
        }
      }
    }
    if(repetition < 10000) OutSample.close();
    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;
      for(Integer_t n = n_min;n <= n_max; ++n) OutMeasure[n].flush();
    }
  }
 
  Finalize(argc, argv);
  return 0;
}

Variable Documentation

filesystem

namespace filesystem = std::experimental::filesystem