4#include <gsl/gsl_spline.h>
22 for (
size_t T = 0; T <
f.efts.size(); T++) {
23 out <<
"EFT " << T <<
'\n';
24 for (
size_t m = 0; m <
f.efts[T].height; m++) {
26 for (
size_t i = 0; i <
f.efts[T].w->width; i++)
27 out <<
' ' <<
f.y(T,m,i)*
f.efts[T].units[i];
41(
Lattice_model *model,
const vector<SingleSiteConstraint*>& cs,
43 EFT(model, flow), constraints(cs), matching(m)
48 tiny_dy(1
e-2), huge_dy(100),
49 max_a_steps(1024), max_iter(100),
50 units_set(false), hybrid(false), scl0(1),
58 for (
auto i: teqidx)
delete constraints[i];
59 for (
auto i: rgeidx)
delete constraints[i];
63(
Lattice_model *model,
const vector<SingleSiteConstraint*>& constraints)
65 add_model(model,
nullptr, constraints);
70 const vector<SingleSiteConstraint*>& constraints)
72 efts.emplace_back(model, constraints, matching,
this);
77 const vector<SingleSiteConstraint*>& upward_constraints,
78 const vector<SingleSiteConstraint*>& downward_constraints)
80 add_model(model,
nullptr, upward_constraints, downward_constraints);
86 const vector<SingleSiteConstraint*>& us,
87 const vector<SingleSiteConstraint*>& ds)
89 using ci_t = vector<SingleSiteConstraint*>::const_iterator;
91 vector<SingleSiteConstraint*> rs(ds.size());
92 reverse_copy(ds.begin(), ds.end(), rs.begin());
94 vector<SingleSiteConstraint*> combined;
98 while (pu != us.cend() || pr != rs.cend()) {
99 auto pun = find(pu, us.cend(), pr != rs.cend() ? *pr :
nullptr);
100 auto prn = find(pr, rs.cend(), pu != us.cend() ? *pu :
nullptr);
102 if (pun == pu && prn == pr) {
104 combined.push_back(*pu);
107 else if ((pr == rs.cend() && pu != us.cend()) ||
108 (prn == rs.cend() && pun != us.cend()))
110 combined.push_back(*pu++);
111 else if ((pu == us.cend() && pr != rs.cend()) ||
112 (pun == us.cend() && prn != rs.cend()))
114 combined.push_back(*pr++);
117 msg <<
"RGFlow<Lattice>::Error: failed to combine upward "
118 "and downward constraints in model" << model->name();
119 throw SetupError(msg.str());
124 stringstream uss, rss, css;
125 for (
auto c: us ) uss <<
' ' << c;
126 for (
auto c: rs ) rss <<
' ' << c;
127 for (
auto c: combined) css <<
' ' << c;
128 VERBOSE_MSG(
" upward constraints:" << uss.str() <<
"\n"
129 "reversed downward constraints:" << rss.str() <<
"\n"
130 " combined constraints:" << css.str());
133 add_model(model, matching, combined);
153 init_profile = guesser;
154 init_profile->
init(
this);
160 throw SetupError(
"RGFlow<Lattice>::Error: EFT tower empty");
165 if (hybrid) rk_stage();
166 else increase_density();
172 size_t max_width = 0;
174 for (
size_t T = 0; T < efts.size(); T++) {
175 efts[T].w->init(
this, T);
177 for (
auto c: efts[T].constraints)
178 c->init(
this, T, height++);
179 if (efts[T].matching) {
180 efts[T].matching->init(
this, T);
183 efts[T].height = height;
185 efts[T].offset = offset;
186 size_t width = efts[T].w->width;
187 if (width > max_width) max_width = width;
188 offset += width * height;
192 row_pool.resize(offset);
193 init_free_row_list();
195 for (
size_t T = 0; T < efts.size(); T++) {
196 for (
size_t m = 0; m < efts[T].height-1; m++) {
197 teqidx.push_back(constraints.size());
199 teq->
init(
this, T, m, 2);
200 constraints.push_back(teq);
202 rgeidx.push_back(constraints.size());
204 rge->
init(
this, T, m, 2);
205 constraints.push_back(rge);
207 for (
auto c: efts[T].constraints)
208 constraints.push_back(c);
209 if (efts[T].matching)
210 constraints.push_back(efts[T].matching);
213 for (
auto c: constraints) c->alloc_rows();
217 KL = KU = 2*max_width;
222 throw MemoryError(
"RGFlow<Lattice>::Error: failed to allocate matrix");
227 enum { END, INIT, RUN } state = INIT;
230 VERBOSE_MSG(
"\n\nentering a-loop with a_steps=" << a_steps);
231 while (state != END) {
238 for (
size_t a_i = 0; a_i <= a_steps; a_i++) {
239 a =
Real(a_i)/a_steps;
241 a_i <<
'/' << a_steps);
246 if (a_steps <= max_a_steps) {
248 "with doubled a_steps=" << a_steps);
252 throw DivergenceError
253 (
"RGFlow<Lattice>::Error: Newton iterations diverged");
256 if (a_steps <= max_a_steps) {
259 "with doubled a_steps=" << a_steps);
278 VERBOSE_MSG(
"\n\nbeginning to increase lattice density");
282 vector<size_t> old_heights(efts.size());
283 vector<size_t> old_offsets(efts.size());
284 for (
size_t T = 0; T < efts.size(); T++) {
285 old_heights[T] = efts[T].height;
286 old_offsets[T] = efts[T].offset;
288 vector<vector<size_t>> site_maps = refine_lattice();
289 stringstream heights;
291 for (
auto&
e: efts) heights <<
' ' <<
e.height;
292 VERBOSE_MSG(
"\n\nentering Newton loop with heights:" <<
298 throw DivergenceError(
"RGFlow<Lattice>::Error: RG flow jumped, "
299 "this cannot happen");
304 for (
size_t T = 0; T < efts.size(); T++)
305 for (
size_t m = 0; m < old_heights[T]; m++)
306 for (
size_t i = 0; i < efts[T].w->width; i++) {
308 fabs(y(T,site_maps[T][m],i) -
309 old_y[old_offsets[T]+m*efts[T].w->width+i]);
310 if (diff > maxdy) maxdy = diff;
313 if (maxdy < tiny_dy) {
324 size_t T_top = efts.size() - 1;
325 Real t_top = x(T_top,efts[T_top].height-1,0);
326 Real t_bot = x(0,0,0);
327 Real min_Dt = t_top - t_bot + 1;
329 for (
size_t T = 0; T < efts.size(); T++) {
330 for (
size_t m = 0; m < efts[T].height-1; m++) {
331 Real Dt = x(T,m+1,0) - x(T,m,0);
332 if (Dt > max_Dt) max_Dt = Dt;
333 if (Dt < min_Dt) min_Dt = Dt;
337 Real spacing = max_Dt > 2*min_Dt ? min_Dt : min_Dt/2;
338 vector<vector<size_t>> site_maps(efts.size());
339 for (
size_t T = 0; T < efts.size(); T++) {
341 site_maps[T].push_back(new_m);
342 for (
size_t m = 0; m < efts[T].height-1; m++) {
343 Real Dt = x(T,m+1,0) - x(T,m,0);
344 size_t q = Dt / spacing + 0.5;
345 site_maps[T].push_back(new_m += q);
357 enable_Runge_Kutta();
362 throw DivergenceError(
"RGFlow<Lattice>::Error: RG flow jumped, "
363 "this cannot happen");
376 vector<vector<bool>> original(efts.size());
377 for (
size_t T = efts.size(); T--; ) {
378 original[T].resize(efts[T].height);
379 for (
auto c: efts[T].constraints)
380 original[T][c->mbegin] =
true;
381 if (efts[T].matching)
382 original[T][efts[T].height-1] = original[T+1][0] =
true;
385 vector<vector<size_t>> site_maps(efts.size());
386 for (
size_t T = 0; T < efts.size(); T++) {
387 site_maps[T].resize(efts[T].height);
388 for (
size_t m = 0, new_m = 0; m < efts[T].height; m++) {
389 site_maps[T][m] = new_m;
390 if (original[T][m]) new_m++;
397 for (
size_t i = 0, T = 0; T < efts.size(); T++)
398 for (
size_t m = 0; m < efts[T].height-1; m++, i++) {
399 constraints[rgeidx[i]]->free_rows();
400 constraints[rgeidx[i]]->deactivate();
401 delete constraints[rgeidx[i]];
403 rkrge->
init(
this, T, m, 2);
404 constraints[rgeidx[i]] = rkrge;
406 for (
auto i: rgeidx) constraints[i]->alloc_rows();
417 for (
size_t i = 0, T = 0; T < efts.size(); T++)
418 for (
size_t m = 0; m < efts[T].height-1; m++, i++) {
419 constraints[rgeidx[i]]->free_rows();
420 constraints[rgeidx[i]]->deactivate();
421 delete constraints[rgeidx[i]];
423 rge->
init(
this, T, m, 2);
424 constraints[rgeidx[i]] = rge;
426 for (
auto i: rgeidx) constraints[i]->alloc_rows();
434 vector<size_t> new_heights(efts.size());
435 vector<size_t> new_offsets(efts.size());
437 for (
size_t T = 0; T < efts.size(); T++) {
438 new_heights[T] = site_maps[T][efts[T].height-1] + 1;
439 new_offsets[T] = offset;
440 offset += efts[T].w->width * new_heights[T];
443 for (
size_t T = 0; T < efts.size(); T++) {
445 maps <<
"EFT " << T <<
":";
446 for (
size_t m = 0; m < efts[T].height; m++)
447 maps <<
' ' << m <<
"->" << site_maps[T][m];
453 for (
size_t T = 0; T < efts.size(); T++) {
454 RVec y0(efts[T].height);
455 for (
size_t m = 0; m < efts[T].height; m++) {
456 size_t new_m = site_maps[T][m];
457 y0[m] = new_y[new_offsets[T] + new_m*efts[T].w->width] = y(T,m,0);
459 for (
size_t m = 0; m < efts[T].height-1; m++) {
460 size_t q = site_maps[T][m+1] - site_maps[T][m];
461 for (
size_t n = 1; n < q; n++) {
462 size_t new_m = site_maps[T][m] + n;
463 new_y[new_offsets[T] + new_m*efts[T].w->width] =
464 ((q-n)*y(T,m,0) + n*y(T,m+1,0)) / q;
467 RVec yi(efts[T].height);
468 for (
size_t i = 1; i < efts[T].w->width; i++) {
469 gsl_interp_accel *acc = gsl_interp_accel_alloc();
470 const gsl_interp_type *itype =
471 efts[T].height > 2 ? gsl_interp_cspline : gsl_interp_linear;
472 gsl_spline *spl = gsl_spline_alloc(itype, efts[T].height);
473 for (
size_t m = 0; m < efts[T].height; m++) yi[m] = y(T,m,i);
474 gsl_spline_init(spl, &y0[0], &yi[0], efts[T].height);
475 for (
size_t m = 0; m < new_heights[T]; m++)
476 new_y[new_offsets[T] + m*efts[T].w->width + i] =
478 (spl, new_y[new_offsets[T] + m*efts[T].w->width], acc);
479 gsl_spline_free(spl);
480 gsl_interp_accel_free(acc);
484 for (
auto c: constraints) c->free_rows();
487 row_pool.resize(offset);
488 init_free_row_list();
489 for (
auto c: constraints) c->relocate(site_maps);
490 for (
size_t T = 0; T < efts.size(); T++) {
491 efts[T].height = new_heights[T];
492 efts[T].offset = new_offsets[T];
494 for (
auto c: constraints) c->alloc_rows();
502 throw MemoryError(
"RGFlow<Lattice>::Error: failed to allocate matrix");
507 for (
size_t T = 0; T < efts.size(); T++) {
508 for (
size_t i = 0; i < efts[T].w->width; i++) {
509 Real miny = y(T,0,i);
511 for (
size_t m = 0; m < efts[T].height; m++) {
512 if (y(T,m,i) < miny) miny = y(T,m,i);
513 if (y(T,m,i) > maxy) maxy = y(T,m,i);
515 efts[T].units[i] = max(maxy-miny,max(fabs(maxy),fabs(miny)));
516 if (efts[T].units[i] == 0) {
517 VERBOSE_MSG(
"no hint about unit of x[" << i <<
"] in EFT "
518 << T <<
", defaulting to 1");
519 efts[T].units[i] = 1;
521 for (
size_t m = 0; m < efts[T].height; m++)
522 y(T,m,i) /= efts[T].units[i];
531 if (multithreading) {
532 threads_begin->wait();
533 (**elementary_constraints.begin())();
537 for (
auto c: elementary_constraints) (*c)();
544 while (threads_begin->wait(), keep_threads) {
552 if (!multithreading)
return;
554 threads_begin =
new boost::barrier(elementary_constraints.size());
555 threads_end =
new boost::barrier(elementary_constraints.size());
557 threads =
new boost::thread_group;
558 for (
auto c: elementary_constraints)
559 if (c != *elementary_constraints.begin())
563 VERBOSE_MSG(
"launched " << threads->size() <<
" subthreads");
568 if (!multithreading)
return;
570 keep_threads =
false;
571 threads_begin->wait();
573 VERBOSE_MSG(
"joined " << threads->size() <<
" subthreads");
575 delete threads_begin;
581 assert(y0.size() == y1.size());
585 while (p < y0.end()) {
586 Real diff = fabs(*p++ - *q++);
587 if (diff > max) max = diff;
593(
size_t T,
size_t m,
size_t span)
595 EqRow *r = free_row_list_head;
596 assert(r !=
nullptr);
598 free_row_list_head = r->next;
599 r->rowSpec = { T, m, span, 0 };
606 r->next = free_row_list_head;
607 free_row_list_head = r;
612 if (free_row_list_head !=
nullptr)
613 throw SetupError(
"RGFlow<Lattice>::Error: constraints not enough");
614 vector<EqRow *> layout(row_pool.size());
615 for (
size_t r = 0; r < layout.size(); r++)
616 layout[r] = &row_pool[r];
617 sort(layout.begin(), layout.end(),
618 [](EqRow *a, EqRow *b) ->
bool {
619 if (a->rowSpec.T < b->rowSpec.T) return true;
620 if (a->rowSpec.T > b->rowSpec.T) return false;
621 if (a->rowSpec.m < b->rowSpec.m) return true;
622 if (a->rowSpec.m > b->rowSpec.m) return false;
623 if (a->rowSpec.n < b->rowSpec.n) return true;
626 for (
size_t r = 0; r < layout.size(); r++)
627 layout[r]->rowSpec.realRow = r;
632 for (
size_t r = 0; r < row_pool.size() - 1; r++)
633 row_pool[r].next = &row_pool[r+1];
634 row_pool.back().next =
nullptr;
635 free_row_list_head = &row_pool[0];
639(
const int& N,
const int& KL,
const int& KU,
const int& NRHS,
640 double *AB,
const int& LDAB,
int *IPIV,
double *
B,
const int& LDB,
int *
INFO);
644 if (!units_set) { set_units(); units_set =
true; }
649 enum { END, INIT } state = INIT;
655 while (iter < max_iter && state != END) {
659 cout << setprecision(15) << scientific;
660 for (
size_t i = 0; i < y_.size(); i++) {
661 for (
size_t j = 0; j < y_.size(); j++) {
663 if (abs(
signed(i)-
signed(j)) < KL) cout << (*A_)(i,j);
669 dgbsv_(N,KL,KU,NRHS,A_->pointer(),LDA,&IPIV[0],&z[0],LDB,&
INFO);
673 msg <<
"RGFlow<Lattice>::Error: failed to solve equations, "
674 "DGBSV returned INFO = " <<
INFO;
675 throw NonInvertibleMatrixError(msg.str());
677 Real maxdy = maxdiff(y_, z);
678 VERBOSE_MSG(
"iter=" << iter <<
" maxdy=" << maxdy);
681 else if (maxdy > huge_dy)
687 VERBOSE_MSG(
"converged after " << iter <<
" Newton iterations");
virtual void init(RGFlow< Lattice > *flow)
virtual void init(RGFlow< Lattice > *flow, size_t theory, size_t site, size_t span_)
void init(RGFlow< Lattice > *flow, size_t theory, size_t site, size_t span)
void init(RGFlow< Lattice > *flow, size_t theory, size_t site, size_t span)
virtual int run_to(double, double eps=-1.0)
No convergence while solving the RGEs.
void add_model(Lattice_model *model, const std::vector< SingleSiteConstraint * > &constraints)
Spectrum generator was not setup correctly.
void() B(bb)) coeff_t &arr
std::complex< double > f(double tau) noexcept
Matching< Lattice > InterTheoryConstraint
std::ostream & operator<<(std::ostream &ostr, const Dynamic_array_view< ElementType > &av)
void sort(Eigen::Array< double, N, 1 > &v)
sorts an Eigen array
void dgbsv_(const int &N, const int &KL, const int &KU, const int &NRHS, double *AB, const int &LDAB, int *IPIV, double *B, const int &LDB, int *INFO)
Integration of ODEs by Runge-Kutta.