102 template <
typename T >
109 const Real& zero_cutoff,
112 int n_min,
int n_max,
118 const Real& delta_cycle,
119 const Real& gamma_factor)
128 const Subset&
s =
A.subset();
131 int N_eig_index = N_eig - 1;
132 int N_eig_minus_one = N_eig_index;
135 psi[
s] = psi_all[N_eig_index];
138 if( N_eig_minus_one > 0 )
151 if( N_eig_minus_one > 0 ) {
165 p[
s] = Apsi - lambda *
psi;
181 QDPIO::cout <<
"Starting Ritz: N_eig=" << N_eig <<
", mu = " <<
mu
182 <<
", g2_0 = " << g2_0 << std::endl;
187 bool minItersDoneP = n_min <= 0;
188 bool maxItersDoneP = n_max <= 0;
192 bool deltaCycleConvP;
198 rsd_a_sq = Rsd_a * Rsd_a;
199 rsd_r_sq = Rsd_r * Rsd_r *
mu *
mu;
205 CGConvP = toBool( g2 < rsd_r_sq ) || toBool(g2 < rsd_a_sq );
211 zeroFoundP = toBool( fabs(
mu) < zero_cutoff );
215 Double delta_cycle_err = delta_cycle;
217 if( Kalk_Sim ==
false ) {
222 convP = minItersDoneP && ( CGConvP || zeroFoundP );
241 convP = minItersDoneP && ( KSConvP || CGConvP || zeroFoundP );
267 if( N_eig_minus_one > 0 ) {
269 GramSchm(Ap, psi_all, N_eig_minus_one,
s);
274 d = innerProductReal(
p, Ap,
s);
287 if( toBool (
a >= s3) )
305 if( toBool( s2 > 0) )
307 s2 = half * (
one+s2);
313 s2 = half * (
one-s2);
333 psi[
s] +=
p * Real(st);
335 Apsi[
s] += Ap * Real(st);
339 Ap[
s] = Apsi -
psi*lambda;
349 minItersDoneP = (
k >= n_min );
350 maxItersDoneP = (
k >= n_max );
353 rsd_a_sq = Rsd_a*Rsd_a;
354 rsd_r_sq = Rsd_r*Rsd_r*
mu*
mu;
360 CGConvP = toBool( g2 < rsd_a_sq ) || toBool( g2 < rsd_r_sq);
366 zeroFoundP = toBool( fabs(
mu) < zero_cutoff );
368 if( Kalk_Sim ==
false ) {
373 convP = minItersDoneP && ( maxItersDoneP || CGConvP || zeroFoundP );
399 deltaCycleConvP =
false;
400 Double g2_g0_ratio = g2 / g2_0;
401 KSConvP = deltaCycleConvP || toBool( g2_g0_ratio <=
Double(gamma_factor)) ;
402 convP = minItersDoneP && ( maxItersDoneP || CGConvP || KSConvP || zeroFoundP );
411 if( N_eig_minus_one > 0 ) {
416 d = sqrt(norm2(
psi,
s));
423 QDPIO::cout <<
"Converged at iter=" <<
k <<
", lambda = " << lambda
424 <<
", rsd | mu | = " << rsd <<
", || g || = "
425 << sqrt(g2) <<
" || x || - 1 = " <<
d << std::endl;
430 QDPIO::cout <<
"KS: gamma = "<< gamma_factor <<
", || g ||^2/|| g_0 ||^2="
432 <<
", delta_cycle_err=" << delta_cycle_err << std::endl;
433 QDPIO::cout <<
"KS: CGConvP = " << CGConvP <<
", KSConvP = " << KSConvP <<
" deltaCycleConvP = " << deltaCycleConvP << std::endl;
440 mu = innerProductReal(
psi, Apsi,
s);
443 QDPIO::cout <<
"Mu-s at convergence: old " << s1 <<
" vs " <<
mu << std::endl;
447 psi_all[N_eig_index][
s] =
psi;
450 Ap[
s] = Apsi -
psi*lambda;
451 final_grad = sqrt(norm2(Ap,
s));
466 if( toBool( ct >
Double(1)) )
469 QDPIO::cout <<
"Restart at iter " <<
k <<
" since beta = " <<
b << std::endl;
482 if(
k % n_renorm == 0 )
486 if( N_eig_minus_one > 0 ) {
490 d = sqrt(norm2(
psi,
s));
502 GramSchm (Apsi, psi_all, N_eig_index,
s);
506 mu = innerProductReal(
psi, Apsi,
s);
510 Ap[
s] = Apsi -
psi * lambda;
516 if( N_eig_minus_one > 0 ) {
531 else if (ProjApsiP && N_eig_minus_one > 0)
534 if( N_eig_minus_one > 0 ) {
546 psi_all[N_eig_index][
s] =
psi;
548 final_grad = sqrt(g2);
549 QDPIO::cerr <<
"too many CG/Ritz iterations: n_count=" <<
n_count
550 <<
", rsd_r =" << sqrt(rsd_r_sq) <<
" rsd_a=" << Rsd_a <<
", ||g||=" << sqrt(g2) <<
", p2=" << p2
551 <<
", lambda" << lambda << std::endl;
561 multi1d<LatticeFermion>& psi_all,
565 const Real& zero_cutoff,
568 int n_min,
int n_max,
574 const Real& delta_cycle,
575 const Real& gamma_factor)
577 Ritz_t(
A, lambda, psi_all, N_eig, Rsd_r, Rsd_a, zero_cutoff,
578 n_renorm, n_min, n_max,
MaxCG,
579 ProjApsiP,
n_count, final_grad,
580 Kalk_Sim, delta_cycle, gamma_factor);
Primary include file for CHROMA library code.
Gramm-Schmidt orthogonolization.
void GramSchm(multi1d< LatticeFermion > &psi, const int Npsi, const multi1d< LatticeFermion > &vec, const int Nvec, const Subset &sub)
Gramm-Schmidt orthogonolization.
void Ritz_t(const LinearOperator< T > &A, Real &lambda, multi1d< T > &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)
Minimizes the Ritz functional with a CG based algorithm.
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
LinOpSysSolverMGProtoClover::T T
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