42 multi1d<Real>& lambda_H,
43 multi2d<LatticeFermion>&
psi,
51 const Real& zero_cutoff,
59 QDPIO::cout <<
"EigSpecArray" << std::endl;
63 const Subset& sub = M.
subset();
65 push(xml_out,
"EigSpecRitzCG");
74 multi1d<Real> resid_rel(n_eig);
75 for(
n = 1,
i = 0;
n <= n_eig;
n++,
i++)
78 lambda_H[
i] = Real(1);
80 QDPIO::cout <<
"Call ritz n=" <<
n << std::endl;
83 Ritz(M, lambda_H[
i],
psi,
n, Rsd_r, Rsd_a, zero_cutoff, n_renorm, n_min,
84 n_max,
MaxCG, ProjApsiP,
n_count, final_grad,
false, Real(1),
93 push(xml_out,
s.str());
97 if( toBool( fabs(lambda_H[
i]) < zero_cutoff ) ) {
98 QDPIO::cout <<
"Evalue["<<
n <<
"] = " << lambda_H[
i] <<
" is considered zero" << std::endl;
102 multi1d<LatticeFermion> D_e(M.
size());
103 multi1d<LatticeFermion> lambda_e(M.
size());
107 for(
int j=0;
j < M.
size(); ++
j)
109 lambda_e[
j][sub] = lambda_H[
i]*
psi[
i][
j];
110 D_e[
j][sub] -= lambda_e[
j];
113 resid_rel[
i] = Real(
r_norm)/lambda_H[
i];
114 QDPIO::cout <<
"Evalue["<<
n<<
"]: eigen_norm = " <<
r_norm <<
" resid_rel = " << resid_rel[
i] << std::endl << std::endl;
119 push(xml_out,
"eValues");
120 write(xml_out,
"lambda", lambda_H);
122 write(xml_out,
"n_cg_tot", n_cg_tot);
132 multi1d<Real>& lambda_H,
133 multi2d<LatticeFermion>&
psi,
141 const Real& gamma_factor,
147 const Real& zero_cutoff,
149 const bool ProjApsiP,
163 if( lambda_H.size() < (n_eig+n_dummy) ) {
164 QDP_error_exit(
"lambda_H is too small to hold n_eig + n_dummy values\n");
168 if(
psi.size2() < (n_eig + n_dummy) ) {
169 QDP_error_exit(
"psi is too small to hold n_eig + n_dummy values\n");
176 int N5 =
psi.size1();
182 Rsd_a, zero_cutoff, ProjApsiP, n_cg_tot, xml_out);
188 int n_working_eig = n_eig+n_dummy;
189 multi1d<Real> lambda_intern(n_working_eig);
190 multi1d<Real> lambda_previous(n_working_eig);
191 multi1d<Real> delta_cycle(n_working_eig);
193 for(
i = 0;
i < n_working_eig;
i++) {
194 lambda_intern[
i] = lambda_H[
i];
195 lambda_previous[
i] = Real(1);
196 delta_cycle[
i] = Real(1);
200 multi1d<Complex> off_diag((n_working_eig)*(n_working_eig-1)/2);
203 multi1d<Real> final_grad(n_working_eig);
213 push(xml_out,
"EigSpecRitzKS");
216 for(
int KS_iter = 0; KS_iter < n_max_KS; KS_iter++)
220 push(xml_out,
"KS_iter");
225 for(ev = 1,
i=0; ev <= n_working_eig; ev++,
i++) {
228 Ritz(M, lambda_intern[
i],
psi, ev, Rsd_r, Rsd_a, zero_cutoff, n_renorm,
229 n_min, n_max,
MaxCG, ProjApsiP,
n_count, final_grad[
i],
true,
230 Real(1), gamma_factor);
236 write(xml_out,
"final_grad", final_grad[
i]);
247 multi1d<LatticeFermion>
tmp(
N5);
248 for(
i=0, ij=0;
i < n_working_eig;
i++) {
250 for(
j=0;
j <
i;
j++) {
252 for(
int n = 1;
n <
N5;
n++) {
265 n_jacob_tot += n_jacob;
267 write(xml_out,
"n_jacob", n_jacob);
275 for(
i=0;
i < n_eig;
i++) {
276 bool zeroReachedP = toBool( fabs(lambda_intern[
i]) < zero_cutoff );
283 convP &= ( toBool( final_grad[
i] < Rsd_a )
284 || toBool( final_grad[
i] < Rsd_r*fabs(lambda_intern[
i]) )
295 for(
i=0, ij=0;
i < n_eig;
i++) {
297 for(
j=0;
j <
i;
j++) {
299 for(
int n=1;
n <
N5;
n++) {
308 write(xml_out,
"final_n_jacob", n_jacob);
309 write(xml_out,
"n_cg_tot", n_cg_tot);
310 write(xml_out,
"n_KS", n_KS);
314 for(
i=0;
i < n_eig;
i++) {
315 lambda_H[
i] = lambda_intern[
i];
318 write(xml_out,
"lambda_H", lambda_H);
326 write(xml_out,
"n_cg_tot", n_cg_tot);
327 write(xml_out,
"n_KS", n_KS);
337 multi1d<Real>& lambda,
339 multi2d<LatticeFermion>& ev_psi,
343 const Real& zero_cutoff,
344 multi1d<bool>& valid_eig,
355 if( n_eig > lambda.size() ) {
359 if( lambda.size() != ev_psi.size2() ) {
364 valid_eig.resize(n_eig);
366 int N5 = ev_psi.size1();
369 multi1d<LatticeFermion>
tmp(
N5);
377 Double lambda_fix_single = innerProductReal(ev_psi[0][0],
tmp[0],
s);
378 for(
int n=1;
n <
N5;
n++) {
379 lambda_fix_single += innerProductReal(ev_psi[0][
n],
tmp[
n],
s);
384 lambda_H_sq = lambda_fix_single*lambda_fix_single;
389 if( toBool( lambda_H_sq < zero_cutoff ) )
394 if( toBool( lambda[0] < zero_cutoff ) ) {
401 valid_eig[0] =
false;
407 delta_lambda = fabs( lambda_H_sq -
Double( lambda[0] ) );
409 convP = toBool( delta_lambda <
Double(Rsd_a) )
410 || toBool( delta_lambda <
Double(Rsd_r)*fabs(
Double(lambda[0])) );
418 valid_eig[0] =
false;
425 lambda[0] = Real(lambda_fix_single);
432 multi1d<Complex> off_diag(n_eig*(n_eig-1)/2 );
433 multi1d<Real> lambda_fix(n_eig);
434 for(
int i = 0;
i < n_eig;
i++) {
435 valid_eig[
i] =
false;
438 for(
int i=0, ij=0;
i < n_eig;
i++) {
440 lambda_fix[
i] = innerProductReal(ev_psi[
i][0],
tmp[0],
s);
441 for(
int n=1;
n <
N5;
n++) {
442 lambda_fix[
i] += innerProductReal(ev_psi[
i][
n],
tmp[
n],
s);
445 for(
int j=0;
j <
i;
j++) {
447 for(
int n=1;
n <
N5;
n++) {
455 n_jacob =
SN_Jacob_Array(ev_psi, n_eig, lambda_fix, off_diag, Rsd_a, 50,
s);
460 for(
int i=0;
i < n_eig;
i++)
464 if( toBool(lambda_H_sq < zero_cutoff ) )
467 for(
int j = n_valid; (
j < n_eig ) && ( valid_eig[
i] ==
false );
j++ )
471 if( toBool( lambda[
j] < zero_cutoff ) ) {
478 Real
tmp = lambda[
j];
479 lambda[
j] = lambda[n_valid];
480 lambda[n_valid] =
tmp;
496 for(
int j = n_valid; (
j < n_eig ) && ( valid_eig[
i] ==
false );
j++ )
498 delta_lambda = fabs( lambda_H_sq -
Double(lambda[
j]) );
500 convP = toBool( delta_lambda <
Double(Rsd_a) )
501 || toBool( delta_lambda <
Double(Rsd_r) * fabs(
Double(lambda[
j]) ) );
509 Real ftmp = lambda[
j];
510 lambda[
j] = lambda[n_valid];
511 lambda[n_valid] = ftmp;
528 for(
int i = 0;
i < n_eig;
i++ ) {
529 lambda[
i] = lambda_fix[
i];
Primary include file for CHROMA library code.
Linear Operator to arrays.
virtual int size() const =0
Expected length of array index.
virtual const Subset & subset() const =0
Return the subset on which the operator acts.
int SN_Jacob_Array(multi2d< 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.
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.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams ¶m)
Writer parameters.
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
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)
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)
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)
push(xml_out,"Condensates")
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
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)
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double