33 const multi1d<Complex>& zeta,
34 const multi1d<Real>& rho,
42 int numroot =
x.size();
43 if( zeta.size() != numroot ) {
44 QDPIO::cerr <<
"zeta size:"<< zeta.size()
45 <<
" is different from x.size:"<<numroot << std::endl;
50 if( rho.size() != numroot) {
51 QDPIO::cerr <<
"rho size:"<< rho.size()
52 <<
" is different from x.size:"<<numroot << std::endl;
56 if(
epsilon.size() != numroot) {
57 QDPIO::cerr <<
"epsilon size:"<<
epsilon.size()
58 <<
" is different from x.size:"<<numroot << std::endl;
64 Real epsilon_sigma =
epsilon[0];
65 for(
int shift = 0; shift < numroot; shift++) {
66 if ( toBool(
epsilon[shift] < epsilon_sigma ) ) {
79 for(
int shift = 0; shift < numroot; shift++) {
87 Real delta = sqrt(norm2(
r));
89 multi1d<Complex> phihat(numroot);
90 multi1d<Complex> tauhat(numroot);
94 for(
int shift = 0; shift < numroot; shift++) {
96 phihat[shift] = Real(1)/delta;
99 tauhat[shift] = delta/rho[shift];
103 LatticeFermion v_old =
zero;
107 multi1d<LatticeFermion> w_old(numroot);
108 multi1d<LatticeFermion>
p(numroot);
114 for(
int shift = 0; shift < numroot; shift++) {
122 Complex gamma = Real(1);
123 Real sigma = Real(1);
126 multi1d<Complex>
phi(numroot);
127 multi1d<Real>
s(numroot);
128 multi1d<Complex> lambda(numroot);
129 multi1d<Complex>
r0(numroot);
130 multi1d<Complex> r1_old(numroot);
131 multi1d<Complex>
c(numroot);
135 for(
int shift = 0; shift < numroot; shift++) {
136 phi[shift] = Real(0);
138 lambda[shift] = Real(0);
140 r1_old[shift] = Real(1);
150 LatticeFermion vtilde;
152 Real ftmp = Real(1)/delta;
159 multi1d<bool> convP(numroot);
161 bool allConvP =
false;
164 for(
int iter = 1; iter <= MaxSUMR && !allConvP; iter++) {
179 sigma = sqrt( Real(1) - real( conj(gamma)*gamma ) );
183 alpha = -gamma*delta;
185 multi1d<Complex> r1hat(numroot);
190 Real abs_sigma_sq = sigma*sigma;
192 multi1d<Complex> r1(numroot);
195 for(
int shift = 0; shift < numroot; shift++) {
198 if( !convP[shift] ) {
201 r0[shift] =
alpha *
phi[shift] +
s[shift]*zeta[shift]/rho[shift];
204 r1hat =
alpha * phihat[shift]+conj(
c[shift])*zeta[shift]/rho[shift];
212 abs_rhat_sq = real(conj(r1hat[shift])*r1hat[shift]);
213 tmp_length = sqrt(abs_rhat_sq + abs_sigma_sq);
215 c[shift] = conj(r1hat[shift])/tmp_length;
216 s[shift] = - sigma / tmp_length;
219 r1[shift] = -
c[shift]*r1hat[shift] + cmplx(
s[shift]*sigma,0);
222 Complex tau = -
c[shift]*tauhat[shift];
225 tauhat[shift] =
s[shift] * tauhat[shift];
228 Complex
eta = tau / r1[shift];
231 Complex
kappa =
r0[shift]/r1_old[shift];
234 r1_old[shift]=r1[shift];
240 LatticeFermion w_minus_v;
241 w_minus_v = w_old[shift] - v_old;
243 LatticeFermion w =
alpha*
p[shift] -
kappa*w_minus_v;
244 p[shift] =
p[shift] - lambda[shift]*w_minus_v;
251 x[shift] =
x[shift] -
eta*w_minus_v;
265 if( toBool( sqrt(real(conj(sigma)*sigma)) > epsilon_sigma ) ) {
273 for(
int shift=0; shift < numroot; shift++) {
276 if( !convP[shift] ) {
279 phi[shift] = -
c[shift]*phihat[shift]
280 +
s[shift]*conj(gamma)/delta;
283 lambda[shift] =
phi[shift]/r1[shift];
286 phihat[shift] =
s[shift]* phihat[shift]
287 + conj(
c[shift]) * conj(gamma)/delta;
296 v =
u + gamma*vtilde;
297 Real
ftmp2 = Real(1)/sigma;
301 Complex gconj = conj(gamma);
302 Complex csigma=cmplx(sigma,0);
303 vtilde = csigma*vtilde + gconj*v;
307 Real ftmp3 = Real(1)/sqrt(norm2(vtilde));
319 for(
int shift=0; shift < numroot; shift++) {
322 if( ! convP[shift] ) {
325 Real taumod = sqrt(real(conj(tauhat[shift])*tauhat[shift]));
327 QDPIO::cout <<
"Iter " << iter <<
": Shift: " << shift
328 <<
" | tauhat |=" << taumod << std::endl;
331 if ( toBool(taumod <
epsilon[shift] ) ) convP[shift] =
true;
334 allConvP &= convP[shift];
351 if(
n_count == MaxSUMR && ! allConvP ) {
352 QDPIO::cout <<
"Solver Nonconvergence Warning " << std::endl;
360 const LatticeFermion&
b,
361 multi1d<LatticeFermion>&
x,
362 const multi1d<Complex>& zeta,
363 const multi1d<Real>& rho,
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 MInvSUMR(const LinearOperator< LatticeFermion > &U, const LatticeFermion &b, multi1d< LatticeFermion > &x, const multi1d< Complex > &zeta, const multi1d< Real > &rho, const multi1d< Real > &epsilon, int MaxSUMR, int &n_count)
LinOpSysSolverMGProtoClover::T T
multi1d< LatticeFermion > s(Ncb)
void MInvSUMR_a(const LinearOperator< T > &U, const T &b, multi1d< T > &x, const multi1d< Complex > &zeta, const multi1d< Real > &rho, const multi1d< Real > &epsilon, int MaxSUMR, int &n_count)
multi1d< LatticeColorMatrix > U