19 #ifdef USE_PARAMETERS_FACTORY
38 #ifdef USE_PARAMETERS_FACTORY
53 const string str_vlevel = params.
get_string(
"verbose_level");
58 string str_sortfield_type;
61 double Enorm_eigen, Vthreshold;
64 err += params.
fetch_string(
"eigensolver_mode", str_sortfield_type);
65 err += params.
fetch_int(
"number_of_wanted_eigenvectors", Nk);
66 err += params.
fetch_int(
"number_of_working_eigenvectors", Np);
67 err += params.
fetch_int(
"maximum_number_of_iteration", Niter_eigen);
68 err += params.
fetch_double(
"convergence_criterion_squared", Enorm_eigen);
69 err += params.
fetch_double(
"threshold_value", Vthreshold);
83 int Niter_eigen,
double Enorm_eigen,
118 int& Nsbt,
int& Nconv,
const Field& b)
132 if (Nk + Np > TDa.size()) {
135 }
else if (Nk + Np > vk.size()) {
140 valarray<double> TDb(Nm);
141 valarray<double> TDa2(Nm);
142 valarray<double> TDb2(Nm);
143 valarray<double> Qt(Nm * Nm);
144 valarray<int> Iconv(Nm);
146 int Nin = vk[0].nin();
147 int Nvol = vk[0].nvol();
148 int Nex = vk[0].nex();
149 valarray<Field> B(Nm);
150 for (
int k = 0; k < Nm; ++k) {
151 B[k].reset(Nin, Nvol, Nex);
172 double vnorm = vk[0] * vk[0];
173 vk[0] = 1.0 / sqrt(vnorm);
177 for (
int k = 0; k < k2; ++k) {
178 step(Nm, k, TDa, TDb, vk, f);
182 for (
int iter = 0; iter < Niter_eigen; ++iter) {
185 int Nm2 = Nm - kconv;
187 for (
int k = k2; k < Nm; ++k) {
188 step(Nm, k, TDa, TDb, vk, f);
194 for (
int k = 0; k < Nm2; ++k) {
195 TDa2[k] = TDa[k + k1 - 1];
196 TDb2[k] = TDb[k + k1 - 1];
199 tqri(TDa2, TDb2, Nm2, Nm, Qt);
206 for (
int ip = k2; ip < Nm; ++ip) {
207 double Dsh = TDa2[ip - kconv];
210 qrtrf(TDa, TDb, Nm, Nm, Qt, Dsh, kmin, kmax);
213 for (
int i = 0; i < (Nk + 1); ++i) {
217 for (
int j = k1 - 1; j < k2 + 1; ++j) {
218 for (
int k = 0; k < Nm; ++k) {
219 B[j] += Qt[k + Nm * j] * vk[k];
223 for (
int j = k1 - 1; j < k2 + 1; ++j) {
228 f *= Qt[Nm - 1 + Nm * (k2 - 1)];
229 f += TDb[k2 - 1] * vk[k2];
231 beta_k = sqrt(beta_k);
236 double beta_r = 1.0 / beta_k;
238 TDb[k2 - 1] = beta_k;
245 tqri(TDa2, TDb2, Nk, Nm, Qt);
246 for (
int k = 0; k < Nk; ++k) {
250 for (
int j = 0; j < Nk; ++j) {
251 for (
int k = 0; k < Nk; ++k) {
252 B[j] += Qt[k + j * Nm] * vk[k];
259 for (
int i = 0; i < Nk; ++i) {
261 double vnum = B[i] * v;
262 double vden = B[i] * B[i];
266 TDa2[i] = vnum / vden;
272 if (vv < Enorm_eigen) {
286 if (Kthreshold > 0) {
291 for (
int i = 0; i < Kdis; ++i) {
292 TDa[i] = TDa2[Iconv[i]];
298 Nsbt = Kdis - Kthreshold;
320 valarray<double>& TDb, valarray<Field>& vk,
328 double alph = vk[k] * w;
333 double beta_r = 1.0 / beta;
334 vk[k + 1] = beta_r * w;
340 w -= TDb[k - 1] * vk[k - 1];
341 double alph = vk[k] * w;
346 double beta_r = 1.0 / beta;
354 if (k < Nm - 1) vk[k + 1] = w;
361 valarray<Field>& vk,
int k)
363 for (
int j = 0; j < k; ++j) {
364 dcomplex prod =
dotc(vk[j], w);
367 prod = cmplx(-1.0, 0.0) * prod;
368 axpy(w, prod, vk[j]);
376 for (
int i = 0; i < Qt.size(); ++i) {
380 for (
int k = 0; k < Nm; ++k) {
381 Qt[k + k * Nm] = 1.0;
388 valarray<double>& TDb,
389 int Nk,
int Nm, valarray<double>& Qt)
391 int Niter = 100 * Nm;
399 for (
int iter = 0; iter < Niter; ++iter) {
401 double dsub = TDa[kmax - 1] - TDa[kmax - 2];
402 double dd = sqrt(dsub * dsub + 4.0 * TDb[kmax - 2] * TDb[kmax - 2]);
403 double Dsh = 0.5 * (TDa[kmax - 2] + TDa[kmax - 1]
404 + fabs(dd) * (dsub / fabs(dsub)));
408 qrtrf(TDa, TDb, Nk, Nm, Qt, Dsh, kmin, kmax);
411 for (
int j = kmax - 1; j >= kmin; --j) {
412 double dds = fabs(TDa[j - 1]) + fabs(TDa[j]);
413 if (fabs(TDb[j - 1]) + dds > dds) {
416 for (
int j = 0; j < kmax - 1; ++j) {
417 double dds = fabs(TDa[j]) + fabs(TDa[j + 1]);
419 if (fabs(TDb[j]) + dds > dds) {
447 valarray<double>& TDb,
448 int Nk,
int Nm, valarray<double>& Qt,
449 double Dsh,
int kmin,
int kmax)
454 double Fden = 1.0 / sqrt((TDa[k] - Dsh) * (TDa[k] - Dsh)
456 double c = (TDa[k] - Dsh) * Fden;
457 double s = -TDb[k] * Fden;
459 double tmpa1 = TDa[k];
460 double tmpa2 = TDa[k + 1];
461 double tmpb = TDb[k];
463 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
464 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
465 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
467 TDb[k + 1] = c * TDb[k + 1];
469 for (
int i = 0; i < Nk; ++i) {
470 double Qtmp1 = Qt[i + Nm * k];
471 double Qtmp2 = Qt[i + Nm * (k + 1)];
472 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
473 Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
477 for (
int k = kmin; k < kmax - 1; ++k) {
478 double Fden = 1.0 / sqrt(x * x + TDb[k - 1] * TDb[k - 1]);
479 double c = TDb[k - 1] * Fden;
480 double s = -x * Fden;
482 double tmpa1 = TDa[k];
483 double tmpa2 = TDa[k + 1];
484 double tmpb = TDb[k];
485 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
486 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
487 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
488 TDb[k - 1] = c * TDb[k - 1] - s * x;
491 TDb[k + 1] = c * TDb[k + 1];
494 for (
int i = 0; i < Nk; ++i) {
495 double Qtmp1 = Qt[i + Nm * k];
496 double Qtmp2 = Qt[i + Nm * (k + 1)];
497 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
498 Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
Bridge::VerboseLevel m_vl
void schmidt_orthogonalization(Field &w, std::valarray< Field > &vk, int k)
void Register_string(const string &, const string &)
virtual const Field mult(const Field &)=0
multiplies fermion operator to a given field and returns the resultant field.
virtual void sort(int, std::valarray< double > &)=0
void general(const char *format,...)
void Register_int(const string &, const int)
Container of Field-type object.
static const std::string class_name
void setUnit_Qt(int Nm, std::valarray< double > &Qt)
int square_non_zero(const double v)
Parameters_Eigensolver_IRLanczos()
void set_parameters(const Parameters ¶ms)
dcomplex dotc(const Field &y, const Field &x)
int fetch_string(const string &key, string &val) const
void paranoiac(const char *format,...)
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
void crucial(const char *format,...)
void solve(std::valarray< double > &TDa, std::valarray< Field > &vk, int &Nsbt, int &Nconv, const Field &b)
static bool Register(const std::string &realm, const creator_callback &cb)
void tqri(std::valarray< double > &TDa, std::valarray< double > &TDb, int Nk, int Nm, std::valarray< double > &Qt)
virtual int converged(double v, double v_thrs)=0
int non_negative(const int v)
void Register_double(const string &, const double)
int fetch_double(const string &key, double &val) const
string get_string(const string &key) const
void qrtrf(std::valarray< double > &TDa, std::valarray< double > &TDb, int Nk, int Nm, std::valarray< double > &Qt, double Dsh, int kmin, int kmax)
int fetch_int(const string &key, int &val) const
static VerboseLevel set_verbose_level(const std::string &str)
void step(int Nm, int k, std::valarray< double > &TDa, std::valarray< double > &TDb, std::valarray< Field > &vk, Field &f)