Fluc_randHOp.cpp File Reference

Classes
struct	statm_t

Functions
void	read_off_memory_status (statm_t &result)
int	main (int argc, char **argv)

Variables
std::string const	USAGE_commom = "Usage: 1.This 2.n_max 3.n_min 4.num_h 5.num_op 6.division No. 7.Output dir "

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
    if(argc != N_ARGS) {
        std::cerr << USAGE_commom << USAGE_sp << std::endl;
        std::exit(EX_USAGE);
    }
 
    Integer_t const MAX     = 100000;  // 1つのファイルに書き込むデータ数の上限
    Integer_t const dim_loc = 2;       // サイト上のヒルベルト空間の次元
    Integer_t const div     = (Integer_t)std::atoll(argv[5]);
    Integer_t const num_h = (Integer_t)std::atoll(argv[3]);  // 局所ハミルトニアンに含まれるサイト数
    Integer_t const num_op = (Integer_t)std::atoll(argv[4]);  // 局所演算子に含まれるサイト数
    Integer_t const n_max   = (Integer_t)std::atoll(argv[1]);  //  全スピンの数
    Integer_t const n_min   = (Integer_t)std::atoll(argv[2]);
    Integer_t const dloc_h  = pow(dim_loc, num_h);  // 局所ハミルトニアンの空間の次元
    Integer_t const dloc_op = pow(dim_loc, num_op);  // 局所演算子の空間の次元
    Integer_t const seed    = 0;
    Integer_t const repMin  = MAX * div;
    Integer_t const repMax  = repMin + MAX - 1;
    std::cout << "# dloc_h=" << dloc_h << ", dloc_op=" << dloc_op << "\n";
 
    Integer_t    dim_sub;                  //全ヒルベルト空間の次元
    Integer_t    ndata, info, failed = 0;  // カウンタ
    double       gE, EnergyRange, OpRange, OpMin;
    double       Fluc, MaxFluc, MaxPos;  // ETHの指標
    double const cutoff = 0.45;
    double const dE = 0.02;  // ミクロカノニカル平均を計算するエネルギーシェルの幅
    double const Emin = 0 + cutoff;
    double const Emax = 1 - cutoff;
 
    std::stringstream buff;
 
    if(n_min > n_max) {
        std::cerr << "Error: n_min(" << n_min << ") > n_max(" << n_max << ")." << std::endl;
        std::exit(EX_USAGE);
    }
    if(!Initialize(argc, argv)) {
        std::cerr << "Error: Initialization failed." << std::endl;
        std::exit(EX_USAGE);
    }
#ifndef NDEBUG
    statm_t Usage;
    std::cerr << "# Successfully initialized." << std::endl;
#endif
 
    //******************** Check for the directory structure ********************
#ifndef NDEBUG
    std::cerr << "# Checking for the directory structure." << std::endl;
#endif
    std::string d_str1(argv[6]);
    d_str1 += "/LocalRandomMatrix/";
    d_str1 += QUOTE(ENSEMBLE);
    std::string d_str2(d_str1);
    buff.str("");
    buff.clear(std::stringstream::goodbit);
    buff << "/" << TYPE << "randHOp" << PRECISION << "_"
         << "Local" << num_h << "_Cut" << (Integer_t)(100 * cutoff) << "_Seed" << seed << "_Min"
         << n_min << "_Max" << n_max;
    d_str1 += buff.str();
    filesystem::create_directories(d_str1);
 
    if(div == 0) {
        buff.str("");
        buff.clear(std::stringstream::goodbit);
        buff << "/" << TYPE << "Structure" << PRECISION << "_"
             << "Local" << num_h << "_Seed" << seed << "_Min" << n_min << "_Max" << n_max;
        d_str2 += buff.str();
        filesystem::create_directories(d_str2);
    }
 
    std::ofstream OutMeasure[n_max + 1];
    for(Integer_t n = n_min; n <= n_max; ++n) {
        buff.str("");
        buff.clear(std::stringstream::goodbit);
        buff << "/MeasureOfETH_N" << n << "_" << div << ".txt";
        std::string filename(d_str1);
        filename += buff.str();
        OutMeasure[n].open(filename);
        if(OutMeasure[n].fail()) {
            std::cerr << "Can't open a file" << filename << "." << std::endl;
            std::exit(EX_CANTCREAT);
        }
        OutMeasure[n] << std::right;
    }
 
    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 ********************
#ifndef NDEBUG
    std::cerr << "# Calculating translation-invariant sectors." << std::endl;
#endif
    Integer_t const          wave_number = 0;
    std::vector<TransSector> Sector(n_max + 1);
    for(Integer_t n = n_max; n >= n_min; --n) {
        Sector[n].initialize(n, wave_number, dim_loc);
#ifdef GPU
        Sector[n].copyToGPU(queue);
#endif
        std::cout << "(n=" << n << ", k=" << wave_number << ") num_EqClass = " << Sector[n].dim()
                  << ", 2^n/n =" << pow(dim_loc, n) / n << "\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 ********************
 
    //********** Allocate CPU memories **********//
#if !defined(NDEBUG) && defined(MAGMA)
    std::cerr << "# Allocating CPU memories." << std::endl;
#endif
    Integer_t const        dim_max = Sector[n_max].dim();
    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);
    matrix<Complex_t>      h(dloc_h, dloc_h);
    matrix<Complex_t>      loc(dloc_op, dloc_op);
#ifndef GPU
    matrix<Complex_t> h_tot(dim_max, dim_max);
    matrix<Complex_t> loc_tot(dim_max, dim_max);
    #define dh_tot   h_tot
    #define dloc_tot loc_tot
#endif
    //********** (END) Allocate CPU memories **********//
 
#ifdef GPU
    //********** Allocate GPU memories **********//
    #ifndef NDEBUG
    std::cerr << "# Allocating GPU memories." << std::endl;
    #endif
    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 **********//
    #ifndef NDEBUG
    std::cerr << "# Determining GPU configuration." << std::endl;
    #endif
    GlobFuncPtr = MatrixElementsInSector;
    struct cudaFuncAttributes attr;
    cudaFuncGetAttributes(&attr, GlobFuncPtr);
    Integer_t nThread = (Integer_t)sqrt(attr.maxThreadsPerBlock);
    Integer_t nBlock  = (Integer_t)dim_max / nThread;
    if(dim_max > nBlock * nThread) nBlock += 1;
    GPUconfig conf(dim3(nBlock, nBlock, 1), dim3(nThread, nThread, 1), 0, queue);
 
    std::cout << "# constSizeBytes     = " << attr.constSizeBytes << "\n"
              << "# maxThreadsPerBlock = " << attr.maxThreadsPerBlock << "\n"
              << "# (dim_max=" << dim_max << ") nBlock=" << nBlock << ", nThread= " << nThread
              << "*" << nThread << " =" << nThread * nThread << std::endl;
    magma_queue_sync(queue);
    //********** (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 ********************
    //***************************************************************************
 
    start = getETtime();
    end   = start;
    for(Integer_t repetition = repMin; repetition <= repMax; ++repetition) {
        // 局所ハミルトニアンと局所物理量をランダムにとる ******************************//
#ifndef NDEBUG
        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);
        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) {
            //****************************** 局所演算子に関する性質を保存するファイルを用意 ******************************//
            buff.str("");
            buff.clear(std::stringstream::goodbit);
            buff << "/Sample_No" << repetition << ".txt";
            std::string filename(d_str2);
            filename += "/RawData";
            filesystem::create_directory(filename);
            filename += buff.str();
            OutSample.open(filename);
            if(OutSample.fail()) {
                std::cerr << "Can't open a file" << filename;
                std::exit(EX_CANTCREAT);
            }
            OutSample << "# dim_loc=" << dim_loc << ", num_h=" << num_h << ", num_op=" << num_op
                      << "\n"
                      << "# wave_number=" << wave_number << "\n"
                      << "# 1.(System size) 2.(Energy density) 3.(Relative Pos) 4.(Exp. value)\n"
                      << std::endl
                      << std::scientific;
            //****************************** (END) 局所演算子に関する性質を保存するファイルを用意 ******************************//
        }
        for(Integer_t n = n_max; n >= n_min; --n) {
#ifndef NDEBUG
            std::cout << "# (rep,n)=(" << repetition << "," << n << ")" << std::endl;
            std::cerr << "# Constructing global matrices in the sector." << std::endl;
#endif
            dim_sub = Sector[n].dim();
            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(conf.queue());
#else
            constructGlobal_h(h_tot, h, num_h, n, dim_loc);
            constructGlobal_h(loc_tot, loc, num_op, n, dim_loc);
#endif
            T_pre += getETtime() - temp_t;
 
#ifndef NDEBUG
            std::cerr << "# Calculating eigenstate expectation values." << std::endl;
#endif
            temp_t = getETtime();
            info   = EigenExpValue(EXPvalue, eigenEnergy, dim_sub, dh_tot, dloc_tot, conf.queue());
            if(info != 0) {
                failed += 1;
                buff.str("");
                buff.clear(std::stringstream::goodbit);
                buff << "/Failed_N" << n << "_No" << repetition << ".txt";
                std::string filename(d_str1);
                filename += buff.str();
                std::ofstream ErrFp;
                ErrFp.open(filename);
                if(!ErrFp.is_open()) {
                    std::cerr << "# Warning: Can't open a file " << filename << ". Continue..."
                              << std::endl;
                }
                else {
                    ErrFp << "# info=" << info
                          << "\n\n"
                             "# seed="
                          << seed << ", div=" << div << ", N=" << n << ", repetition=" << repetition
                          << "\n"
                             "# Local Hamiltonian\n";
                    h.print(ErrFp, dloc_h, dloc_h);
                    ErrFp << "\n";
                    ErrFp << "# Local operator\n";
                    loc.print(ErrFp, dloc_op, dloc_op);
                    ErrFp << "\n";
                    if(Sector[n].dim() < 110) {
                        ErrFp << "# Total Hamiltonian\n";
                        dh_tot.print(ErrFp, dim_sub, dim_sub, queue);
                        ErrFp << "# Total operator\n";
                        dloc_tot.print(ErrFp, dim_sub, dim_sub, queue);
                    }
                    ErrFp.close();
                }
                continue;
            }
#ifndef NDEBUG
            std::cerr << "# Calculating the spectral range of the observable." << std::endl;
#endif
            OpRange = SpectralRange(OpMin, dim_sub, dloc_tot);
            T_diag += getETtime() - temp_t;
 
#ifndef NDEBUG
            std::cerr << "# Calculating measures of the ETH." << std::endl;
#endif
            temp_t = getETtime();
            for(Integer_t i = 0; i < dim_sub; ++i) EXPvalue[i] /= OpRange;
            EnergyRange = eigenEnergy[dim_sub - 1] - eigenEnergy[0];
            gE          = eigenEnergy[0];
            for(Integer_t i = 0; i < dim_sub; ++i)
                eigenEnergy[i] = (eigenEnergy[i] - gE) / EnergyRange;
 
            MCAverage.resize(dim);
            NdataInShell.resize(dim);
            std::fill(MCAverage.begin(), MCAverage.end(), 0.0);
            std::fill(NdataInShell.begin(), NdataInShell.end(), 0.0);
            for(Integer_t i = 0; i < dim; ++i) {
                MCAverage[i] = MicroCanonicalAverage(NdataInShell[i], eigenEnergy[i], dE, dim_sub,
                                                     eigenEnergy, EXPvalue, i);
            }
 
            ndata = MeasureOfETH(Fluc, MaxFluc, MaxPos, Emin, Emax, dim_sub, eigenEnergy, EXPvalue,
                                 MCAverage);
            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)
                              << eigenEnergy[i] * EnergyRange + gE << " " << std::setw(13)
                              << eigenEnergy[i] << " " << std::setw(13) << EXPvalue[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;
}

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

USAGE_commom

std::string const USAGE_commom = "Usage: 1.This 2.n_max 3.n_min 4.num_h 5.num_op 6.division No. 7.Output dir "