Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
eigensolver_IRLanczos.cpp
Go to the documentation of this file.
1 
14 #include "eigensolver_IRLanczos.h"
15 
16 using std::valarray;
17 using std::string;
18 
19 #ifdef USE_PARAMETERS_FACTORY
20 #include "parameters_factory.h"
21 #endif
22 
23 //- parameter entries
24 namespace {
25  void append_entry(Parameters& param)
26  {
27  param.Register_string("eigensolver_mode", "NULL");
28  param.Register_int("number_of_wanted_eigenvectors", 0);
29  param.Register_int("number_of_working_eigenvectors", 0);
30  param.Register_int("maximum_number_of_iteration", 0);
31  param.Register_double("convergence_criterion_squared", 0.0);
32  param.Register_double("threshold_value", 0.0);
33 
34  param.Register_string("verbose_level", "NULL");
35  }
36 
37 
38 #ifdef USE_PARAMETERS_FACTORY
39  bool init_param = ParametersFactory::Register("Eigensolver.IRLanczos", append_entry);
40 #endif
41 }
42 //- end
43 
44 //- parameters class
46 //- end
47 
48 const std::string Eigensolver_IRLanczos::class_name = "Eigensolver_IRLanczos";
49 
50 //====================================================================
52 {
53  const string str_vlevel = params.get_string("verbose_level");
54 
55  m_vl = vout.set_verbose_level(str_vlevel);
56 
57  //- fetch and check input parameters
58  string str_sortfield_type;
59  int Nk, Np;
60  int Niter_eigen;
61  double Enorm_eigen, Vthreshold;
62 
63  int err = 0;
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);
70 
71  if (err) {
72  vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
73  abort();
74  }
75 
76 
77  set_parameters(Nk, Np, Niter_eigen, Enorm_eigen, Vthreshold);
78 }
79 
80 
81 //====================================================================
83  int Niter_eigen, double Enorm_eigen,
84  double Vthreshold)
85 {
86  //- print input parameters
87  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
88  vout.general(m_vl, " Nk = %d\n", Nk);
89  vout.general(m_vl, " Np = %d\n", Np);
90  vout.general(m_vl, " Niter_eigen = %d\n", Niter_eigen);
91  vout.general(m_vl, " Enorm_eigen = %16.8e\n", Enorm_eigen);
92  vout.general(m_vl, " Vthreshold = %16.8e\n", Vthreshold);
93 
94  //- range check
95  int err = 0;
98  err += ParameterCheck::non_negative(Niter_eigen);
99  err += ParameterCheck::square_non_zero(Enorm_eigen);
100  // NB. Vthreshold == 0 is allowed.
101 
102  if (err) {
103  vout.crucial(m_vl, "%s: parameter range check failed.\n", class_name.c_str());
104  abort();
105  }
106 
107  //- store values
108  m_Nk = Nk;
109  m_Np = Np;
110  m_Niter_eigen = Niter_eigen;
111  m_Enorm_eigen = Enorm_eigen;
112  m_Vthreshold = Vthreshold;
113 }
114 
115 
116 //====================================================================
117 void Eigensolver_IRLanczos::solve(valarray<double>& TDa, valarray<Field>& vk,
118  int& Nsbt, int& Nconv, const Field& b)
119 {
120  int Nk = m_Nk;
121  int Np = m_Np;
122  int Niter_eigen = m_Niter_eigen;
123  double Enorm_eigen = m_Enorm_eigen;
124  double Vthreshold = m_Vthreshold;
125 
126 
127  Nconv = -1;
128  Nsbt = 0;
129 
130  int Nm = Nk + Np;
131 
132  if (Nk + Np > TDa.size()) {
133  vout.crucial(m_vl, "%s: valarray TDa is too small.\n", class_name.c_str());
134  abort();
135  } else if (Nk + Np > vk.size()) {
136  vout.crucial(m_vl, "%s: valarray vk is too small.\n", class_name.c_str());
137  abort();
138  }
139 
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);
145 
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);
152  }
153 
154  Field f(vk[0]);
155  Field v(vk[0]);
156 
157  vout.general(m_vl, " Nk = %d Np = %d\n", Nk, Np);
158  vout.general(m_vl, " Nm = %d\n", Nm);
159  vout.general(m_vl, " size of TDa = %d\n", TDa.size());
160  vout.general(m_vl, " size of vk = %d\n", vk.size());
161 
162  int k1 = 1;
163  int k2 = Nk;
164  int kconv = 0;
165 
166  int Kdis = 0;
167  int Kthreshold = 0;
168  double beta_k;
169 
170  //- Set initial vector
171  vk[0] = 1.0;
172  double vnorm = vk[0] * vk[0];
173  vk[0] = 1.0 / sqrt(vnorm);
174  // (uniform vector)
175 
176  //- Initial Nk steps
177  for (int k = 0; k < k2; ++k) {
178  step(Nm, k, TDa, TDb, vk, f);
179  }
180 
181  //- Restarting loop begins
182  for (int iter = 0; iter < Niter_eigen; ++iter) {
183  vout.paranoiac(m_vl, "\n iter=%d\n", iter);
184 
185  int Nm2 = Nm - kconv;
186 
187  for (int k = k2; k < Nm; ++k) {
188  step(Nm, k, TDa, TDb, vk, f);
189  }
190 
191  f *= TDb[Nm - 1];
192 
193  //- getting eigenvalues
194  for (int k = 0; k < Nm2; ++k) {
195  TDa2[k] = TDa[k + k1 - 1];
196  TDb2[k] = TDb[k + k1 - 1];
197  }
198  setUnit_Qt(Nm, Qt);
199  tqri(TDa2, TDb2, Nm2, Nm, Qt);
200 
201  //- sorting
202  m_sort->sort(Nm, TDa2);
203 
204  //- Implicitly shifted QR transformations
205  setUnit_Qt(Nm, Qt);
206  for (int ip = k2; ip < Nm; ++ip) {
207  double Dsh = TDa2[ip - kconv];
208  int kmin = k1;
209  int kmax = Nm;
210  qrtrf(TDa, TDb, Nm, Nm, Qt, Dsh, kmin, kmax);
211  }
212 
213  for (int i = 0; i < (Nk + 1); ++i) {
214  B[i] = 0.0;
215  }
216 
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];
220  }
221  }
222 
223  for (int j = k1 - 1; j < k2 + 1; ++j) {
224  vk[j] = B[j];
225  }
226 
227  //- Compressed vector f and beta(k2)
228  f *= Qt[Nm - 1 + Nm * (k2 - 1)];
229  f += TDb[k2 - 1] * vk[k2];
230  beta_k = f * f;
231  beta_k = sqrt(beta_k);
232 
233  vout.paranoiac(m_vl, " beta(k) = %20.14f\n", beta_k);
234 
235 
236  double beta_r = 1.0 / beta_k;
237  vk[k2] = beta_r * f;
238  TDb[k2 - 1] = beta_k;
239 
240  //- Convergence test
241  TDa2 = TDa;
242  TDb2 = TDb;
243  setUnit_Qt(Nm, Qt);
244 
245  tqri(TDa2, TDb2, Nk, Nm, Qt);
246  for (int k = 0; k < Nk; ++k) {
247  B[k] = 0.0;
248  }
249 
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];
253  }
254  }
255 
256  Kdis = 0;
257  Kthreshold = 0;
258 
259  for (int i = 0; i < Nk; ++i) {
260  v = m_fopr->mult(B[i]);
261  double vnum = B[i] * v;
262  double vden = B[i] * B[i];
263 
264  // vout.paranoiac(m_vl, " vden = %20.14e\n",vden);
265 
266  TDa2[i] = vnum / vden;
267  v -= TDa2[i] * B[i];
268  double vv = v * v;
269 
270  vout.paranoiac(m_vl, " %4d %18.14f %18.14e\n", i, TDa2[i], vv);
271 
272  if (vv < Enorm_eigen) {
273  Iconv[Kdis] = i;
274  ++Kdis;
275  // if(fabs(TDa2[i]) > Vthreshold){
276  if (m_sort->converged(TDa2[i], Vthreshold)) {
277  ++Kthreshold;
278  }
279  }
280  } // i-loop end
281 
282 
283  vout.paranoiac(m_vl, " #modes converged: %d\n", Kdis);
284 
285 
286  if (Kthreshold > 0) {
287  // (there is a converged eigenvalue larger than Vthreshold.)
288  Nconv = iter;
289 
290  //- Sorting
291  for (int i = 0; i < Kdis; ++i) {
292  TDa[i] = TDa2[Iconv[i]];
293  vk[i] = B[Iconv[i]];
294  }
295 
296  m_sort->sort(Kdis, TDa, vk);
297 
298  Nsbt = Kdis - Kthreshold;
299 
300  vout.general(m_vl, "\n Converged:\n");
301  vout.general(m_vl, " Nconv = %d\n", Nconv);
302  vout.general(m_vl, " beta(k) = %20.14e\n", beta_k);
303  vout.general(m_vl, " Kdis = %d\n", Kdis);
304  vout.general(m_vl, " Nsbt = %d\n", Nsbt);
305 
306  return;
307  }
308  } // end of iter loop
309 
310 
311  if (Nconv == -1) {
312  vout.crucial(m_vl, "%s: NOT converged.\n", class_name.c_str());
313  abort();
314  }
315 }
316 
317 
318 //====================================================================
319 void Eigensolver_IRLanczos::step(int Nm, int k, valarray<double>& TDa,
320  valarray<double>& TDb, valarray<Field>& vk,
321  Field& w)
322 {
323  if (k >= Nm) {
324  vout.crucial(m_vl, "Eigensolver_IRLanczos: k is larger than Nm.\n");
325  abort();
326  } else if (k == 0) { // Initial step
327  w = m_fopr->mult(vk[k]);
328  double alph = vk[k] * w;
329 
330  w -= alph * vk[k];
331  double beta = w * w;
332  beta = sqrt(beta);
333  double beta_r = 1.0 / beta;
334  vk[k + 1] = beta_r * w;
335 
336  TDa[k] = alph;
337  TDb[k] = beta;
338  } else { // Iteration step
339  w = m_fopr->mult(vk[k]);
340  w -= TDb[k - 1] * vk[k - 1];
341  double alph = vk[k] * w;
342 
343  w -= alph * vk[k];
344  double beta = w * w;
345  beta = sqrt(beta);
346  double beta_r = 1.0 / beta;
347  w *= beta_r;
348 
349  TDa[k] = alph;
350  TDb[k] = beta;
351 
352  schmidt_orthogonalization(w, vk, k);
353 
354  if (k < Nm - 1) vk[k + 1] = w;
355  }
356 }
357 
358 
359 //====================================================================
361  valarray<Field>& vk, int k)
362 {
363  for (int j = 0; j < k; ++j) {
364  dcomplex prod = dotc(vk[j], w);
365  // w.daxpy(-prod, vk[j]);
366  // prod = - prod;
367  prod = cmplx(-1.0, 0.0) * prod;
368  axpy(w, prod, vk[j]);
369  }
370 }
371 
372 
373 //====================================================================
374 void Eigensolver_IRLanczos::setUnit_Qt(int Nm, valarray<double>& Qt)
375 {
376  for (int i = 0; i < Qt.size(); ++i) {
377  Qt[i] = 0.0;
378  }
379 
380  for (int k = 0; k < Nm; ++k) {
381  Qt[k + k * Nm] = 1.0;
382  }
383 }
384 
385 
386 //====================================================================
387 void Eigensolver_IRLanczos::tqri(valarray<double>& TDa,
388  valarray<double>& TDb,
389  int Nk, int Nm, valarray<double>& Qt)
390 {
391  int Niter = 100 * Nm;
392  int kmin = 1;
393  int kmax = Nk;
394  // (these parameters should be tuned)
395 
396 
397  int Nconv = -1;
398 
399  for (int iter = 0; iter < Niter; ++iter) {
400  //- determination of 2x2 leading submatrix
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)));
405  // (Dsh: shift)
406 
407  //- transformation
408  qrtrf(TDa, TDb, Nk, Nm, Qt, Dsh, kmin, kmax);
409 
410  //- Convergence criterion (redef of kmin and 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) {
414  kmax = j + 1;
415 
416  for (int j = 0; j < kmax - 1; ++j) {
417  double dds = fabs(TDa[j]) + fabs(TDa[j + 1]);
418 
419  if (fabs(TDb[j]) + dds > dds) {
420  kmin = j + 1;
421 
422  break;
423  }
424  }
425 
426  break;
427  }
428 
429  if (j == kmin) {
430  Nconv = iter;
431  vout.paranoiac(m_vl, " converged at iter = %d\n", Nconv);
432 
433  return;
434  }
435  } // end of j loop
436  } // end of iter loop
437 
438  if (Nconv == -1) {
439  vout.crucial(m_vl, "%s: QL method NOT converged, Niter = %d.\n", class_name.c_str(), Niter);
440  abort();
441  }
442 }
443 
444 
445 //====================================================================
446 void Eigensolver_IRLanczos::qrtrf(valarray<double>& TDa,
447  valarray<double>& TDb,
448  int Nk, int Nm, valarray<double>& Qt,
449  double Dsh, int kmin, int kmax)
450 {
451  int k = kmin - 1;
452  double x;
453 
454  double Fden = 1.0 / sqrt((TDa[k] - Dsh) * (TDa[k] - Dsh)
455  + TDb[k] * TDb[k]);
456  double c = (TDa[k] - Dsh) * Fden;
457  double s = -TDb[k] * Fden;
458 
459  double tmpa1 = TDa[k];
460  double tmpa2 = TDa[k + 1];
461  double tmpb = TDb[k];
462 
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;
466  x = -s * TDb[k + 1];
467  TDb[k + 1] = c * TDb[k + 1];
468 
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;
474  }
475 
476  //- Givens transformations
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;
481 
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;
489  if (k != kmax - 2) {
490  x = -s * TDb[k + 1];
491  TDb[k + 1] = c * TDb[k + 1];
492  }
493 
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;
499  }
500  }
501 }
502 
503 
504 //====================================================================
505 //============================================================END=====
Bridge::VerboseLevel m_vl
Definition: eigensolver.h:48
BridgeIO vout
Definition: bridgeIO.cpp:207
void schmidt_orthogonalization(Field &w, std::valarray< Field > &vk, int k)
void Register_string(const string &, const string &)
Definition: parameters.cpp:352
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,...)
Definition: bridgeIO.cpp:38
void Register_int(const string &, const int)
Definition: parameters.cpp:331
Container of Field-type object.
Definition: field.h:37
Class for parameters.
Definition: parameters.h:40
static const std::string class_name
void setUnit_Qt(int Nm, std::valarray< double > &Qt)
int square_non_zero(const double v)
Definition: checker.cpp:41
void set_parameters(const Parameters &params)
dcomplex dotc(const Field &y, const Field &x)
Definition: field.cpp:98
int fetch_string(const string &key, string &val) const
Definition: parameters.cpp:194
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:62
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:193
void crucial(const char *format,...)
Definition: bridgeIO.cpp:26
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)
Definition: checker.cpp:21
void Register_double(const string &, const double)
Definition: parameters.cpp:324
int fetch_double(const string &key, double &val) const
Definition: parameters.cpp:124
string get_string(const string &key) const
Definition: parameters.cpp:85
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
Definition: parameters.cpp:141
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:191
void step(int Nm, int k, std::valarray< double > &TDa, std::valarray< double > &TDb, std::valarray< Field > &vk, Field &f)