MaxETHMeasure_FewBodyObservables.cpp File Reference

Classes
class	IOManager

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

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
    constexpr size_t nargs = 7;
    constexpr size_t nopts = 2;
    size_t arg_count = nargs;
    if(argc < nargs + HamiltonianEnsemble::nargs()) {
        std::cerr << "Usage: 0.(This) 1.(LMax) 2.(LMin) 3.(hRepMin) 4.(hRepMax) 5.(shellWidth) "
                     "6.(OutDir)\n";
        std::cerr << "       " << HamiltonianEnsemble::usage(arg_count) << "\n";
        arg_count += HamiltonianEnsemble::nargs();
        // std::cerr << "       " << ObservableEnsemble::usage(arg_count) << "\n";
        std::cerr << "       " << "Opt1.(MMax) Opt2.(MMin)" << "\n";
        std::cerr << "argc = " << argc << std::endl;
        std::exit(EXIT_FAILURE);
    }
    size_t const LMax       = std::atoi(argv[1]);
    size_t const LMin       = std::atoi(argv[2]);
    size_t const hRepMin    = std::atoi(argv[3]);
    size_t const hRepMax    = std::atoi(argv[4]);
    double const shellWidth = std::atof(argv[5]);
    IOManager    data_writer(argv[6]);
    std::cout << "outDir = " << data_writer.outDir() << std::endl;
 
    arg_count += HamiltonianEnsemble::nargs();
    std::cout << "arg_count = " << arg_count << ", argc = " << argc << std::endl;
    size_t const MMax = (arg_count+0 < argc ? std::atoi(argv[arg_count+0]) : LMax);
    size_t const MMin = (arg_count+1 < argc ? std::atoi(argv[arg_count+1]) : 1);
    std::cout << "MMax=" << MMax << ", MMin=" << MMin << std::endl;
 
#ifdef MAGMA
    CHECK(magma_init());
#endif
 
    arg_count = nargs;
    auto&  hEnsemble = HamiltonianEnsemble::get_instance(argv + arg_count);
    arg_count += HamiltonianEnsemble::nargs();
    std::vector<std::vector<ObservableEnsemble>> obsEnsembles(LMax + 1);
    for(size_t L = LMin; L <= LMax; ++L) {
        obsEnsembles[L].resize(L + 1);
        for(size_t m = MMin; m <= std::min(MMax,L); ++m) { obsEnsembles[L][m] = ObservableEnsemble(L, m); }
    }
 
    FuncETHmeasure< SubSpace<ManyBodySpinSpace, Complex_t<Real_t>> > calculateETHmeasure;
 
    hEnsemble.discard(hRepMin);
 
    double start, t_int, end, temp_t;
    double T_pre = 0, T_diag = 0, T_meas = 0, T_IO = 0;
    start = getETtime();
    end   = start;
    for(size_t hRep = hRepMin; hRep <= hRepMax; ++hRep) {
        // Sample a Hamiltonian
        auto loc_hamiltonian = hEnsemble.sample();
 
        for(size_t L = LMin; L <= LMax; ++L) {
            std::cout << "hRep = " << hRep << ", L = " << L << std::endl;
            temp_t           = getETtime();
            auto hamiltonian = loc_hamiltonian.construct_globalOp(L);
            T_pre += getETtime() - temp_t;
 
            // Diagonalize a Hamiltonian to get eigenvalues and eigenvectors
            temp_t = getETtime();
            HamiltonianEnsemble::EigenSolver eigSolver(std::move(hamiltonian));
            T_diag += getETtime() - temp_t;
 
            auto eigRange
                = *std::max_element(eigSolver.eigenvalues().begin(), eigSolver.eigenvalues().end())
                  - *std::min_element(eigSolver.eigenvalues().begin(),
                                      eigSolver.eigenvalues().end());
            MicroCanonicalAverage const MCAverage(eigSolver.eigenvalues().template cast<double>(),
                                                  eigRange * shellWidth);
 
            Eigen::MatrixXd ETHmeasure
                = Eigen::MatrixXd::Zero(eigSolver.eigenvectors().cols(), L + 1);
            Eigen::VectorXd colVec;
 
            temp_t = getETtime();
            for(size_t m = MMin; m <= std::min(MMax,L); ++m) {
                auto& opEnsemble = obsEnsembles[L][m];
                auto& stateSpace = HamiltonianEnsemble::stateSpace(L);
                auto& opSpace    = opEnsemble.operatorSpace();
                std::cout << "m=" << m << ", opSpace.dim() = " << opSpace.dim() << std::endl;
                calculateETHmeasure(colVec, eigSolver.eigenvectors(), stateSpace, opSpace,
                                    MCAverage);
                ETHmeasure.col(m) += colVec;
            }
            T_meas += getETtime() - temp_t;
 
            double maxDiff;
            {  // Verifying the results
                auto theorySum = Eigen::VectorXd::Ones(MCAverage.dim())
                                 - MCAverage.dimShell().cast<double>().cwiseInverse();
                maxDiff = (ETHmeasure.rowwise().sum() - theorySum).maxCoeff();
                if(maxDiff > 1.0e-4) {
                    std::cout << "ETHmeasure.rowwise().sum():\n"
                              << ETHmeasure.rowwise().sum() << "\n"
                              << std::endl;
                    std::cout << "theorySum:\n" << theorySum << "\n" << std::endl;
                    std::cout << "ETHmeasure:\n" << ETHmeasure << "\n" << std::endl;
                    std::cerr << "Error : maxDiff = " << maxDiff << " is too large." << std::endl;
                    continue;
                    // std::exit(EXIT_SUCCESS);
                }
            }
 
            temp_t = getETtime();
            data_writer.save_ResultsToFile(eigSolver.eigenvalues(), ETHmeasure, hRep, L,
                                           MCAverage.dimShell(), maxDiff);
            T_IO += getETtime() - temp_t;
        }
 
        if(hRep % 10 == 9 || hRep == hRepMax) {
            t_int = end;
            end   = getETtime();
            std::cerr << "(total=" << std::setw(6) << hRep + 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_meas=" << std::setprecision(6) << std::setw(8) << T_meas << "("
                      << std::setprecision(1) << 100 * T_meas / (end - start) << "%)"
                      << ", T_IO=" << std::setprecision(6) << std::setw(8) << T_IO << "("
                      << std::setprecision(1) << 100 * T_IO / (end - start) << "%)" << std::endl;
        }
    }
 
#ifdef MAGMA
    CHECK(magma_finalize());
#endif
    return EXIT_SUCCESS;
}