72 {
73 if(argc != 9) {
74 std::cerr << "Usage: 1.(This) 2.(dimLoc) 3.(Lmax) 4.(Lmin) 5.(momentum) 6.(shellWidth) "
75 "7.(repMin) 8.(repMax) 9.(OutDir)\n"
76 << "argc = " << argc << std::endl;
77 std::exit(EXIT_FAILURE);
78 }
79#ifdef MAGMA
81
82
83
84 std::size_t heapSize, stackSize;
85 cuCHECK(cudaDeviceGetLimit(&stackSize, cudaLimitStackSize));
86 cuCHECK(cudaDeviceGetLimit(&heapSize, cudaLimitMallocHeapSize));
87 std::cout << "Max stack size = " << stackSize << " bytes \n"
88 << "Max heap size = " << heapSize << " bytes \n"
89 << std::endl;
90#endif
91
92 int const dimLoc = std::atoi(argv[1]);
93 int const Lmax = std::atoi(argv[2]);
94 int const Lmin = std::atoi(argv[3]);
95 int const k = std::atoi(argv[4]);
96 double const shellWidth = std::atof(argv[5]);
97 int const repMin = std::atoi(argv[6]);
98 int const repMax = std::atoi(argv[7]);
102 std::ofstream OutFs;
103
104 double EnergyRange, groundEnergy, diff;
107
108
109
110
111
112 int dimTot;
113 {
115 dimTot = transSector.dim();
116 }
117 std::cout << "\t dimTot = " << dimTot << "\n\n\n" << std::endl;
118
119 Eigen::MatrixXd ETHmeasure;
120 Eigen::VectorXd eigEnergy, theorySum, colVec;
121
123#ifdef GPU
125 Eigen::Matrix<Complex_t<Real_t>, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> >
126 eigenSolver(dimTot);
127#else
129#endif
130
132
133 double start, t_int, end, temp_t;
134 double T_diag = 0, T_post = 0, T_pre = 0, T_IO = 0;
135
136 ensemble.discardSamples(
repMin);
138
140 end = start;
141 for(
auto repetition =
repMin; repetition <=
repMax; ++repetition) {
144 {
145 std::stringstream
buff;
146 buff <<
"/RawData/Sample_No" << repetition;
147 outDirName +=
buff.str();
148 outDirName = std::regex_replace(outDirName, std::regex("//"), "/");
149 if(!std::filesystem::exists(outDirName)
150 && !std::filesystem::create_directories(outDirName)) {
151 std::cerr << "Error: Failed to create directories " << outDirName << std::endl;
152 std::exit(EXIT_FAILURE);
153 };
154 }
156
157 ensemble.getNewSample();
159
160 for(int L = Lmax; L >= Lmin; --L) {
162
165
166 eigenSolver.inMatrix().resize(transSector.dim(), transSector.dim());
167 eigEnergy.resize(transSector.dim());
168 ETHmeasure = Eigen::MatrixXd::Zero(transSector.dim(), L + 1).eval();
169 diff = ETHmeasure.norm();
170 if(diff > 1.0e-4) {
171 std::cerr << "Error : failed to initialize ETHmeasure: diff = " << diff
172 << " is too large." << std::endl;
173 std::exit(EXIT_FAILURE);
174 }
175 debug_print(
"eigenSolver.inMatrix():\n" << eigenSolver.inMatrix());
176 ensemble.constructHamiltonian(eigenSolver.inMatrix(), transSector);
177
179
180
182 { eigenSolver.compute(); }
184
185
186
187
188
189
190 eigEnergy = eigenSolver.eigenvalues().template cast<double>();
191 EnergyRange = eigEnergy[eigEnergy.size() - 1] - eigEnergy[0];
192 groundEnergy = eigEnergy[0];
193
194 std::for_each(eigEnergy.begin(), eigEnergy.end(),
195 [&](auto& x) { x = (x - groundEnergy) / EnergyRange; });
197
199
200 std::cout << "Sample No. " << repetition << ", L = " << L << std::endl;
201 for(int m = 1; m <= L; ++m) {
202 std::cout << "m=" << m << std::endl;
204 calculateETHmeasure(colVec, eigenSolver.eigenvectors(), transSector,
207 ETHmeasure.col(m) += colVec;
208 }
209
210 std::cout << std::endl;
211 theorySum = Eigen::VectorXd::Ones(MCaverage.dim())
212 - MCaverage.dimShell().cast<double>().cwiseInverse();
213 double maxDiff = (ETHmeasure.rowwise().sum() - theorySum).maxCoeff();
214 if(maxDiff > 1.0e-4) {
215 std::cout << "ETHmeasure.rowwise().sum():\n"
216 << ETHmeasure.rowwise().sum() << "\n"
217 << std::endl;
218 std::cout << "theorySum:\n" << theorySum << "\n" << std::endl;
219 std::cout << "ETHmeasure:\n" << ETHmeasure << "\n" << std::endl;
220 std::cerr << "Error : maxDiff = " << maxDiff << " is too large." << std::endl;
221
222 }
223
224
226 {
228 std::stringstream
buff(
"");
229 buff <<
"/RestrictionToMBodyOperators" << PRECISION <<
"_N" << L <<
".txt";
230 std::string filename(outDirName);
231 filename +=
buff.str();
232 OutFs.open(filename);
233 if(!OutFs) {
234 std::cerr << "Error: Can't open a file (" << filename << ")" << std::endl;
235 std::exit(EXIT_FAILURE);
236 }
237 OutFs << std::right << std::showpos << std::scientific << std::setprecision(6);
238 OutFs << "# energyRange= " << EnergyRange << "\n"
239 << "# groundEnergy= " << groundEnergy << "\n"
240 << "# shellWidth= " << shellWidth << "\n"
241 << "# diffWithTheory= " << diff << "\n"
242 << "# 1.(Normalized energy) 2.(dimShell) 2.(ETHmeasure) "
243 << "\n\n";
244 for(auto j = 0; j < eigEnergy.size(); ++j) {
245 OutFs << eigEnergy[j] << " " << MCaverage.dimShell(j);
246 for(auto m = 1; m <= L; ++m) { OutFs << " " << ETHmeasure(j, m); }
247 OutFs << std::endl;
248 }
249 OutFs.close();
250 }
252 }
253
254 if(repetition % 10 == 9 || repetition ==
repMax) {
255 t_int = end;
257 std::cerr << "(total=" << std::setw(6) << repetition + 1
258 << "): timeINT=" << std::setprecision(6) << std::setw(8) << (end - t_int)
259 << ", timeTOT=" << std::setprecision(6) << std::setw(8) << (end - start)
260 << ", T_construct=" << std::setprecision(6) << std::setw(10) << T_pre << "("
261 << std::setprecision(1) << 100 * T_pre / (end - start) << "%)"
262 << ", T_diag=" << std::setprecision(6) << std::setw(8) << T_diag << "("
263 << std::setprecision(1) << 100 * T_diag / (end - start) << "%)"
264 << ", T_process=" << std::setprecision(6) << std::setw(8) << T_post << "("
265 << std::setprecision(1) << 100 * T_post / (end - start) << "%)"
266 << ", T_IO=" << std::setprecision(6) << std::setw(8) << T_IO << "("
267 << std::setprecision(1) << 100 * T_IO / (end - start) << "%)" << std::endl;
268 }
269 }
270 return EXIT_SUCCESS;
271}
double getETtime()
Definition EnergySpectrum.c:14
Definition ETHmeasure.hpp:133
Definition HilbertSpace.hpp:32
static MagmaController & getInstance()
Definition ObjectOnGPU_EigenMatrix_test.cu:14
Definition HilbertSpace.hpp:423
Definition MatrixUtils.hpp:20
Calculate the microcanonical averages with respect to a given sorted vector 'eigVal'.
Definition MicroCanonicalAverage.hpp:25
Definition ObjectOnGPU.cuh:149
Definition Ensemble_old.hpp:6
Translation invariant sector of a many-body Hilbert space.
Definition TransSector.hpp:19
Definition OperatorSpace.hpp:430
cuCHECK(cudaFuncGetAttributes(&attr, MatrixElementsInSector))
debug_print("# Determining GPU configuration.")
baseDirName
Definition setVariablesForEnsemble.cpp:50
Integer_t const repMin
Definition setVariablesForEnsemble.cpp:31
Integer_t const repMax
Definition setVariablesForEnsemble.cpp:32
std::stringstream buff("")