5 #error "THIS ROUTINE IS NOT FULLY CONVERTED"
77 const multi1d<Real>& shifts,
78 const multi1d<Real>&
RsdMR,
84 const Subset& sub =
A.subset();
86 if (shifts.size() !=
RsdMR.size())
88 QDPIO::cerr <<
"MInvMR: number of shifts and residuals must match" << std::endl;
92 int n_shift = shifts.size();
96 QDPIO::cerr <<
"MInvMR: You must supply at least 1 mass: mass.size() = "
97 << n_shift << std::endl;
107 multi1d<Real> ass(Nmass);
108 multi1d<Complex> rhos(Nmass);
111 multi1d<Boolean> convsP(Nmass);
112 multi1d<Real> rsdmr_sq(Nmass);
124 T Ar; moveToFastMemoryHint(Ar);
125 T chi_internal; moveToFastMemoryHint(chi_internal);
126 chi_internal[sub] =
chi;
127 moveToFastMemoryHint(
psi,
true);
129 FlopCounter flopcount;
136 rsdmr_sq = localNorm2(
RsdMR);
141 for(
int i= 0;
i < n_shift; ++
i) {
148 r[sub] = chi_internal;
155 for(
k = 1;
k <= MaxMR && ! convP ; ++
k)
158 cp = norm2(
r, sub); flopcount.addSiteFlops(4*Nc*Ns,sub);
162 A(Ar,
r,
PLUS); flopcount.addFlops(
A.nFlops());
163 Ar[sub] +=
r *
mass[
isz]; flopcount.addSiteFlops(4*Nc*Ns,sub);
166 DComplex
c =
innerProduct(Ar,
r, sub); flopcount.addSiteFlops(4*Nc*Ns,sub);
169 d = norm2(Ar, sub); flopcount.addSiteFlops(4*Nc*Ns,sub);
176 for(
s = 0;
s < Nmass; ++
s)
178 if (convsP[
s])
continue;
182 atmp += adj(
a) * asq;
183 asq = localNorm2(atmp);
188 as = rhos[
s] * atmp2;
189 atmp2 = rhos[
s] * atmp;
192 ass[
s] = localNorm2(as);
209 for(
s = 0;
s < Nmass; ++
s)
211 if (convsP[
s])
continue;
220 Real cs = Real(ass[
s]) *
cp;
224 PRINTF(
"MInvMR: k = %d s = %d r = %g d = %g\n",
k,
s,sqrt(cs),sqrt(
d));
225 FLUSH_WRITE_NAMELIST(stdout);
236 for(
int s = 0;
s < Nmass; ++
s)
242 cp = norm2(Ar, sub); flopcount.addSiteFlops(4*Nc*Ns,sub);
243 QDPIO::cout <<
"MInvMR (conv): s = " <<
s <<
" r = " << sqrt(
cp) <<
" m = " <<
mass[
s] << std::endl;
259 const LatticeFermion&
chi,
260 multi1d<LatticeFermion>&
psi,
261 const multi1d<Real>& shifts,
262 const multi1d<Real>&
RsdMR,
273 multi1d<LatticeColorMatrix>,
274 multi1d<LatticeColorMatrix> >& M,
275 const LatticeFermion&
chi,
276 multi1d<LatticeFermion>&
psi,
277 const multi1d<Real>& shifts,
278 const multi1d<Real>&
RsdMR,
Differentiable Linear Operator.
void GramSchm(multi1d< LatticeFermion > &psi, const int Npsi, const multi1d< LatticeFermion > &vec, const int Nvec, const Subset &sub)
Gramm-Schmidt orthogonolization.
void MInvMR_a(const LinearOperator< T > &A, const T &chi, multi1d< T > &psi, const multi1d< Real > &shifts, const multi1d< Real > &RsdMR, int MaxMR, int &n_count)
Multishift Conjugate-Gradient (CG1) algorithm for a Linear Operator.
void MInvMR(const LinearOperator< LatticeFermion > &M, const LatticeFermion &chi, multi1d< LatticeFermion > &psi, const multi1d< Real > &shifts, const multi1d< Real > &RsdMR, int MaxMR, int &n_count)
Multishift Conjugate-Gradient (CG1) algorithm for a Linear Operator.
Multishift Conjugate-Gradient algorithm for a Linear Operator.
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
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)
LinOpSysSolverMGProtoClover::T T
multi1d< LatticeFermion > chi(Ncb)
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
PRINTF("TrLnSolve: n_count = %d, RsdMR = %16.8e, MaxRsd = %16.8e\n", n_count, RsdMR, MaxRsd)