jscaux
/
ABACUS
2
0
Fork 0
Code Issues 1 Pull-Requests Releases Wiki Aktivität
ABACUS/dev/ODSLF/ODSLF.cc

1429 Zeilen
47 KiB
C++
Originalformat Permalink Blame Verlauf

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: src/1DSLF/1DSLF.cc
Purpose: defines functions for 1d spinless fermion classes.
***********************************************************/
#include "ABACUS.h"
using namespace std;
namespace ABACUS {
//***************************************************************************************************
// Function definitions: class ODSLF_Base
ODSLF_Base::ODSLF_Base () : Mdown(0), Nrap(Vect<int>()), Nraptot(0), Ix2_infty(Vect<DP>()), Ix2_max(Vect<int>()) {}
ODSLF_Base::ODSLF_Base (const ODSLF_Base& RefBase) // copy constructor
: Mdown(RefBase.Mdown), Nrap(Vect<int>(RefBase.Nrap.size())), Nraptot(RefBase.Nraptot),
Ix2_infty(Vect<DP>(RefBase.Nrap.size())), Ix2_max(Vect<int>(RefBase.Nrap.size())),
id(RefBase.id)
{
for (int i = 0; i < Nrap.size(); ++i) {
Nrap[i] = RefBase.Nrap[i];
Ix2_infty[i] = RefBase.Ix2_infty[i];
Ix2_max[i] = RefBase.Ix2_max[i];
}
}
ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, int M)
: Mdown(M), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_max(Vect<int>(RefChain.Nstrings))
{
for (int i = 0; i < RefChain.Nstrings; ++i) Nrap[i] = 0;
Nrap[0] = M;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
// The id of this is zero by definition
id = 0LL;
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
}
ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities)
: Mdown(0), Nrap(Nrapidities), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_max(Vect<int>(RefChain.Nstrings)), id (0LL)
{
// Check consistency of Nrapidities vector with RefChain
if (RefChain.Nstrings != Nrapidities.size())
ABACUSerror("Incompatible Nrapidities vector used in ODSLF_Base constructor.");
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
Mdown = Mcheck;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
// Compute id
id += Nrapidities[0];
long long int factor = 100000LL;
for (int i = 1; i < RefChain.Nstrings; ++i) {
id += factor * Nrapidities[i];
factor *= 100LL;
}
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
}
ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, long long int id_ref)
: Mdown(0), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_max(Vect<int>(RefChain.Nstrings)), id (id_ref)
{
// Build Nrapidities vector from id_ref
long long int factor = pow_ulli (10LL, 2* RefChain.Nstrings + 1);
long long int id_eff = id_ref;
for (int i = 0; i < RefChain.Nstrings - 1; ++i) {
Nrap[RefChain.Nstrings - 1 - i] = id_eff/factor;
id_eff -= factor * Nrap[RefChain.Nstrings - 1 - i];
factor /= 100LL;
}
Nrap[0] = id_eff;
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
Mdown = Mcheck;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
}
ODSLF_Base& ODSLF_Base::operator= (const ODSLF_Base& RefBase)
{
if (this != & RefBase) {
Mdown = RefBase.Mdown;
Nrap = RefBase.Nrap;
Nraptot = RefBase.Nraptot;
Ix2_infty = RefBase.Ix2_infty;
Ix2_max = RefBase.Ix2_max;
id = RefBase.id;
}
return(*this);
}
bool ODSLF_Base::operator== (const ODSLF_Base& RefBase)
{
bool answer = (Nrap == RefBase.Nrap);
return (answer);
}
bool ODSLF_Base::operator!= (const ODSLF_Base& RefBase)
{
bool answer = (Nrap != RefBase.Nrap);
return (answer);
}
void ODSLF_Base::Compute_Ix2_limits (const Heis_Chain& RefChain)
{
if ((RefChain.Delta > 0.0) && (RefChain.Delta < 1.0)) {
// Compute the Ix2_infty numbers
DP sum1 = 0.0;
DP sum2 = 0.0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0.0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum2 = 0.0;
sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 :
2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));
for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));
sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2;
}
Ix2_infty[j] = (1.0/PI) * fabs(RefChain.Nsites * 2.0
* atan(tan(0.25 * PI * (1.0 + RefChain.par[j])
- 0.5 * RefChain.Str_L[j] * RefChain.anis)) - sum1);
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Reject formally infinite rapidities (i.e. if Delta is root of unity)
// For ODSLF: parity of quantum numbers is complicated:
// if (N-1 + M_j + (M + 1) n_j + (1 - nu_j) N/2) is even, q# are integers,
// so Ix2_max must be even (conversely, if odd, then odd).
if ((Ix2_max[j] + RefChain.Nsites - 1 + Nrap[j] + (Nraptot + 1) * RefChain.Str_L[j]
+ (RefChain.par[j] == 1 ? 0 : RefChain.Nsites)) % 2) Ix2_max[j] -= 1;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
}
} // if XXZ gapless
else if (RefChain.Delta == 1.0) {
// Compute the Ix2_infty numbers
int sum1 = 0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
}
Ix2_infty[j] = (RefChain.Nsites + 1.0 - sum1); // to get counting right...
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Give the correct parity to Ix2_max
// If Nrap is even, Ix2_max must be odd. If odd, then even.
//if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1; DEACTIVATED FOR ODSLF !
// If Ix2_max equals Ix2_infty, we reduce it by 2:
if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
}
} // if XXX AFM
else if (RefChain.Delta > 1.0) {
// Compute the Ix2_infty numbers
int sum1 = 0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
}
Ix2_infty[j] = (RefChain.Nsites - 1 + 2 * RefChain.Str_L[j] - sum1);
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Give the correct parity to Ix2_max
// If Nrap is even, Ix2_max must be odd. If odd, then even.
//if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1; DEACTIVATED FOR ODSLF !
// If Ix2_max equals Ix2_infty, we reduce it by 2:
//if (Ix2_max[j] == Ix2_infty[j]) Ix2_max[j] -= 2;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
}
} // if XXZ_gpd
}
void ODSLF_Base::Scan_for_Possible_Types (Vect<long long int>& possible_type_id, int& nfound,
int base_level, Vect<int>& Nexcitations)
{
if (base_level == 0) {
// We set all sectors 0, 1, 2, 3
int Nexc_max_0 = (Nrap[0] + 1)/2; // holes on R
int Nexc_max_1 = Nrap[0]/2;
int Nexc_max_2 = (Ix2_max[0] - Nrap[0] + 1)/2;
int Nexc_max_3 = Nexc_max_2;
// Count the number of excitations in higher levels:
int Nexc_above = 0;
for (int i = 4; i < Nexcitations.size(); ++i) Nexc_above += Nexcitations[i];
for (int nphpairs = 0; nphpairs <= NEXC_MAX_HEIS/2 - Nexc_above; ++nphpairs)
for (int nexc0 = 0; nexc0 <= ABACUS::min(Nexc_max_0, nphpairs); ++nexc0) {
int nexc1 = nphpairs - nexc0;
for (int nexc2 = 0; nexc2 <= ABACUS::min(Nexc_max_2, nphpairs); ++nexc2) {
int nexc3 = nphpairs - nexc2;
if (nexc1 >= 0 && nexc1 <= Nexc_max_1 && nexc3 >= 0 && nexc3 <= Nexc_max_3) {
// acceptable type found
Nexcitations[0] = nexc0; Nexcitations[1] = nexc1; Nexcitations[2] = nexc2; Nexcitations[3] = nexc3;
possible_type_id[nfound] = Ix2_Offsets_type_id (Nexcitations);
nfound++;
}
}
}
}
else { // base_level > 0
// We scan over the possible distributions of the excitations on the L and R sides:
// Count the number of slots available on R:
int Nexc_R_max = (Ix2_max[base_level] + 2)/2;
int Nexc_L_max = (Ix2_max[base_level] + 1)/2;
for (int Nexcleft = 0; Nexcleft <= Nrap[base_level]; ++Nexcleft)
if (Nexcleft <= Nexc_L_max && Nrap[base_level] - Nexcleft <= Nexc_R_max) {
Nexcitations[2*base_level + 2] = Nrap[base_level] - Nexcleft;
Nexcitations[2*base_level + 3] = Nexcleft;
Scan_for_Possible_Types (possible_type_id, nfound, base_level - 1, Nexcitations);
}
}
}
Vect<long long int> ODSLF_Base::Possible_Types ()
{
// Given a base, this returns a vector of possible types
Vect<long long int> possible_type_id (0LL, 1000);
int nfound = 0;
Vect<int> Nexcitations (0, 2*Nrap.size() + 2);
Scan_for_Possible_Types (possible_type_id, nfound, Nrap.size() - 1, Nexcitations);
// Copy results into new vector:
Vect<long long int> possible_type_id_found (nfound);
for (int i = 0; i < nfound; ++i)
possible_type_id_found[i] = possible_type_id[i];
return(possible_type_id_found);
}
//***************************************************************************************************
// Function definitions: class Ix2_Config
ODSLF_Ix2_Config::ODSLF_Ix2_Config () {}
ODSLF_Ix2_Config::ODSLF_Ix2_Config (const Heis_Chain& RefChain, int M)
: Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), Ix2(new int*[1]) // single type of string here
{
Ix2[0] = new int[M];
for (int j = 0; j < M; ++j) {
Ix2[0][j] = -(M - 1) + 2*j;
// Use simplification: Nrap[0] = M = Nraptot, Str_L[0] = 1, par = 1, so
if ((Ix2[0][j] + RefChain.Nsites) % 2) Ix2[0][j] -= 1;
}
}
ODSLF_Ix2_Config::ODSLF_Ix2_Config (const Heis_Chain& RefChain, const ODSLF_Base& base)
: Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), Ix2(new int*[RefChain.Nstrings])
{
Ix2[0] = new int[base.Nraptot];
for (int i = 1; i < RefChain.Nstrings; ++i) Ix2[i] = Ix2[i-1] + base[i-1];
// And now for the other string types...
for (int i = 0; i < RefChain.Nstrings; ++i) {
for (int j = 0; j < base[i]; ++j) Ix2[i][j] = -(base[i] - 1) + 2*j + 1 - (base[i] % 2);
}
// ODSLF: put back correct parity of quantum nr:
for (int i = 0; i < RefChain.Nstrings; ++i)
for (int j = 0; j < base[i]; ++j)
if ((Ix2[i][j] + RefChain.Nsites - 1 + base.Nrap[i] + (base.Nraptot + 1) * RefChain.Str_L[i] + (RefChain.par[i] == 1 ? 0 : RefChain.Nsites)) % 2)
Ix2[i][j] -= 1;
}
ODSLF_Ix2_Config& ODSLF_Ix2_Config::operator= (const ODSLF_Ix2_Config& RefConfig)
{
if (this != &RefConfig) {
Nstrings = RefConfig.Nstrings;
Nrap = RefConfig.Nrap;
Nraptot = RefConfig.Nraptot;
if (Ix2 != 0) {
delete[] (Ix2[0]);
delete[] (Ix2);
}
Ix2 = new int*[Nstrings];
Ix2[0] = new int[Nraptot];
for (int i = 1; i < Nstrings; ++i) Ix2[i] = Ix2[i-1] + Nrap[i-1];
for (int i = 0; i < Nstrings; ++i)
for (int j = 0; j < Nrap[i]; ++j) Ix2[i][j] = RefConfig.Ix2[i][j];
}
return(*this);
}
ODSLF_Ix2_Config::~ODSLF_Ix2_Config()
{
if (Ix2 != 0) {
delete[] (Ix2[0]);
delete[] Ix2;
}
}
std::ostream& operator<< (std::ostream& s, const ODSLF_Ix2_Config& RefConfig)
{
s << &RefConfig << "\t" << RefConfig.Nstrings << "\t" << RefConfig.Nraptot << endl;
for (int i = 0; i < RefConfig.Nstrings; ++i) s << RefConfig.Nrap[i] << "\t";
s << endl;
for (int i = 0; i < RefConfig.Nstrings; ++i) {
s << "i " << i << "\t";
for (int j = 0; j < RefConfig.Nrap[i]; ++j) s << RefConfig.Ix2[i][j] << "\t";
s << endl;
}
return(s);
}
//***************************************************************************************************
// Function definitions: class Lambda
ODSLF_Lambda::ODSLF_Lambda () {}
ODSLF_Lambda::ODSLF_Lambda (const Heis_Chain& RefChain, int M)
: Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), lambda(new DP*[1]) // single type of string here
{
lambda[0] = new DP[M];
for (int j = 0; j < M; ++j) lambda[0][j] = 0.0;
}
ODSLF_Lambda::ODSLF_Lambda (const Heis_Chain& RefChain, const ODSLF_Base& base)
: Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), lambda(new DP*[RefChain.Nstrings])
{
lambda[0] = new DP[base.Nraptot];
for (int i = 1; i < RefChain.Nstrings; ++i) lambda[i] = lambda[i-1] + base[i-1];
for (int i = 0; i < RefChain.Nstrings; ++i) {
for (int j = 0; j < base[i]; ++j) lambda[i][j] = 0.0;
}
}
ODSLF_Lambda& ODSLF_Lambda::operator= (const ODSLF_Lambda& RefLambda)
{
if (this != &RefLambda) {
Nstrings = RefLambda.Nstrings;
Nrap = RefLambda.Nrap;
Nraptot = RefLambda.Nraptot;
if (lambda != 0) {
delete[] (lambda[0]);
delete[] (lambda);
}
lambda = new DP*[Nstrings];
lambda[0] = new DP[Nraptot];
for (int i = 1; i < Nstrings; ++i) lambda[i] = lambda[i-1] + Nrap[i-1];
for (int i = 0; i < Nstrings; ++i)
for (int j = 0; j < Nrap[i]; ++j) lambda[i][j] = RefLambda.lambda[i][j];
}
return(*this);
}
ODSLF_Lambda::~ODSLF_Lambda()
{
if (lambda != 0) {
delete[] (lambda[0]);
delete[] lambda;
}
}
//***************************************************************************************************
// Function definitions: class ODSLF_Ix2_Offsets
ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets () : base(), Tableau(Vect<Young_Tableau>()), type_id(0LL), id(0LL), maxid(0LL) {};
ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset) // copy constructor
: base(RefOffset.base), Tableau(Vect<Young_Tableau> (2 * base.Nrap.size() + 2)), type_id(RefOffset.type_id), id(RefOffset.id), maxid(RefOffset.maxid)
{
for (int i = 0; i < 2 * base.Nrap.size() + 2; ++i) Tableau[i] = RefOffset.Tableau[i];
}
ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, long long int req_type_id)
// sets all tableaux to empty ones, with nparticles(req_type_id) at each level
{
// Build nparticles vector from req_type_id
Vect<int> nparticles(0, 2* RefBase.Nrap.size() + 2);
long long int factor = pow_ulli (10LL, nparticles.size() - 1);
long long int id_eff = req_type_id;
for (int i = 0; i < nparticles.size(); ++i) {
nparticles[nparticles.size() - 1 - i] = id_eff/factor;
id_eff -= factor * nparticles[nparticles.size() - 1 - i];
factor /= 10LL;
}
// Check if we've got the right vector...
long long int idcheck = Ix2_Offsets_type_id (nparticles);
if (idcheck != req_type_id) ABACUSerror("idcheck != req_type_id in Ix2_Offsets constructor.");
(*this) = ODSLF_Ix2_Offsets(RefBase, nparticles);
}
ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, Vect<int> nparticles)
// sets all tableaux to empty ones, with nparticles at each level
: base(RefBase), Tableau(Vect<Young_Tableau> (2 * base.Nrap.size() + 2)),
type_id(Ix2_Offsets_type_id (nparticles)), id(0LL), maxid(0LL)
{
// Checks on nparticles vector:
if (nparticles.size() != 2 * base.Nrap.size() + 2)
ABACUSerror("Wrong nparticles.size in Ix2_Offsets constructor.");
if (nparticles[3] + nparticles[2] != nparticles[0] + nparticles[1]) {
cout << nparticles[0] << "\t" << nparticles[1] << "\t" << nparticles[2] << "\t" << nparticles[3] << endl;
ABACUSerror("Wrong Npar[0-3] in Ix2_Offsets constructor.");
}
for (int base_level = 1; base_level < base.Nrap.size(); ++ base_level)
if (base.Nrap[base_level] != nparticles[2*base_level + 2] + nparticles[2*base_level + 3]) {
cout << base_level << "\t" << base.Nrap[base_level] << "\t" << nparticles[2*base_level + 2]
<< "\t" << nparticles[2*base_level + 3] << endl;
ABACUSerror("Wrong Nrap[] in Ix2_Offsets constructor.");
}
// nparticles[0,1]: number of holes on R and L side in GS interval
if (nparticles[0] > (base.Nrap[0] + 1)/2)
ABACUSerror("nparticles[0] too large in Ix2_Offsets constructor.");
if (nparticles[1] > base.Nrap[0]/2)
ABACUSerror("nparticles[1] too large in Ix2_Offsets constructor.");
// nparticles[2,3]: number of particles of type 0 on R and L side out of GS interval
if (nparticles[2] > (base.Ix2_max[0] - base.Nrap[0] + 1)/2)
ABACUSerror("nparticles[2] too large in Ix2_Offsets constructor.");
if (nparticles[3] > (base.Ix2_max[0] - base.Nrap[0] + 1)/2)
ABACUSerror("nparticles[3] too large in Ix2_Offsets constructor.");
for (int base_level = 1; base_level < base.Nrap.size(); ++ base_level)
if ((nparticles[2*base_level + 2] > 0 && nparticles[2*base_level + 2]
> (base.Ix2_max[base_level] + 2)/2) // ODSLF modif
|| (nparticles[2*base_level + 3] > 0
&& nparticles[2*base_level + 3] > (base.Ix2_max[base_level] + 1)/2)) // ODSLF modif
{
cout << base_level << "\t" << base.Ix2_max[base_level] << "\t" << base.Ix2_infty[base_level]
<< "\t" << nparticles[2*base_level + 2] << "\t"
<< (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2
<< "\t" << nparticles[2*base_level + 3] << "\t"
<< (base.Ix2_max[base_level] - (base.Nrap[base_level] % 2) - 1)/2
<< "\t" << (nparticles[2*base_level + 2] > 0) << "\t"
<< (nparticles[2*base_level + 2] > (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2)
<< "\t" << (nparticles[2*base_level + 3] > 0) << "\t"
<< (nparticles[2*base_level + 3] > base.Ix2_max[base_level] + 1
- (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2)
<< endl;
ABACUSerror("nparticles too large in Ix2_Offsets constructor.");
}
// Check sum of rapidities
// Holes in GS interval
Tableau[0] = Young_Tableau(nparticles[0], (base.Nrap[0] + 1)/2 - nparticles[0]);
Tableau[1] = Young_Tableau(nparticles[1], base.Nrap[0]/2 - nparticles[1], Tableau[0]);
// Particles of type 0 out of GS interval
Tableau[2] = Young_Tableau(nparticles[2], (base.Ix2_max[0] - base.Nrap[0] + 1)/2 - nparticles[2], Tableau[0]);
Tableau[3] = Young_Tableau(nparticles[3], (base.Ix2_max[0] - base.Nrap[0] + 1)/2 - nparticles[3], Tableau[2]);
// Tableaux of index i = 2,...: data about string type i/2-1.
for (int base_level = 1; base_level < base.Nrap.size(); ++base_level) {
Tableau[2*base_level + 2] =
Young_Tableau(nparticles[2*base_level + 2],
(base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2))/2 + 1
- nparticles[2*base_level + 2], Tableau[2]);
Tableau[2*base_level + 3] =
Young_Tableau(nparticles[2*base_level + 3],
(base.Ix2_max[base_level] - (base.Nrap[base_level] % 2) - 1)/2 + 1
- nparticles[2*base_level + 3], Tableau[3]);
}
maxid = 1LL;
for (int i = 0; i < nparticles.size(); ++i) {
maxid *= Tableau[i].maxid + 1LL;
}
maxid -= 1LL;
}
ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets::operator= (const ODSLF_Ix2_Offsets& RefOffset)
{
if (this != &RefOffset) {
base = RefOffset.base;
Tableau = RefOffset.Tableau;
type_id = RefOffset.type_id;
id = RefOffset.id;
maxid = RefOffset.maxid;
}
return(*this);
}
bool ODSLF_Ix2_Offsets::operator<= (const ODSLF_Ix2_Offsets& RefOffsets)
{
// Check whether all nonzero tableau row lengths in RefOffsets
// are <= than those in *this
bool answer = true;
for (int level = 0; level < 4; ++level) { // check fundamental level only
// First check whether all rows which exist in both tableaux satisfy rule:
for (int tableau_level = 0; tableau_level
< ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows); ++tableau_level)
if (Tableau[level].Row_L[tableau_level] > RefOffsets.Tableau[level].Row_L[tableau_level])
answer = false;
// Now check whether there exist extra rows violating rule:
for (int tableau_level = ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows);
tableau_level < Tableau[level].Nrows; ++tableau_level)
if (Tableau[level].Row_L[tableau_level] > 0) answer = false;
}
return(answer);
}
bool ODSLF_Ix2_Offsets::operator>= (const ODSLF_Ix2_Offsets& RefOffsets)
{
// Check whether all nonzero tableau row lengths in RefOffsets
// are >= than those in *this
bool answer = true;
for (int level = 0; level < 4; ++level) { // check fundamental level only
// First check whether all rows which exist in both tableaux satisfy rule:
for (int tableau_level = 0; tableau_level
< ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows); ++tableau_level)
if (Tableau[level].Row_L[tableau_level] < RefOffsets.Tableau[level].Row_L[tableau_level])
answer = false;
// Now check whether there exist extra rows violating rule:
for (int tableau_level = ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows);
tableau_level < RefOffsets.Tableau[level].Nrows; ++tableau_level)
if (RefOffsets.Tableau[level].Row_L[tableau_level] > 0) answer = false;
}
return(answer);
}
void ODSLF_Ix2_Offsets::Set_to_id (long long int idnr)
{
// The idnr of the Offset is given by
// sub_id[0] + (total number of tableaux of type 0) * (sub_id[1] + (total number of tableaux of type 1) * (sub_id[2] + ...
// + total number of tableaux of type (2*base.Nrap.size()) * sub_id[2*base.Nrap.size() + 1]
if (idnr > maxid) {
cout << idnr << "\t" << maxid << endl;
ABACUSerror("idnr too large in offsets.Set_to_id.");
}
id = idnr;
Vect<long long int> sub_id(0LL, 2*base.Nrap.size() + 2);
long long int idnr_eff = idnr;
long long int temp_prod = 1LL;
Vect<long long int> result_choose(2*base.Nrap.size() + 2);
for (int i = 0; i <= 2*base.Nrap.size(); ++i) {
result_choose[i] = Tableau[i].maxid + 1LL;
temp_prod *= result_choose[i];
}
for (int i = 2*base.Nrap.size() + 1; i > 0; --i) {
sub_id[i] = idnr_eff/temp_prod;
idnr_eff -= sub_id[i] * temp_prod;
temp_prod /= result_choose[i-1];
}
sub_id[0] = idnr_eff; // what's left goes to the bottom...
for (int i = 0; i <= 2*base.Nrap.size() + 1; ++i) {
if ((Tableau[i].Nrows * Tableau[i].Ncols == 0) && (sub_id[i] != 0))
ABACUSerror("index too large in offset.Set_to_id.");
if (Tableau[i].id != sub_id[i]) Tableau[i].Set_to_id(sub_id[i]);
}
Compute_type_id ();
return;
}
void ODSLF_Ix2_Offsets::Compute_type_id ()
{
type_id = 0LL;
for (int i = 0; i < 2*base.Nrap.size() + 2; ++i) {
Tableau[i].Compute_id();
type_id += Tableau[i].Nrows * pow_ulli(10LL, i);
}
}
void ODSLF_Ix2_Offsets::Compute_id ()
{
long long int prod_maxid = 1LL;
id = 0LL;
for (int i = 0; i < 2*base.Nrap.size() + 2; ++i) {
Tableau[i].Compute_id();
id += Tableau[i].id * prod_maxid;
prod_maxid *= Tableau[i].maxid + 1LL;
}
}
bool ODSLF_Ix2_Offsets::Add_Boxes_From_Lowest (int Nboxes, bool odd_sectors)
{
// Tries to add Nboxes to Young tableau in given sector parity.
// If Nboxes is too large to be accommodated, `false' is returned.
ODSLF_Ix2_Offsets offsets_init(*this); // keep a copy in case it fails
if (Nboxes < 0) ABACUSerror("Can't use Nboxes < 0 in Add_Boxes_From_Lowest");
else if (Nboxes == 0) return(true); // nothing to do
int Nboxes_put = 0;
int working_sector = odd_sectors;
// We start from the bottom up.
while (Nboxes_put < Nboxes && working_sector < 2*base.Nrap.size() + 2) {
Nboxes_put += Tableau[working_sector].Add_Boxes_From_Lowest(Nboxes - Nboxes_put);
working_sector += 2;
}
Compute_id();
if (Nboxes_put < Nboxes) (*this) = offsets_init; // reset !
else if (Nboxes_put > Nboxes) {
cout << Nboxes_put << "\t" << Nboxes << endl;
ABACUSerror("Putting too many boxes in Ix2_Offsets::Add_Boxes_From_Lowest.");
}
return(Nboxes_put == Nboxes);
}
//***************************************************************************************************
// Function definitions: class ODSLF_Ix2_Offsets_List
ODSLF_Ix2_Offsets_List::ODSLF_Ix2_Offsets_List() : ndef(0), Offsets(Vect<ODSLF_Ix2_Offsets>(NBASESMAX)) {};
ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets_List::Return_Offsets (ODSLF_Base& RefBase, Vect<int> nparticles)
{
long long int Ix2_Offsets_type_id_ref = Ix2_Offsets_type_id (nparticles);
int n = 0;
while (n < ndef && !(Ix2_Offsets_type_id_ref == Offsets[n].type_id && RefBase == Offsets[n].base)) n++;
if (n == ndef) {
Offsets[n] = ODSLF_Ix2_Offsets (RefBase, nparticles);
ndef++;
}
if (ndef >= NBASESMAX) ABACUSerror("Ix2_Offsets_List: need bigger Offsets vector.");
return(Offsets[n]);
}
ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets_List::Return_Offsets (ODSLF_Base& RefBase, long long int req_type_id)
{
int n = 0;
while (n < ndef && !(req_type_id == Offsets[n].type_id && RefBase == Offsets[n].base)) n++;
if (n == ndef) {
Offsets[n] = ODSLF_Ix2_Offsets (RefBase, req_type_id);
ndef++;
}
if (ndef >= NBASESMAX) ABACUSerror("Ix2_Offsets_List: need bigger Offsets vector.");
return(Offsets[n]);
}
//***************************************************************************************************
// Function definitions: class ODSLF_Bethe_State
ODSLF_Bethe_State::ODSLF_Bethe_State ()
: chain(Heis_Chain()), base(ODSLF_Base()), offsets(ODSLF_Ix2_Offsets()), Ix2(ODSLF_Ix2_Config(chain, 1)),
lambda(ODSLF_Lambda(chain, 1)), BE(ODSLF_Lambda(chain, 1)), diffsq(0.0), conv(0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0), base_id(0LL), type_id(0LL), id(0LL), maxid(0LL), nparticles(0)
{
};
ODSLF_Bethe_State::ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState) // copy constructor
: chain(RefState.chain), base(RefState.base), offsets(RefState.offsets),
Ix2(ODSLF_Ix2_Config(RefState.chain, RefState.base)),
lambda(ODSLF_Lambda(RefState.chain, RefState.base)), BE(ODSLF_Lambda(RefState.chain, RefState.base)),
diffsq(RefState.diffsq), conv(RefState.conv), iter(RefState.iter), iter_Newton(RefState.iter_Newton),
E(RefState.E), iK(RefState.iK), K(RefState.K), lnnorm(RefState.lnnorm),
base_id(RefState.base_id), type_id(RefState.type_id), id(RefState.id), maxid(RefState.maxid),
nparticles(RefState.nparticles)
{
// copy arrays into new ones
for (int j = 0; j < RefState.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < RefState.base[j]; ++j) {
Ix2[j][alpha] = RefState.Ix2[j][alpha];
lambda[j][alpha] = RefState.lambda[j][alpha];
}
}
}
ODSLF_Bethe_State::ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState, long long int type_id_ref)
: chain(RefState.chain), base(RefState.base), offsets(ODSLF_Ix2_Offsets(RefState.base, type_id_ref)),
Ix2(ODSLF_Ix2_Config(RefState.chain, RefState.base)), lambda(ODSLF_Lambda(RefState.chain, RefState.base)),
BE(ODSLF_Lambda(RefState.chain, RefState.base)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0),
base_id(RefState.base_id), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
{
(*this).Set_Ix2_Offsets (offsets);
}
ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain, int M)
: chain(RefChain), base(RefChain, M), offsets(base, 0LL), Ix2(ODSLF_Ix2_Config(RefChain, M)),
lambda(ODSLF_Lambda(RefChain, M)),
BE(ODSLF_Lambda(RefChain, M)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0),
base_id(0LL), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
{
}
ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain, const ODSLF_Base& RefBase)
: chain(RefChain), base(RefBase), offsets(RefBase, 0LL),
Ix2(ODSLF_Ix2_Config(RefChain, RefBase)), lambda(ODSLF_Lambda(RefChain, RefBase)),
BE(ODSLF_Lambda(RefChain, RefBase)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0),
base_id(RefBase.id), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
{
// Check that the number of rapidities is consistent with Mdown
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += base[i] * RefChain.Str_L[i];
if (RefBase.Mdown != Mcheck)
ABACUSerror("Wrong M from Nrapidities input vector, in ODSLF_Bethe_State constructor.");
}
ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain,
long long int base_id_ref, long long int type_id_ref)
: chain(RefChain), base(ODSLF_Base(RefChain, base_id_ref)), offsets(base, type_id_ref),
Ix2(ODSLF_Ix2_Config(RefChain, base)), lambda(ODSLF_Lambda(RefChain, base)),
BE(ODSLF_Lambda(RefChain, base)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0),
base_id(base_id_ref), type_id(type_id_ref), id(0LL), maxid(offsets.maxid), nparticles(0)
{
}
void ODSLF_Bethe_State::Set_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset) // sets the Ix2 to given offsets
{
if (base != RefOffset.base)
ABACUSerror("Offset with incompatible base used in ODSLF_Bethe_State::Set_Ix2_Offsets.");
// For each base_level, we make sure that the Ix2 are ordered: Ix2[0][j] < Ix2[0][k] if j < k.
// Level 2: particles in R part outside GS interval
for (int j = 0; j < RefOffset.Tableau[2].Nrows; ++j) {
Ix2[0][base.Nrap[0] - 1 - j]
= base.Nrap[0] - 1 + 2* RefOffset.Tableau[2].Nrows - 2*j + 2*RefOffset.Tableau[2].Row_L[j];
}
nparticles = RefOffset.Tableau[2].Nrows;
// Level 0: holes in R part of GS interval
for (int j = 0; j < RefOffset.Tableau[0].Ncols; ++j) {
Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[2].Nrows - j] = (base.Nrap[0] + 1) % 2
+ 2*RefOffset.Tableau[0].Ncols - 2 - 2*j + 2*RefOffset.Tableau[0].Col_L[j];
}
nparticles += RefOffset.Tableau[0].Ncols;
// Level 1: holes in L part of GS interval
for (int j = 0; j < RefOffset.Tableau[1].Ncols; ++j) {
Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[0].Ncols - RefOffset.Tableau[2].Nrows - j]
= -1 - (base.Nrap[0] % 2) - 2* j - 2*RefOffset.Tableau[1].Col_L[RefOffset.Tableau[1].Ncols - 1 - j];
}
nparticles += RefOffset.Tableau[1].Ncols;
// Level 3: particles in L part outside GS interval
for (int j = 0; j < RefOffset.Tableau[3].Nrows; ++j) {
Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[0].Ncols - RefOffset.Tableau[1].Ncols - RefOffset.Tableau[2].Nrows - j]
= -(base.Nrap[0] + 1 + 2*j + 2*RefOffset.Tableau[3].Row_L[RefOffset.Tableau[3].Nrows - 1 - j]);
}
nparticles += RefOffset.Tableau[3].Nrows;
// Subsequent levels: particles on R and L
for (int base_level = 1; base_level < base.Nrap.size(); ++base_level) {
for (int j = 0; j < RefOffset.Tableau[2*base_level + 2].Nrows; ++j)
Ix2[base_level][base.Nrap[base_level] - 1 - j]
= ((base.Nrap[base_level] + 1) % 2) + 2*(RefOffset.Tableau[2*base_level + 2].Nrows - 1) -2*j
+ 2*RefOffset.Tableau[2*base_level + 2].Row_L[j];
nparticles += RefOffset.Tableau[2*base_level + 2].Nrows;
for (int j = 0; j < RefOffset.Tableau[2*base_level + 3].Nrows; ++j)
Ix2[base_level][base.Nrap[base_level] - 1 - RefOffset.Tableau[2*base_level + 2].Nrows - j]
= -1 - ((base.Nrap[base_level]) % 2) -2*j
- 2*RefOffset.Tableau[2*base_level + 3].Row_L[RefOffset.Tableau[2*base_level + 3].Nrows - 1 - j];
nparticles += RefOffset.Tableau[2*base_level + 3].Nrows;
}
// ODSLF: put back correct parity of quantum nr:
for (int j = 0; j < base.Nrap.size(); ++j)
for (int alpha = 0; alpha < base[j]; ++alpha)
if ((Ix2[j][alpha] + chain.Nsites - 1 + base.Nrap[j] + (base.Nraptot + 1) * chain.Str_L[j]
+ (chain.par[j] == 1 ? 0 : chain.Nsites)) % 2)
Ix2[j][alpha] -= 1;
type_id = RefOffset.type_id;
id = RefOffset.id;
maxid = RefOffset.maxid;
return;
}
void ODSLF_Bethe_State::Set_to_id (long long int id_ref)
{
offsets.Set_to_id (id_ref);
(*this).Set_Ix2_Offsets (offsets);
}
void ODSLF_Bethe_State::Set_to_id (long long int id_ref, ODSLF_Bethe_State& RefState)
{
offsets.Set_to_id (id_ref);
(*this).Set_Ix2_Offsets (offsets);
}
bool ODSLF_Bethe_State::Check_Symmetry ()
{
// Checks whether the I's are symmetrically distributed.
bool symmetric_state = true;
int arg, test1, test3;
for (int j = 0; j < chain.Nstrings; ++j) {
test1 = 0;
test3 = 0;
for (int alpha = 0; alpha < base[j]; ++alpha) {
// since Ix2 = N is same as Ix2 = -N by periodicity, this is symmetric.
arg = (Ix2[j][alpha] != chain.Nsites) ? Ix2[j][alpha] : 0;
test1 += arg;
test3 += arg * arg * arg; // to make sure that all I's are symmetrical...
}
if (!(symmetric_state && (test1 == 0) && (test3 == 0))) symmetric_state = false;
}
return symmetric_state;
}
void ODSLF_Bethe_State::Compute_diffsq ()
{
DP maxterm = 0.0;
int jmax, alphamax;
diffsq = 0.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
diffsq += pow(BE[j][alpha], 2.0);
if (pow(BE[j][alpha], 2.0)/chain.Nsites > maxterm) {
jmax = j;
alphamax = alpha;
maxterm = pow(BE[j][alpha], 2.0)/chain.Nsites;
}
}
diffsq /= DP(chain.Nsites);
}
void ODSLF_Bethe_State::Iterate_BAE ()
{
for (int j = 0; j < chain.Nstrings; ++j) {
for (int alpha = 0; alpha < base[j]; ++alpha)
{
lambda[j][alpha] = Iterate_BAE (j, alpha);
}
}
iter++;
(*this).Compute_BE();
(*this).Compute_diffsq();
}
void ODSLF_Bethe_State::BAE_smackdown (DP max_allowed)
{
// Re-solves for all rapidities lambda[j][alpha] such that BE[j][alpha]^2/N > max_allowed.
// Assumes that BE[][] is up-to-date.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha)
if (pow(BE[j][alpha], 2.0)/chain.Nsites > max_allowed) (*this).Solve_BAE (j, alpha, max_allowed, 100);
}
void ODSLF_Bethe_State::Solve_BAE_smackdown (DP max_allowed, int maxruns)
{
int runs_done = 0;
(*this).Compute_BE();
(*this).Compute_diffsq();
while (diffsq > chain.prec && diffsq > max_allowed && runs_done < maxruns) {
(*this).BAE_smackdown (max_allowed);
(*this).Compute_BE();
(*this).Compute_diffsq();
runs_done++;
}
}
void ODSLF_Bethe_State::Iterate_BAE_Newton ()
{
// does one step of a Newton method on the rapidities...
// Assumes that BE[j][alpha] have been computed
Vect_CX dlambda (0.0, base.Nraptot); // contains delta lambda computed from Newton's method
SQMat_CX Gaudin (0.0, base.Nraptot);
Vect_INT indx (base.Nraptot);
ODSLF_Lambda lambda_ref(chain, base);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
(*this).Build_Reduced_Gaudin_Matrix (Gaudin);
int index = 0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
dlambda[index] = - BE[j][alpha] * chain.Nsites;
index++;
}
DP d;
try {
ludcmp_CX (Gaudin, indx, d);
lubksb_CX (Gaudin, indx, dlambda);
}
catch (Divide_by_zero) {
diffsq = log(-1.0); // reset to nan to stop Newton method
return;
}
// Regularize dlambda: max step is +-1.0 to prevent rapidity flying off into outer space.
for (int i = 0; i < base.Nraptot; ++i) if (fabs(real(dlambda[i])) > 1.0) dlambda[i] = 0.0;
index = 0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
lambda[j][alpha] = lambda_ref[j][alpha] + real(dlambda[index]);
index++;
}
(*this).Compute_BE();
(*this).Iterate_BAE();
(*this).Compute_diffsq();
// if we've converged, calculate the norm here, since the work has been done...
if (diffsq < chain.prec) {
lnnorm = 0.0;
for (int j = 0; j < base.Nraptot; j++) lnnorm += log(abs(Gaudin[j][j]));
}
iter_Newton++;
return;
}
void ODSLF_Bethe_State::Solve_BAE (int j, int alpha, DP req_prec, int itermax)
{
// Finds the root lambda[j][alpha] to precision req_prec
DP prec_obtained = 1.0;
int niter_here = 0;
while (prec_obtained > req_prec && niter_here < itermax) {
lambda[j][alpha] = (*this).Iterate_BAE (j, alpha);
(*this).Compute_BE (j, alpha);
prec_obtained = pow(BE[j][alpha], 2.0)/chain.Nsites;
niter_here++;
}
}
void ODSLF_Bethe_State::Find_Rapidities (bool reset_rapidities)
{
// This function finds the rapidities of the eigenstate
lnnorm = -100.0; // sentinel value, recalculated if Newton method used in the last step of iteration.
ODSLF_Lambda lambda_ref(chain, base);
diffsq = 1.0;
if (reset_rapidities)
(*this).Set_Free_lambdas();
(*this).Compute_BE();
iter = 0;
iter_Newton = 0;
// Start with conventional iterations
DP iter_prec = 1.0e-2/DP(chain.Nsites);
clock_t interp_start;
clock_t interp_stop;
clock_t Newton_start;
clock_t Newton_stop;
int iter_interp_start, iter_interp_stop, iter_Newton_start, iter_Newton_stop;
DP diffsq_interp_start, diffsq_interp_stop, diffsq_Newton_start, diffsq_Newton_stop;
while (diffsq > chain.prec && iter < 300 && iter_Newton < 20) {
do {
interp_start = clock();
iter_interp_start = iter;
diffsq_interp_start = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_interp (iter_prec, 10);
interp_stop = clock();
iter_interp_stop = iter;
diffsq_interp_stop = diffsq;
iter_prec = diffsq * 0.001;
} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_interp_stop/diffsq_interp_start < 0.001);
// If we haven't even reached iter_prec, try normal iterations without extrapolations
if (diffsq > iter_prec) {
(*this).Solve_BAE_smackdown (0.1 * diffsq, 1);
}
// Now try Newton
Newton_start = clock();
iter_Newton_start = iter_Newton;
diffsq_Newton_start = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_Newton (chain.prec, 5);
Newton_stop = clock();
iter_Newton_stop = iter_Newton;
diffsq_Newton_stop = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_smackdown (0.1 * diffsq, 1);
if (diffsq > iter_prec) (*this).Solve_BAE_straight_iter (0.1 * diffsq, 50);
// If we haven't converged here, retry with more precision in interp method...
if (is_nan(diffsq)) {
(*this).Set_Free_lambdas(); // good start if we've messed up with Newton
diffsq = 1.0;
}
iter_prec *= 1.0e-4;
iter_prec = ABACUS::max(iter_prec, chain.prec);
}
// Check convergence:
conv = (diffsq < chain.prec && (*this).Check_Rapidities()) ? 1 : 0;
return;
}
void ODSLF_Bethe_State::Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton)
{
// This function attempts to get convergence diffsq <= Newton_prec in at most max_iter_Newton steps.
// The results are accepted if diffsq has decreased, otherwise the lambda's are reset to original values, defined as...
ODSLF_Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
int iter_done_here = 0;
while ((diffsq > Newton_prec) && (diffsq < 10.0) && (!is_nan(diffsq)) && (iter_done_here < max_iter_Newton)) {
(*this).Iterate_BAE_Newton();
iter_done_here++;
(*this).Iterate_BAE();
}
if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
diffsq = diffsq_ref;
}
return;
}
void ODSLF_Bethe_State::Solve_BAE_straight_iter (DP interp_prec, int max_iter_interp)
{
// This function attempts to get convergence diffsq <= interp_prec in at most max_iter_interp steps.
// If this fails, the lambda's are reset to original values, defined as...
ODSLF_Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
diffsq = 1.0;
int iter_done_here = 0;
while ((diffsq > interp_prec) && (max_iter_interp > iter_done_here)) {
(*this).Iterate_BAE();
iter_done_here++;
}
if ((diffsq > diffsq_ref)) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
diffsq = diffsq_ref;
}
}
void ODSLF_Bethe_State::Solve_BAE_interp (DP interp_prec, int max_iter_interp)
{
// This function attempts to get convergence diffsq <= interp_prec in at most max_iter_interp steps.
// If this fails, the lambda's are reset to original values, defined as...
ODSLF_Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
diffsq = 1.0;
int iter_done_here = 0;
ODSLF_Lambda lambda1(chain, base);
ODSLF_Lambda lambda2(chain, base);
ODSLF_Lambda lambda3(chain, base);
ODSLF_Lambda lambda4(chain, base);
ODSLF_Lambda lambda5(chain, base);
while ((diffsq > interp_prec) && (max_iter_interp > iter_done_here)) {
(*this).Iterate_BAE();
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda1[j][alpha] = lambda[j][alpha];
(*this).Iterate_BAE();
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda2[j][alpha] = lambda[j][alpha];
(*this).Iterate_BAE();
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda3[j][alpha] = lambda[j][alpha];
(*this).Iterate_BAE();
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda4[j][alpha] = lambda[j][alpha];
iter_done_here += 4;
// now extrapolate the result to infinite number of iterations...
Vect_DP rap(0.0, 4);
Vect_DP oneoverP(0.0, 4);
DP deltalambda = 0.0;
for (int i = 0; i < 4; ++i) oneoverP[i] = 1.0/(1.0 + i*i);
for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) {
rap[0] = lambda1[j][alpha];
rap[1] = lambda2[j][alpha];
rap[2] = lambda3[j][alpha];
rap[3] = lambda4[j][alpha];
polint (oneoverP, rap, 0.0, lambda[j][alpha], deltalambda);
}
(*this).Iterate_BAE();
}
if ((diffsq >= diffsq_ref)) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
diffsq = diffsq_ref;
}
return;
}
void ODSLF_Bethe_State::Compute_lnnorm ()
{
if (true || lnnorm == -100.0) { // else norm already calculated by Newton method
// Actually, compute anyway to increase accuracy !
SQMat_CX Gaudin_Red(base.Nraptot);
(*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
complex<DP> lnnorm_check = lndet_LU_CX_dstry(Gaudin_Red);
lnnorm = real(lnnorm_check);
}
return;
}
void ODSLF_Bethe_State::Compute_Momentum ()
{
int sum_Ix2 = 0;
DP sum_M = 0.0;
for (int j = 0; j < chain.Nstrings; ++j) {
sum_M += 0.5 * (1.0 + chain.par[j]) * base[j];
for (int alpha = 0; alpha < base[j]; ++alpha) {
sum_Ix2 += Ix2[j][alpha];
}
}
iK = (chain.Nsites/2) * int(sum_M + 0.1) - (sum_Ix2/2); // + 0.1: for safety...
while (iK >= chain.Nsites) iK -= chain.Nsites;
while (iK < 0) iK += chain.Nsites;
K = PI * sum_M - PI * sum_Ix2/chain.Nsites;
while (K >= 2.0*PI) K -= 2.0*PI;
while (K < 0.0) K += 2.0*PI;
return;
}
void ODSLF_Bethe_State::Compute_All (bool reset_rapidities) // solves BAE, computes E, K and lnnorm
{
(*this).Find_Rapidities (reset_rapidities);
if (conv == 1) {
(*this).Compute_Energy ();
(*this).Compute_Momentum ();
(*this).Compute_lnnorm ();
}
return;
}
bool ODSLF_Bethe_State::Boost_Momentum (int iKboost)
{
if (iKboost == 0) return(true); // done
ODSLF_Ix2_Offsets offsets_here = offsets;
bool success = false;
if (iKboost < 0)
success = offsets_here.Add_Boxes_From_Lowest(-iKboost, 0); // add boxes in even sectors to decrease iK
else if (iKboost > 0)
success = offsets_here.Add_Boxes_From_Lowest(iKboost, 1); // add boxes in odd sectors to increase iK
if (success) (*this).Set_Ix2_Offsets(offsets_here);
return(success);
}
std::ostream& operator<< (std::ostream& s, const ODSLF_Bethe_State& state)
{
// sends all the state data to output stream
s << endl << "******** Chain with Nsites = " << state.chain.Nsites
<< " Mdown (nr fermions) = " << state.base.Mdown
<< ": eigenstate with base_id " << state.base_id << ", type_id " << state.type_id
<< " id " << state.id << " maxid " << state.offsets.maxid << endl
<< "E = " << state.E << " K = " << state.K << " iK = " << state.iK
<< " lnnorm = " << state.lnnorm << endl
<< "conv = " << state.conv << " iter = " << state.iter
<< " iter_Newton = " << state.iter_Newton << "\tdiffsq " << state.diffsq << endl;
for (int j = 0; j < state.chain.Nstrings; ++j) {
if (state.base.Nrap[j] > 0) {
s << "Type " << j << " Str_L = " << state.chain.Str_L[j]
<< " par = " << state.chain.par[j] << " M_j = " << state.base.Nrap[j]
<< " Ix2_infty = " << state.base.Ix2_infty[j] << " Ix2_max = " << state.base.Ix2_max[j] << endl;
Vect_INT qnumbers(state.base.Nrap[j]);
Vect_DP rapidities(state.base.Nrap[j]);
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) {
qnumbers[alpha] = state.Ix2[j][alpha];
rapidities[alpha] = state.lambda[j][alpha];
}
qnumbers.QuickSort();
rapidities.QuickSort();
s << "Ix2 quantum numbers: " << endl;
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << qnumbers[alpha] << " ";
s << endl;
s << "Rapidities: " << endl;
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << rapidities[alpha] << " ";
s << endl;
}
}
s << endl;
return s;
}
} // namespace ABACUS
In neuem Issue referenzieren Git Blame ansehen Permalink kopieren