RestrictionToMBodyOperators_MeanOp.cpp File Reference

Classes
class	MagmaController

Typedefs
using	SubHilbertSpace = TransSector< double >

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

Typedef Documentation

SubHilbertSpace

using SubHilbertSpace = TransSector<double>

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
    if(argc != 13) {
        std::cerr << "Usage: 1.(This) 2.(dimLoc) 3.(LMax) 4.(LMin) 4.(mMax) 5.(mMin) 6.(momentum) "
                     "7.(shellWidth) "
                     "8.(repMin) 9.(repMax) 10.(opRepMin) 11.(opRepMax) 12.(OutDir)\n"
                  << "argc = " << argc << std::endl;
        std::exit(EXIT_FAILURE);
    }
    int const    dimLoc     = std::atoi(argv[1]);
    int const    LMax       = std::atoi(argv[2]);
    int const    LMin       = std::atoi(argv[3]);
    int const    mMax       = std::atoi(argv[4]);
    int const    mMin       = std::atoi(argv[5]);
    int const    k          = std::atoi(argv[6]);
    double const shellWidth = std::atof(argv[7]);
    int const    repMin     = std::atoi(argv[8]);
    int const    repMax     = std::atoi(argv[9]);
    int const    opRepMin   = std::atoi(argv[10]);
    int const    opRepMax   = std::atoi(argv[11]);
    std::string  baseDirName(argv[12]);
    baseDirName += "/TypicalityOfETH/";
    baseDirName += QUOTE(ENSEMBLE);
    std::ofstream OutFs;
 
#ifdef MAGMA
    MagmaController::getInstance();
 
    // cuCHECK(cudaDeviceSetLimit(cudaLimitMallocHeapSize, size_t(4) * size_t(1024 * 1024 * 1024)));
    // Setting the heap size too large causes a "double free or corruption" error on some of the destructObject_kernel calls.
    std::size_t heapSize, stackSize;
    cuCHECK(cudaDeviceGetLimit(&stackSize, cudaLimitStackSize));
    cuCHECK(cudaDeviceGetLimit(&heapSize, cudaLimitMallocHeapSize));
    std::cout << "Max stack size = " << stackSize << " bytes \n"
              << "Max heap size = " << heapSize << " bytes \n"
              << std::endl;
#endif
 
    double                  EnergyRange, groundEnergy, diff;
    HilbertSpace<int> const H_loc(dimLoc);
    RandomMatrixEnsemble    ensemble(H_loc);  //GUEgenerator<>       GUE(0, 1);
                                              // Requirements for class RandomMatrixEnsemble
                                              // -method discardSamples(int)
                                              // -method getNewSample(void)
                                              // -method constructHamiltonian(int)
 
    int dimTot;
    {
        SubHilbertSpace const transSector(k, LMax, dimLoc);
        dimTot = transSector.dim();
    }
    std::cout << "\t dimTot = " << dimTot << "\n\n\n" << std::endl;
 
    Eigen::MatrixXd ETHmeasure;
    Eigen::VectorXd eigEnergy, theorySum, colVec;
 
    FuncETHmeasure< SubSpace<ManyBodySpinSpace, Complex_t<Real_t>> > calculateETHmeasure;
#ifdef GPU
    MatrixUtils::EigenSolver< ObjectOnGPU<
        Eigen::Matrix<Complex_t<Real_t>, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> >
        eigenSolver(dimTot);
#else
    MatrixUtils::EigenSolver< Eigen::MatrixXcd > eigenSolver(dimTot);
#endif
 
    debug_print("baseDirName = " << baseDirName);
 
    double start, t_int, end, temp_t;
    double T_diag = 0, T_post = 0, T_pre = 0, T_IO = 0;
 
    ensemble.discardSamples(repMin);  //GUE.discard(repMin - 1);  // discardSamples(repMin-1);
    debug_print(repMin << " samples are discarded.");
 
    start = getETtime();
    end   = start;
    for(auto repetition = repMin; repetition <= repMax; ++repetition) {
        temp_t = getETtime();
        std::string outDirName(baseDirName);
        {
            std::stringstream buff;
            buff << "/RawData/Hamiltonian_No" << repetition;
            outDirName += buff.str();
            outDirName = std::regex_replace(outDirName, std::regex("//"), "/");
            if(!std::filesystem::exists(outDirName)
               && !std::filesystem::create_directories(outDirName)) {
                std::cerr << "Error: Failed to create directories " << outDirName << std::endl;
                std::exit(EXIT_FAILURE);
            };
        }
        debug_print("repetition = " << repetition);
 
        ensemble.getNewSample();
        T_IO += getETtime() - temp_t;
 
        for(int L = LMax; L >= LMin; --L) {
            temp_t = getETtime();
            // { // Construct a sector of a many-body Hilbert space and a Hamiltonian within that sector
            ManyBodySpinSpace const H_tot(L, H_loc);
            SubHilbertSpace const   transSector(k, H_tot);
 
            eigenSolver.inMatrix().resize(transSector.dim(), transSector.dim());
            eigEnergy.resize(transSector.dim());
            ETHmeasure = Eigen::MatrixXd::Zero(transSector.dim(), L + 1).eval();
            diff       = ETHmeasure.norm();
            if(diff > 1.0e-4) {
                std::cerr << "Error : failed to initialize ETHmeasure: diff = " << diff
                          << " is too large." << std::endl;
                std::exit(EXIT_FAILURE);
            }
            debug_print("eigenSolver.inMatrix():\n" << eigenSolver.inMatrix());
            ensemble.constructHamiltonian(eigenSolver.inMatrix(), transSector);
            // }
            T_pre += getETtime() - temp_t;
 
            // Diagonalize the Hamiltonian to obtain eigenvalues and eigenvectors
            temp_t = getETtime();
            { eigenSolver.compute(); }
            T_diag += getETtime() - temp_t;
            // debug_print("eigenSolver.eigenvalues(): After diagonalization\n"
            //             << eigenSolver.eigenvalues());
            // debug_print("eigenSolver.eigenvectors(): After diagonalization\n"
            //             << eigenSolver.eigenvectors());
 
            // {
            eigEnergy    = eigenSolver.eigenvalues();
            EnergyRange  = eigEnergy[eigEnergy.size() - 1] - eigEnergy[0];
            groundEnergy = eigEnergy[0];
            // Normalize eigenenergies
            std::for_each(eigEnergy.begin(), eigEnergy.end(),
                          [&](auto& x) { x = (x - groundEnergy) / EnergyRange; });
            debug_print("eigEnergy\n" << eigEnergy);
 
            MicroCanonicalAverage const MCaverage(eigEnergy, shellWidth);
 
            temp_t = getETtime();
            std::cout << "Sample No. " << repetition << ", L = " << L << std::endl;
            for(int opRep = opRepMin; opRep <= opRepMax; ++opRep) {
                for(int m = mMin; m <= mMax; ++m) {
                    std::cout << "opRep=" << opRep << ", m=" << m << std::endl;
                    Eigen::MatrixX<Complex_t<double>> testOp
                        = transSector.basis().adjoint() * GRME[m]() * transSector.basis();
 
                    calculateETHmeasure(colVec, eigenSolver.eigenvectors(), transSector,
                                        mBodyOperatorSpace(m, H_tot), MCaverage);
                    ETHmeasure.col(m) += colVec;
                }
                temp_t = getETtime();
                {
                    debug_print("# Writing results to a file.");
                    std::stringstream buff("");
                    buff << "/OpSample_No" << opRep;
                    std::string filename(outDirName + buff.str());
                    if(!std::filesystem::exists(filename)
                       && !std::filesystem::create_directories(filename)) {
                        std::cerr << "Error: Failed to create directories: " << filename
                                  << std::endl;
                        std::exit(EXIT_FAILURE);
                    };
 
                    buff << "/StrongETHMeasure" << PRECISION << "_N" << L << ".txt";
                    filename += buff.str();
                    OutFs.open(filename);
                    if(!OutFs) {
                        std::cerr << "Error: Can't open a file (" << filename << ")" << std::endl;
                        std::exit(EXIT_FAILURE);
                    }
                    OutFs << std::right << std::showpos << std::scientific << std::setprecision(6);
                    OutFs << "# energyRange= " << EnergyRange << "\n"
                          << "# groundEnergy= " << groundEnergy << "\n"
                          << "# shellWidth= " << shellWidth << "\n"
                          << "# diffWithTheory= " << diff << "\n"
                          << "# 1.(Normalized energy) 2.(dimShell) 2.(Distinguishability measure) "
                          << "\n\n";
                    for(auto j = 0; j < eigEnergy.size(); ++j) {
                        OutFs << eigEnergy[j] << " " << MCaverage.dimShell(j);
                        for(auto m = 1; m <= L; ++m) { OutFs << " " << ETHmeasure(j, m); }
                        OutFs << std::endl;
                    }
                    OutFs.close();
                }
                T_IO += getETtime() - temp_t;
 
                std::cout << std::endl;
            }
            T_post += getETtime() - temp_t;
 
            // std::cout << std::endl;
            // theorySum = Eigen::VectorXd::Ones(MCaverage.dim())
            //             - MCaverage.dimShell().cast<double>().cwiseInverse();
            // double 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;
            // }
        }
 
        if(repetition % 10 == 9 || repetition == repMax) {
            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) << "%)"
                      << ", T_IO=" << std::setprecision(6) << std::setw(8) << T_IO << "("
                      << std::setprecision(1) << 100 * T_IO / (end - start) << "%)" << std::endl;
        }
    }
    return EXIT_SUCCESS;
}