CHROMA
eig_spec.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Compute low lying eigenvalues of the hermitian H
3  */
4 
5 
6 #include "chromabase.h"
7 
8 #include <sstream>
9 
10 #include "meas/eig/eig_spec.h"
11 #include "meas/eig/ritz.h"
12 #include "meas/eig/sn_jacob.h"
13 
14 namespace Chroma {
15 
16 //! Compute low lying eigenvalues of the hermitian H
17 /*!
18  * \ingroup eig
19  *
20  * Compute low lying eigenvalues of the hermitian H
21  * using the Ritz functional minimization routine,
22  * if desired in the Kalkreuter-Simma algorithm
23  * H The operator (Read)
24  * lambda_H Eigenvalues (Modify)
25  * psi Eigenvectors (Modify)
26  * N_eig No of eigenvalues (Read)
27  * N_min Minimal CG iterations (Read)
28  * N_max Maximal CG iterations (Read)
29  * Kalk_Sim Use Kalkreuter-Simma criterion (Read)
30  * n_renorm Renormalize every n_renorm iter. (Read)
31  * N_KS_max Max number of Kalkreuter-Simma iterations (Read)
32  * RsdR_r (relative) residue (Read)
33  * Cv_fact "Convergence factor" required (Read)
34  * NCG_tot Total number of CG iter (Write)
35  * n_KS total number of Kalkreuter-Simma iterations (Write)
36  * n_valid number of valid eigenvalues (Write)
37  * ProjApsiP flag for projecting A.psi (Read)
38  */
39 
40 
41 void EigSpecRitzCG(const LinearOperator<LatticeFermion>& M, // Herm pos def operator
42  multi1d<Real>& lambda_H, // E-values
43  multi1d<LatticeFermion>& psi, // E-vectors
44  int n_eig, // No of e-values to find
45  int n_renorm, // renorm frequency
46  int n_min, // minimum iters / e_value
47  int MaxCG, // Max no of CG iters
48  const Real& Rsd_r, // relative residuum of each
49  // e-value
50  const Real& Rsd_a, // absolute residuum
51  const Real& zero_cutoff, // if ev slips below this
52  // we consider it zero
53  const bool ProjApsiP, // Project in Ritz?
54 
55  int& n_cg_tot, // Total no of CG iters
56  XMLWriter& xml_out // Diagnostics
57  )
58 {
59  START_CODE();
60 
61  const Subset& sub = M.subset(); // Subset over which M acts
62 
63  push(xml_out, "EigSpecRitzCG");
64 
65 
66  n_cg_tot = 0;
67  int n_count = 0;
68  int n_max = MaxCG+1;
69 
70  int n, i;
71  Real final_grad;
72 
73  multi1d<Real> resid_rel(n_eig);
74  for(n = 1, i = 0; n <= n_eig; n++, i++) {
75 
76  // Initialise lambda[i] = 1
77  lambda_H[i] = Real(1);
78 
79  // Do the Ritz
80  Ritz(M, lambda_H[i], psi, n, Rsd_r, Rsd_a, zero_cutoff, n_renorm, n_min,
81  n_max, MaxCG, ProjApsiP, n_count, final_grad, false, Real(1),
82  Real(1));
83 
84  // Add n_count
85  n_cg_tot += n_count;
86 
87  // Check e-value
88  std::ostringstream s;
89  s << "eval" << n;
90  push(xml_out, s.str());
91  write(xml_out, "n_count", n_count);
92  pop(xml_out);
93 
94  if( toBool( fabs(lambda_H[i]) < zero_cutoff ) )
95  {
96  QDPIO::cout << "Evalue["<< n << "] = " << lambda_H[i] << " is considered zero" << std::endl;
97  }
98  else
99  {
100  LatticeFermion D_e;
101  LatticeFermion lambda_e;
102  M(D_e, psi[i], PLUS);
103  lambda_e[sub] = lambda_H[i]*psi[i];
104  D_e[sub] -= lambda_e;
105  Double r_norm = sqrt(norm2(D_e,sub));
106  resid_rel[i] = Real(r_norm)/lambda_H[i];
107  QDPIO::cout << "Evalue["<<n<<"]: eigen_norm = " << r_norm << " resid_rel = " << resid_rel[i] << std::endl << std::endl;
108  }
109  }
110 
111  // All EV-s done. Dump-em
112  push(xml_out, "eValues");
113  write(xml_out, "lambda", lambda_H);
114  pop(xml_out);
115  write(xml_out, "n_cg_tot", n_cg_tot);
116 
117  pop(xml_out); // EigSpecRitz
118 
119 
120  END_CODE();
121 }
122 
123 
124 void EigSpecRitzKS(const LinearOperator<LatticeFermion>& M, // Herm pos def operator
125  multi1d<Real>& lambda_H, // E-values
126  multi1d<LatticeFermion>& psi, // E-vectors
127  int n_eig, // no of eig wanted
128  int n_dummy, // No of Dummy e-vals to use
129 
130  int n_renorm, // renorm frequency
131  int n_min, // minimum iters / e_value
132  int n_max, // max iters / e_value
133  int n_max_KS, // max KS cycles
134  const Real& gamma_factor, // the KS gamma factor
135 
136  int MaxCG, // Max no of CG iters
137  const Real& Rsd_r, // relative residuum of each
138  // e-value
139  const Real& Rsd_a, // Absolute residuum
140  const Real& zero_cutoff, // if e-value slips below this, we
141  // consider it zero
142  const bool ProjApsiP, // Project in Ritz?
143 
144  int& n_cg_tot, // Total no of CG iters
145  int& n_KS, // Total no of KS cycles
146  int& n_jacob_tot,
147  XMLWriter& xml_out // Diagnostics
148  )
149 {
150  START_CODE();
151 
152  const Subset& s = M.subset(); // Subset over which M acts
153 
154  // Sanity Checks:
155  // Make sure lambda_H is large enough
156  if( lambda_H.size() < (n_eig+n_dummy) ) {
157  QDP_error_exit("lambda_H is too small to hold n_eig + n_dummy values\n");
158  }
159 
160  // Make sure psi is large enough
161  if( psi.size() < (n_eig + n_dummy) ) {
162  QDP_error_exit("psi is too small to hold n_eig + n_dummy values\n");
163  }
164 
165  if( n_eig <=0 ) {
166  QDP_error_exit("n_eig must be > 0. it is %d\n", n_eig);
167  }
168 
169 
170  if( n_eig ==1 ) {
171  // if n_eig is one, KS algorithm is not applicable. We revert to the
172  // Normal CG method
173  EigSpecRitzCG(M, lambda_H, psi, n_eig, n_renorm, n_min, MaxCG, Rsd_r,
174  Rsd_a, zero_cutoff, ProjApsiP, n_cg_tot, xml_out);
175  return;
176  }
177 
178  // Internal lambda's
179  // lambda_H holds initial guesses. Copy these over.
180  int n_working_eig = n_eig+n_dummy;
181  multi1d<Real> lambda_intern(n_working_eig);
182  multi1d<Real> lambda_previous(n_working_eig);
183  multi1d<Real> delta_cycle(n_working_eig);
184  int i;
185  for(i = 0; i < n_working_eig; i++) {
186  lambda_intern[i] = lambda_H[i];
187  lambda_previous[i] = Real(1);
188  delta_cycle[i] = Real(1);
189  }
190 
191  // Off diag elements of hermitian matrix to diagonalise
192  multi1d<Complex> off_diag((n_working_eig)*(n_working_eig-1)/2);
193 
194 
195  multi1d<Real> final_grad(n_working_eig);
196 
197  int n_count = 0;
198  n_cg_tot = 0;
199  n_KS = 0;
200  n_jacob_tot = 0;
201 
202  int ev, j, ij;
203  int n_jacob;
204 
205  push(xml_out, "EigSpecRitzKS");
206 
207 
208  for( int KS_iter = 0; KS_iter < n_max_KS; KS_iter++) {
209  n_KS++;
210 
211  push(xml_out, "KS_iter");
212 
213  // KS Step 1: for each k=0, n-1 in succession compute
214  // approximations to the eigenvectors of M by performing
215  // only a certain number of CG iterations (governed by gamma_factor)
216  for(ev = 1, i=0; ev <= n_working_eig; ev++, i++) {
217 
218  // Do the Ritz for all working eigenvalues
219  Ritz(M, lambda_intern[i], psi, ev, Rsd_r, Rsd_a, zero_cutoff, n_renorm,
220  n_min, n_max, MaxCG, ProjApsiP, n_count, final_grad[i], true,
221  Real(1), gamma_factor);
222 
223  // Count the CG cycles
224  n_cg_tot += n_count;
225  push(xml_out, "ev");
226  write(xml_out, "n_count", n_count);
227  write(xml_out, "final_grad", final_grad[i]);
228  pop(xml_out);
229 
230  }
231 
232  // Construct the hermitian matrix M to diagonalise.
233  // Note only the off diagonal elements are needed
234  // the lambda_intern[i] are the diagonal elements
235 
236  // M_ij = (psi_j, M psi_i )
237  //
238  LatticeFermion tmp;
239  for(i=0, ij=0; i < n_working_eig; i++) {
240  M(tmp, psi[i], PLUS);
241  for(j=0; j < i; j++) {
242  off_diag[ij] = innerProduct(psi[j], tmp);
243  ij++;
244  }
245  }
246 
247 
248 
249  // Now diagonalise it, rotate evecs, and sort
250  //
251  // Jacobi at the moment works with the absolute error
252  n_jacob = SN_Jacob(psi, n_working_eig, lambda_intern, off_diag, Rsd_a, 50, s);
253  n_jacob_tot += n_jacob;
254 
255  write(xml_out, "n_jacob", n_jacob);
256 
257 
258  // Now we should check convergence
259  // Pessimistic but safe criterion || g ||/lambda < omega
260  // but can also be converged if we consider something to have a zero ev
261  // we only check convergence of the wanted n_eig
262  bool convP = true;
263  for(i=0; i < n_eig; i++) {
264  bool zeroReachedP = toBool( fabs(lambda_intern[i]) < zero_cutoff );
265 
266  // This ev is converged if either
267  // i) Pessimistic absolute error is reached
268  // ii) Pessimistic relative error is reached
269  // iii) if the ev is below the zero cutoff
270  //
271  convP &= ( toBool( final_grad[i] < Rsd_a )
272  || toBool( final_grad[i] < Rsd_r*fabs(lambda_intern[i]) )
273  || zeroReachedP );
274  }
275 
276  pop(xml_out); // KS_iter
277 
278  // All the wanted e-values converged
279  if( convP ) {
280 
281  // Do final Jacobi without the dummies.
282  // Make the matrix
283  for(i=0, ij=0; i < n_eig; i++) {
284  M(tmp, psi[i], PLUS);
285  for(j=0; j < i; j++) {
286  off_diag[ij] = innerProduct(psi[j], tmp);
287  ij++;
288  }
289  }
290 
291  // Diagonalise, rotate, sort
292  n_jacob = SN_Jacob(psi, n_eig, lambda_intern, off_diag, Rsd_a, 50, s);
293  write(xml_out, "final_n_jacob", n_jacob);
294  write(xml_out, "n_cg_tot", n_cg_tot);
295  write(xml_out, "n_KS", n_KS);
296 
297 
298  // Copy lambda_intern back into lambda_H and return
299  for(i=0; i < n_eig; i++) {
300  lambda_H[i] = lambda_intern[i];
301  }
302 
303  write(xml_out, "lambda_H", lambda_H);
304  pop(xml_out); // EigSpecRitzKS
305  END_CODE();
306  return;
307  }
308 
309  }
310 
311  write(xml_out, "n_cg_tot", n_cg_tot);
312  write(xml_out, "n_KS", n_KS);
313  pop(xml_out); //EigSpecRitzKS
314 
315  // If we reached here then we have done more than n_max KS
316  QDP_error_exit("n_max_KS reached with no convergence");
317  END_CODE();
318 }
319 
320 
321 void fixMMev2Mev(const LinearOperator<LatticeFermion>& M, // The Op to fix to
322  multi1d<Real>& lambda, // The Evals of M^{dag}M on input
323  // The Evals of M on output
324  multi1d<LatticeFermion>& ev_psi, // The Evecs
325  const int n_eig, // The no of evals/evecs
326  const Real& Rsd_r, // Relative error for validity
327  const Real& Rsd_a, // Absolute error for validity
328  const Real& zero_cutoff, // Zero cutoff
329  multi1d<bool>& valid_eig, // Validity mask (Write)
330  int& n_valid, // No of valids (Write)
331  int& n_jacob // How many Jacobis were done
332  )
333 {
334  START_CODE();
335 
336  const Subset& s = M.subset(); // Subset over which M acts
337 
338  // Sanity checking
339  if( n_eig > lambda.size() ) {
340  QDP_error_exit("n_eig greater than size of lambda array\n");
341  }
342 
343  if( lambda.size() != ev_psi.size() ) {
344  QDP_error_exit("lambda and ev_psi arrays must have same size\n");
345  }
346 
347  // Make the valid eig the right size
348  valid_eig.resize(n_eig);
349 
350  // Temporaries
351  LatticeFermion tmp;
352  Double lambda_H_sq;
353  Double delta_lambda;
354  bool convP;
355  // We are all set -- lets do it
356 
357  if( n_eig == 1 ) { // only one eigenvalue
358  M(tmp, ev_psi[0], PLUS);
359  Double lambda_fix_single = innerProductReal(ev_psi[0], tmp);
360 
361  // No diagonalisation needed -- only 1 eval
362 
363  // We square lambda_fix_single.
364  lambda_H_sq = lambda_fix_single*lambda_fix_single;
365 
366 
367  // Different validity criteria, depending on whether lambda_fix_single^2
368  // is less than our zero cutoff
369  if( toBool( lambda_H_sq < zero_cutoff ) ) {
370 
371 
372  // Yes, the square of our estimate for lambda is less than the cutoff
373  // Check whether the original lambda was also less than the cutoff
374  // if so, our ev is considered a valid zero e-value
375  if( toBool( lambda[0] < zero_cutoff ) ) {
376  valid_eig[0] = true;
377  n_valid = 1;
378  }
379  else {
380  // No, the square of our estimate was not zero. Mismatch
381  // BTW what else can we do here?
382  valid_eig[0] = false;
383  n_valid = 1;
384  }
385  }
386  else {
387  // Check relative error or absolute error is satisfied
388  delta_lambda = fabs( lambda_H_sq - Double( lambda[0] ) );
389 
390  convP = toBool( delta_lambda < Double(Rsd_a) )
391  || toBool( delta_lambda < Double(Rsd_r)*fabs(Double(lambda[0])) );
392 
393  if ( convP ) {
394  // Condition is satisfied
395  valid_eig[0] = true;
396  n_valid = 1;
397  } else {
398  // Condition is not satisfied
399  valid_eig[0] = false;
400  n_valid = 0;
401  }
402  }
403 
404  // Whatever happened lambda_fix_single is our estimate for
405  // lambda[0]
406  lambda[0] = Real(lambda_fix_single);
407  return;
408  } // end if (n_eig == 1)
409 
410 
411  // Interesting case multiple ev-s
412  // Construct new ev-s and off diagonal matrix
413  multi1d<Complex> off_diag(n_eig*(n_eig-1)/2 );
414  multi1d<Real> lambda_fix(n_eig);
415  for(int i = 0; i < n_eig; i++) {
416  valid_eig[i] = false;
417  }
418 
419  for(int i=0, ij=0; i < n_eig; i++) {
420  M(tmp, ev_psi[i], PLUS);
421  lambda_fix[i] = innerProductReal(ev_psi[i], tmp);
422 
423  for(int j=0; j < i; j++) {
424  off_diag[ij] = innerProduct(ev_psi[j], tmp);
425  ij++;
426  }
427  }
428 
429  // Diagonalise and sort according to size
430  n_jacob = SN_Jacob(ev_psi, n_eig, lambda_fix, off_diag, Rsd_a, 50, s);
431 
432  // Now the tricky business of matching up with ev-s of H^2
433  n_valid = 0;
434 
435  for(int i=0; i < n_eig; i++) {
436 
437  lambda_H_sq = Double(lambda_fix[i])*Double(lambda_fix[i]);
438 
439 
440  if( toBool(lambda_H_sq < zero_cutoff ) ) {
441 
442  // Now compare with all the invalid ev-s
443  for(int j = n_valid; ( j < n_eig ) && ( valid_eig[i] == false ); j++ ) {
444 
445  // If we find an original zero lambda, we match it with this
446  // and consider this matched
447  if( toBool( lambda[j] < zero_cutoff ) ) {
448 
449  // lambda_fix[i] matches lambda[j]
450  valid_eig[i] = true;
451 
452  // swap lambda[j] with lambda[valid_eig];
453  // so we don't search it again
454  Real tmp = lambda[j];
455  lambda[j] = lambda[n_valid];
456  lambda[n_valid] = tmp;
457 
458  // Increase the validity count
459  n_valid++;
460  }
461 
462  // Else compare with next j
463  }
464  // At this point either valid_eig[i] is true, or we were unable to match
465  // in which case valid_eig[i] is false from initialisation
466  // in any case time for next_eig
467 
468  }
469  else {
470  // Now compare with all the invalid ev-s
471  for(int j = n_valid; ( j < n_eig ) && ( valid_eig[i] == false ); j++ ) {
472  delta_lambda = fabs( lambda_H_sq - Double(lambda[j]) );
473 
474  convP = toBool( delta_lambda < Double(Rsd_a) )
475  || toBool( delta_lambda < Double(Rsd_r) * fabs( Double(lambda[j]) ) );
476  if( convP ) {
477 
478  // lambda_fix[i] matches lambda[j]
479  valid_eig[i] = true;
480 
481  // swap lambda[j] with lambda[valid_eig];
482  // so we don't search it again
483  Real tmp = lambda[j];
484  lambda[j] = lambda[n_valid];
485  lambda[n_valid] = tmp;
486 
487  // Increase the validity count
488  n_valid++;
489  }
490 
491  // Else compare with next j
492  }
493  // At this point either valid_eig[i] is true, or we were unable to match
494  // in which case valid_eig[i] is false from initialisation
495  // in any case time for next_eig
496  }
497  } // We are done
498 
499  // Now we can overwrite lambda with lambda_fix
500  // A desirable feature would be to move the invalid ev-s to the end
501  // but I cannot be bothered
502  for(int i = 0; i < n_eig; i++ ) {
503  lambda[i] = lambda_fix[i];
504  }
505 
506  END_CODE();
507 
508  // we are done
509  return;
510 }
511 
512 } // end namespace Chroma
Primary include file for CHROMA library code.
virtual const Subset & subset() const =0
Return the subset on which the operator acts.
int SN_Jacob(multi1d< LatticeFermion > &psi, const int N_eig, multi1d< Real > &lambda, multi1d< Complex > &off_diag, Real tolerance, int N_max, const Subset &sub)
Single-node Jacobi rotation.
Definition: sn_jacob.cc:226
void EigSpecRitzCG(const LinearOperator< LatticeFermion > &M, multi1d< Real > &lambda_H, multi1d< LatticeFermion > &psi, int n_eig, int n_renorm, int n_min, int MaxCG, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, const bool ProjApsiP, int &n_cg_tot, XMLWriter &xml_out)
Compute low lying eigenvalues of the hermitian H.
Definition: eig_spec.cc:41
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
void fixMMev2Mev(const LinearOperator< LatticeFermion > &M, multi1d< Real > &lambda, multi1d< LatticeFermion > &ev_psi, const int n_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, multi1d< bool > &valid_eig, int &n_valid, int &n_jacob)
Definition: eig_spec.cc:321
void EigSpecRitzKS(const LinearOperator< LatticeFermion > &M, multi1d< Real > &lambda_H, multi1d< LatticeFermion > &psi, int n_eig, int n_dummy, int n_renorm, int n_min, int n_max, int n_max_KS, const Real &gamma_factor, int MaxCG, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, const bool ProjApsiP, int &n_cg_tot, int &n_KS, int &n_jacob_tot, XMLWriter &xml_out)
Definition: eig_spec.cc:124
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
int n_count
Definition: invbicg.cc:78
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
push(xml_out,"Condensates")
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
Definition: pbg5p_w.cc:32
int i
Definition: pbg5p_w.cc:55
@ PLUS
Definition: chromabase.h:45
void Ritz(const LinearOperator< LatticeFermion > &A, Real &lambda, multi1d< LatticeFermion > &psi_all, int N_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, int n_renorm, int n_min, int n_max, int MaxCG, bool ProjApsiP, int &n_count, Real &final_grad, bool Kalk_Sim, const Real &delta_cycle, const Real &gamma_factor)
Definition: ritz.cc:559
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
START_CODE()
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
Double r_norm
Definition: pade_trln_w.cc:86