32 const string str_vlevel = params.
get_string(
"verbose_level");
37 string str_sortfield_type;
40 double Enorm_eigen, Vthreshold;
43 err += params.
fetch_string(
"eigensolver_mode", str_sortfield_type);
44 err += params.
fetch_int(
"number_of_wanted_eigenvectors", Nk);
45 err += params.
fetch_int(
"number_of_working_eigenvectors", Np);
46 err += params.
fetch_int(
"maximum_number_of_iteration", Niter_eigen);
47 err += params.
fetch_double(
"convergence_criterion_squared", Enorm_eigen);
48 err += params.
fetch_double(
"threshold_value", Vthreshold);
55 set_parameters(str_sortfield_type, Nk, Np, Niter_eigen, Enorm_eigen, Vthreshold);
62 int Niter_eigen,
double Enorm_eigen,
75 int Niter_eigen,
double Enorm_eigen,
110 int& Nsbt,
int& Nconv,
const Field& b)
124 if (Nk + Np > TDa.size()) {
127 }
else if (Nk + Np > vk.size()) {
132 std::vector<double> TDb(Nm);
133 std::vector<double> TDa2(Nm);
134 std::vector<double> TDb2(Nm);
135 std::vector<double> Qt(Nm * Nm);
136 std::vector<int> Iconv(Nm);
138 int Nin = vk[0].nin();
139 int Nvol = vk[0].nvol();
140 int Nex = vk[0].nex();
141 std::vector<Field> B(Nm);
142 for (
int k = 0; k < Nm; ++k) {
143 B[k].reset(Nin, Nvol, Nex);
164 double vnorm =
dot(vk[0], vk[0]);
165 vk[0].set(1.0 / sqrt(vnorm));
169 for (
int k = 0; k < k2; ++k) {
170 step(Nm, k, TDa, TDb, vk, f);
174 for (
int iter = 0; iter < Niter_eigen; ++iter) {
177 int Nm2 = Nm - kconv;
179 for (
int k = k2; k < Nm; ++k) {
180 step(Nm, k, TDa, TDb, vk, f);
184 scal(f, TDb[Nm - 1]);
187 for (
int k = 0; k < Nm2; ++k) {
188 TDa2[k] = TDa[k + k1 - 1];
189 TDb2[k] = TDb[k + k1 - 1];
192 tqri(TDa2, TDb2, Nm2, Nm, Qt);
199 for (
int ip = k2; ip < Nm; ++ip) {
200 double Dsh = TDa2[ip - kconv];
203 qrtrf(TDa, TDb, Nm, Nm, Qt, Dsh, kmin, kmax);
206 for (
int i = 0; i < (Nk + 1); ++i) {
210 for (
int j = k1 - 1; j < k2 + 1; ++j) {
211 for (
int k = 0; k < Nm; ++k) {
212 axpy(B[j], Qt[k + Nm * j], vk[k]);
216 for (
int j = k1 - 1; j < k2 + 1; ++j) {
221 scal(f, Qt[Nm - 1 + Nm * (k2 - 1)]);
222 axpy(f, TDb[k2 - 1], vk[k2]);
225 beta_k = sqrt(beta_k);
230 double beta_r = 1.0 / beta_k;
232 scal(vk[k2], beta_r);
233 TDb[k2 - 1] = beta_k;
240 tqri(TDa2, TDb2, Nk, Nm, Qt);
241 for (
int k = 0; k < Nk; ++k) {
245 for (
int j = 0; j < Nk; ++j) {
246 for (
int k = 0; k < Nk; ++k) {
247 axpy(B[j], Qt[k + j * Nm], vk[k]);
254 for (
int i = 0; i < Nk; ++i) {
256 double vnum =
dot(B[i], v);
257 double vden =
dot(B[i], B[i]);
261 TDa2[i] = vnum / vden;
262 axpy(v, -TDa2[i], B[i]);
264 double vv =
dot(v, v);
268 if (vv < Enorm_eigen) {
282 if (Kthreshold > 0) {
287 for (
int i = 0; i < Kdis; ++i) {
288 TDa[i] = TDa2[Iconv[i]];
293 for (
int i = 0; i < Kdis; ++i) {
294 vk[i] = B[Iconv[idx[i]]];
297 Nsbt = Kdis - Kthreshold;
319 std::vector<double>& TDb, std::vector<Field>& vk,
327 double alph =
dot(vk[k], w);
329 axpy(w, -alph, vk[k]);
331 double beta =
dot(w, w);
333 double beta_r = 1.0 / beta;
335 scal(vk[k + 1], beta_r);
341 axpy(w, -TDb[k - 1], vk[k - 1]);
343 double alph =
dot(vk[k], w);
345 axpy(w, -alph, vk[k]);
347 double beta =
dot(w, w);
349 double beta_r = 1.0 / beta;
357 if (k < Nm - 1) vk[k + 1] = w;
364 std::vector<Field>& vk,
int k)
366 for (
int j = 0; j < k; ++j) {
367 dcomplex prod =
dotc(vk[j], w);
368 prod *= cmplx(-1.0, 0.0);
369 axpy(w, prod, vk[j]);
377 for (
int i = 0; i < Qt.size(); ++i) {
381 for (
int k = 0; k < Nm; ++k) {
382 Qt[k + k * Nm] = 1.0;
389 std::vector<double>& TDb,
390 int Nk,
int Nm, std::vector<double>& Qt)
392 int Niter = 100 * Nm;
400 for (
int iter = 0; iter < Niter; ++iter) {
402 double dsub = TDa[kmax - 1] - TDa[kmax - 2];
403 double dd = sqrt(dsub * dsub + 4.0 * TDb[kmax - 2] * TDb[kmax - 2]);
404 double Dsh = 0.5 * (TDa[kmax - 2] + TDa[kmax - 1]
405 + fabs(dd) * (dsub / fabs(dsub)));
409 qrtrf(TDa, TDb, Nk, Nm, Qt, Dsh, kmin, kmax);
412 for (
int j = kmax - 1; j >= kmin; --j) {
413 double dds = fabs(TDa[j - 1]) + fabs(TDa[j]);
414 if (fabs(TDb[j - 1]) + dds > dds) {
417 for (
int j = 0; j < kmax - 1; ++j) {
418 double dds = fabs(TDa[j]) + fabs(TDa[j + 1]);
420 if (fabs(TDb[j]) + dds > dds) {
448 std::vector<double>& TDb,
449 int Nk,
int Nm, std::vector<double>& Qt,
450 double Dsh,
int kmin,
int kmax)
455 double Fden = 1.0 / sqrt((TDa[k] - Dsh) * (TDa[k] - Dsh)
457 double c = (TDa[k] - Dsh) * Fden;
458 double s = -TDb[k] * Fden;
460 double tmpa1 = TDa[k];
461 double tmpa2 = TDa[k + 1];
462 double tmpb = TDb[k];
464 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
465 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
466 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
468 TDb[k + 1] = c * TDb[k + 1];
470 for (
int i = 0; i < Nk; ++i) {
471 double Qtmp1 = Qt[i + Nm * k];
472 double Qtmp2 = Qt[i + Nm * (k + 1)];
473 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
474 Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
478 for (
int k = kmin; k < kmax - 1; ++k) {
479 double Fden = 1.0 / sqrt(x * x + TDb[k - 1] * TDb[k - 1]);
480 double c = TDb[k - 1] * Fden;
481 double s = -x * Fden;
483 double tmpa1 = TDa[k];
484 double tmpa2 = TDa[k + 1];
485 double tmpb = TDb[k];
486 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
487 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
488 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
489 TDb[k - 1] = c * TDb[k - 1] - s * x;
492 TDb[k + 1] = c * TDb[k + 1];
495 for (
int i = 0; i < Nk; ++i) {
496 double Qtmp1 = Qt[i + Nm * k];
497 double Qtmp2 = Qt[i + Nm * (k + 1)];
498 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
499 Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
void scal(Field &x, const double a)
scal(x, a): x = a * x
Bridge::VerboseLevel m_vl
void qrtrf(std::vector< double > &TDa, std::vector< double > &TDb, int Nk, int Nm, std::vector< double > &Qt, double Dsh, int kmin, int kmax)
void detailed(const char *format,...)
double dot(const Field &y, const Field &x)
void general(const char *format,...)
void sort(std::vector< double > &v)
sort an array of values; v is sorted on exit.
Container of Field-type object.
int fetch_double(const string &key, double &value) const
void solve(std::vector< double > &TDa, std::vector< Field > &vk, int &Nsbt, int &Nconv, const Field &b)
void tqri(std::vector< double > &TDa, std::vector< double > &TDb, int Nk, int Nm, std::vector< double > &Qt)
int fetch_string(const string &key, string &value) const
void setUnit_Qt(int Nm, std::vector< double > &Qt)
static const std::string class_name
int square_non_zero(const double v)
void set_parameters(const Parameters ¶ms)
std::vector< int > sort_index(std::vector< double > &v)
sort an array and return list of index; v is sorted on exit.
dcomplex dotc(const Field &y, const Field &x)
int fetch_int(const string &key, int &value) 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 step(int Nm, int k, std::vector< double > &TDa, std::vector< double > &TDb, std::vector< Field > &vk, Field &f)
void crucial(const char *format,...)
void schmidt_orthogonalization(Field &w, std::vector< Field > &vk, int k)
virtual void mult(Field &, const Field &)=0
multiplies fermion operator to a given field (2nd argument)
int non_negative(const int v)
string get_string(const string &key) const
bool comp(const double lhs, const double rhs)
call sort condition.
static VerboseLevel set_verbose_level(const std::string &str)