52 LatticeFermion
r =
b - zeta*
x;
58 Real delta = sqrt(norm2(
r));
61 Complex phihat = Real(1)/delta;
64 Complex tauhat = delta/rho;
66 LatticeFermion w_old,
p, v_old;
76 Complex
phi = Real(0);
78 Complex lambda = Real(0);
82 Complex r1_old = Real(1);
83 Complex gamma = Real(1);
88 LatticeFermion v, vtilde;
89 Real ftmp = Real(1)/delta;
97 Real taumod = sqrt(real(conj(tauhat)*tauhat));
98 for(
int iter = 1; iter <= MaxSUMR && !convP; iter++) {
103 Real inner_tol =
epsilon/taumod;
115 sigma = sqrt( Real(1) - real( conj(gamma)*gamma ) );
118 Complex
alpha = - gamma * delta;
124 Complex r1hat =
alpha*phihat+conj(
c)*zeta/rho;
133 Real abs_rhat_sq = real(conj(r1hat)*r1hat);
134 Real abs_sigma_sq = sigma*sigma;
135 Real tmp_length = sqrt(abs_rhat_sq + abs_sigma_sq);
137 c = conj(r1hat)/tmp_length;
138 s = - sigma / tmp_length;
141 Complex r1 = -
c*r1hat + cmplx(
s*sigma,0);
144 Complex tau = -
c*tauhat;
150 Complex
eta = tau / r1;
162 LatticeFermion w_minus_v;
163 w_minus_v = w_old - v_old;
166 p =
p - lambda*w_minus_v;
173 x =
x -
eta*w_minus_v;
184 if( toBool( sqrt(real(conj(sigma)*sigma)) >
epsilon) ) {
190 phi = -
c*phihat +
s*conj(gamma)/delta;
196 phihat =
s* phihat + conj(
c) * conj(gamma)/delta;
202 v =
u + gamma*vtilde;
203 Real
ftmp2 = Real(1)/sigma;
207 Complex gconj = conj(gamma);
208 Complex csigma=cmplx(sigma,0);
209 vtilde = csigma*vtilde + gconj*v;
213 Real ftmp3 = Real(1)/sqrt(norm2(vtilde));
220 taumod = sqrt(real(conj(tauhat)*tauhat));
221 if ( toBool(taumod <
epsilon) ) convP =
true;
223 QDPIO::cout <<
"Iter " << iter <<
": | tauhat |=" << taumod << std::endl;
237 if(
n_count >= MaxSUMR && ! convP ) {
238 QDPIO::cout <<
"Solver Nonconvergence Warning " << std::endl;
245 const LatticeFermion&
b,
Primary include file for CHROMA library code.
Conjugate-Gradient algorithm for a generic Linear Operator.
int epsilon(int i, int j, int k)
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
static const LatticeInteger & alpha(const int dim)
Asqtad Staggered-Dirac operator.
static multi1d< LatticeColorMatrix > u
void InvRelSUMR(const LinearOperator< LatticeFermion > &U, const LatticeFermion &b, LatticeFermion &x, const Complex &zeta, const Real &rho, const Real &epsilon, int MaxSUMR, int &n_count)
LinOpSysSolverMGProtoClover::T T
void InvRelSUMR_a(const LinearOperator< T > &U, const T &b, T &x, const Complex &zeta, const Real &rho, const Real &epsilon, int MaxSUMR, int &n_count)
multi1d< LatticeFermion > s(Ncb)
multi1d< LatticeColorMatrix > U