28 const std::string str_vlevel = params.
get_string(
"verbose_level");
33 std::string str_sortfield_type;
36 double Enorm_eigen, Vthreshold;
39 err += params.
fetch_string(
"eigensolver_mode", str_sortfield_type);
40 err += params.
fetch_int(
"number_of_wanted_eigenvectors", Nk);
41 err += params.
fetch_int(
"number_of_working_eigenvectors", Np);
42 err += params.
fetch_int(
"maximum_number_of_iteration", Niter_eigen);
43 err += params.
fetch_double(
"convergence_criterion_squared", Enorm_eigen);
44 err += params.
fetch_double(
"threshold_value", Vthreshold);
51 set_parameters(str_sortfield_type, Nk, Np, Niter_eigen, Enorm_eigen, Vthreshold);
57 const int Nk,
const int Np,
58 const int Niter_eigen,
const double Enorm_eigen,
59 const double Vthreshold)
71 const int Niter_eigen,
double Enorm_eigen,
106 int& Nsbt,
int& Nconv,
const Field& b)
117 const int Nm = Nk + Np;
119 if (Nk + Np > TDa.size()) {
122 }
else if (Nk + Np > vk.size()) {
127 std::vector<double> TDb(Nm);
128 std::vector<double> TDa2(Nm);
129 std::vector<double> TDb2(Nm);
130 std::vector<double> Qt(Nm * Nm);
131 std::vector<int> Iconv(Nm);
133 const int Nin = vk[0].nin();
134 const int Nvol = vk[0].nvol();
135 const int Nex = vk[0].nex();
136 std::vector<Field> B(Nm);
137 for (
int k = 0; k < Nm; ++k) {
138 B[k].reset(Nin, Nvol, Nex);
161 const double vnorm =
dot(vk[0], vk[0]);
162 vk[0].set(1.0 / sqrt(vnorm));
166 for (
int k = 0; k < k2; ++k) {
167 step(Nm, k, TDa, TDb, vk, f);
171 for (
int iter = 0; iter < Niter_eigen; ++iter) {
174 int Nm2 = Nm - kconv;
176 for (
int k = k2; k < Nm; ++k) {
177 step(Nm, k, TDa, TDb, vk, f);
181 scal(f, TDb[Nm - 1]);
184 for (
int k = 0; k < Nm2; ++k) {
185 TDa2[k] = TDa[k + k1 - 1];
186 TDb2[k] = TDb[k + k1 - 1];
189 tqri(TDa2, TDb2, Nm2, Nm, Qt);
196 for (
int ip = k2; ip < Nm; ++ip) {
197 double Dsh = TDa2[ip - kconv];
200 qrtrf(TDa, TDb, Nm, Nm, Qt, Dsh, kmin, kmax);
203 for (
int i = 0; i < (Nk + 1); ++i) {
207 for (
int j = k1 - 1; j < k2 + 1; ++j) {
208 for (
int k = 0; k < Nm; ++k) {
209 axpy(B[j], Qt[k + Nm * j], vk[k]);
213 for (
int j = k1 - 1; j < k2 + 1; ++j) {
218 scal(f, Qt[Nm - 1 + Nm * (k2 - 1)]);
219 axpy(f, TDb[k2 - 1], vk[k2]);
221 double beta_k =
dot(f, f);
222 beta_k = sqrt(beta_k);
227 const double beta_r = 1.0 / beta_k;
229 scal(vk[k2], beta_r);
230 TDb[k2 - 1] = beta_k;
237 tqri(TDa2, TDb2, Nk, Nm, Qt);
238 for (
int k = 0; k < Nk; ++k) {
242 for (
int j = 0; j < Nk; ++j) {
243 for (
int k = 0; k < Nk; ++k) {
244 axpy(B[j], Qt[k + j * Nm], vk[k]);
251 for (
int i = 0; i < Nk; ++i) {
253 const double vnum =
dot(B[i], v);
254 const double vden =
dot(B[i], B[i]);
258 TDa2[i] = vnum / vden;
259 axpy(v, -TDa2[i], B[i]);
261 const double vv =
dot(v, v);
265 if (vv < Enorm_eigen) {
279 if (Kthreshold > 0) {
284 for (
int i = 0; i < Kdis; ++i) {
285 TDa[i] = TDa2[Iconv[i]];
290 for (
int i = 0; i < Kdis; ++i) {
291 vk[i] = B[Iconv[idx[i]]];
294 Nsbt = Kdis - Kthreshold;
316 std::vector<double>& TDa, std::vector<double>& TDb,
317 std::vector<Field>& vk,
Field& w)
324 const double alph =
dot(vk[k], w);
326 axpy(w, -alph, vk[k]);
328 double beta =
dot(w, w);
331 const double beta_r = 1.0 / beta;
333 scal(vk[k + 1], beta_r);
339 axpy(w, -TDb[k - 1], vk[k - 1]);
341 const double alph =
dot(vk[k], w);
343 axpy(w, -alph, vk[k]);
345 double beta =
dot(w, w);
348 const double beta_r = 1.0 / beta;
356 if (k < Nm - 1) vk[k + 1] = w;
363 const std::vector<Field>& vk,
const int k)
365 for (
int j = 0; j < k; ++j) {
366 dcomplex prod =
dotc(vk[j], w);
367 prod *= cmplx(-1.0, 0.0);
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 std::vector<double>& TDb,
389 const int Nk,
const int Nm, std::vector<double>& Qt)
391 const 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 std::vector<double>& TDb,
448 const int Nk,
const int Nm, std::vector<double>& Qt,
449 const double Dsh,
const int kmin,
const int kmax)
451 const int k = kmin - 1;
453 const double Fden = 1.0 / sqrt((TDa[k] - Dsh) * (TDa[k] - Dsh)
455 const double c = (TDa[k] - Dsh) * Fden;
456 const double s = -TDb[k] * Fden;
458 const double tmpa1 = TDa[k];
459 const double tmpa2 = TDa[k + 1];
460 const double tmpb = TDb[k];
462 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
463 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
464 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
465 double x = -s * TDb[k + 1];
466 TDb[k + 1] = c * TDb[k + 1];
468 for (
int i = 0; i < Nk; ++i) {
469 double Qtmp1 = Qt[i + Nm * k];
470 double Qtmp2 = Qt[i + Nm * (k + 1)];
471 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
472 Qt[i + Nm * (k + 1)] = s * Qtmp1 + c * Qtmp2;
476 for (
int k = kmin; k < kmax - 1; ++k) {
477 double Fden = 1.0 / sqrt(x * x + TDb[k - 1] * TDb[k - 1]);
478 double c = TDb[k - 1] * Fden;
479 double s = -x * Fden;
481 double tmpa1 = TDa[k];
482 double tmpa2 = TDa[k + 1];
483 double tmpb = TDb[k];
484 TDa[k] = c * c * tmpa1 + s * s * tmpa2 - 2.0 * c * s * tmpb;
485 TDa[k + 1] = s * s * tmpa1 + c * c * tmpa2 + 2.0 * c * s * tmpb;
486 TDb[k] = c * s * (tmpa1 - tmpa2) + (c * c - s * s) * tmpb;
487 TDb[k - 1] = c * TDb[k - 1] - s * x;
490 TDb[k + 1] = c * TDb[k + 1];
493 for (
int i = 0; i < Nk; ++i) {
494 double Qtmp1 = Qt[i + Nm * k];
495 double Qtmp2 = Qt[i + Nm * (k + 1)];
496 Qt[i + Nm * k] = c * Qtmp1 - s * Qtmp2;
497 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 setUnit_Qt(const int Nm, std::vector< double > &Qt)
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)
int fetch_string(const string &key, string &value) const
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 qrtrf(std::vector< double > &TDa, std::vector< double > &TDb, const int Nk, const int Nm, std::vector< double > &Qt, const double Dsh, const int kmin, const int kmax)
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
void crucial(const char *format,...)
virtual void mult(Field &, const Field &)=0
multiplies fermion operator to a given field (2nd argument)
int non_negative(const int v)
void step(const int Nm, const int k, std::vector< double > &TDa, std::vector< double > &TDb, std::vector< Field > &vk, Field &f)
void tqri(std::vector< double > &TDa, std::vector< double > &TDb, const int Nk, const int Nm, std::vector< double > &Qt)
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)
void schmidt_orthogonalization(Field &w, const std::vector< Field > &vk, const int k)