CHROMA
syssolver_mdagm_richardson_multiprec_clover.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! \file
3  * \brief Solve a MdagM*psi=chi linear system by BiCGStab
4  */
5 
6 #ifndef __syssolver_mdagm_richardson_multiprec_h__
7 #define __syssolver_mdagm_richardson_multiprec_h__
8 #include "chroma_config.h"
9 
10 #include "handle.h"
11 #include "state.h"
12 #include "syssolver.h"
13 #include "linearop.h"
14 #include "lmdagm.h"
22 
23 #include <string>
24 
25 namespace Chroma
26 {
27 
28  //! Richardson system solver namespace
29  namespace MdagMSysSolverRichardsonCloverEnv
30  {
31  //! Register the syssolver
32  bool registerAll();
33  }
34 
35 
36 
37  //! Solve a system using Richardson iteration.
38  /*! \ingroup invert
39  */
40 
41  class MdagMSysSolverRichardsonClover : public MdagMSystemSolver<LatticeFermion>
42  {
43  public:
44  typedef LatticeFermion T;
45  typedef LatticeColorMatrix U;
46  typedef multi1d<LatticeColorMatrix> Q;
47 
48  typedef LatticeFermionF TF;
49  typedef LatticeColorMatrixF UF;
50  typedef multi1d<LatticeColorMatrixF> QF;
51 
52  typedef LatticeFermionD TD;
53  typedef LatticeColorMatrixD UD;
54  typedef multi1d<LatticeColorMatrixD> QD;
55 
56  //! Constructor
57  /*!
58  * \param M_ Linear operator ( Read )
59  * \param invParam inverter parameters ( Read )
60  */
62  Handle< FermState<T,Q,Q> > state_,
63  const SysSolverRichardsonCloverParams& invParam_) :
64  A(A_), invParam(invParam_)
65  {
66 
67  // Get the Links out of the state and convertto single.
68  QF links_single; links_single.resize(Nd);
69  QD links_double; links_double.resize(Nd);
70 
71 
72  const Q& links = state_->getLinks();
73  for(int mu=0; mu < Nd; mu++) {
74  links_single[mu] = links[mu];
75  links_double[mu] = links[mu];
76  }
77 
78 
79  // Links single hold the possibly stouted links
80  // with gaugeBCs applied...
81  // Now I need to create a single prec state...
82  fstate_single = new PeriodicFermState<TF,QF,QF>(links_single);
83  fstate_double = new PeriodicFermState<TD,QD,QD>(links_double);
84 
85  // Make single precision M
88 
89  std::istringstream is( invParam_.innerSolverParams.xml );
90  XMLReader paramtop(is);
91 
92  DInv = TheMdagMFermFSystemSolverFactory::Instance().createObject( invParam_.innerSolverParams.id, paramtop,
93  invParam_.innerSolverParams.path,
95  M_single);
96 
97 
98  }
99 
100  //! Destructor is automatic
102 
103  //! Return the subset on which the operator acts
104  const Subset& subset() const {return A->subset();}
105 
106  //! Solver the linear system
107  /*!
108  * \param psi solution ( Modify )
109  * \param chi source ( Read )
110  * \return syssolver results
111  */
113  {
115 
116  START_CODE();
117  StopWatch swatch;
118  swatch.start();
120 
121  TD psi_d = psi;
122  TD chi_d = chi;
123 
125  *MM_double,
126  chi_d,
127  psi_d,
130  res);
131 
132  psi = psi_d;
133  swatch.stop();
134  double time = swatch.getTimeInSeconds();
135  {
136  T r;
137  r[A->subset()]=chi;
138  T tmp, tmp1;
139  (*A)(tmp, psi, PLUS);
140  (*A)(tmp1, tmp, MINUS);
141  r[A->subset()] -= tmp1;
142  res.resid = sqrt(norm2(r, A->subset())/norm2(chi, A->subset()));
143  }
144 
145  QDPIO::cout << "MULTIPREC_RICHARDSON_SOLVER: " << res.n_count << " iterations. Rsd = " << res.resid << " Relative Rsd = " << res.resid/sqrt(norm2(chi,A->subset())) << std::endl;
146  QDPIO::cout << "MULTIPREC_RICHARDSON_SOLVER_TIME: "<<time<< " sec" << std::endl;
147 
148  END_CODE();
149  return res;
150  }
151 
152  //! Solve the linear system starting with a chrono guess
153  /*!
154  * \param psi solution (Write)
155  * \param chi source (Read)
156  * \param predictor a chronological predictor (Read)
157  * \return syssolver results
158  */
159 
161  AbsChronologicalPredictor4D<T>& predictor) const
162  {
163 
164  START_CODE();
165 
166  // This solver uses InvCG2, so A is just the matrix.
167  // I need to predict with A^\dagger A
168  {
170  predictor(psi, (*MdagM), chi);
171  }
172  // Do solve
173  SystemSolverResults_t res=(*this)(psi,chi);
174 
175  // Store result
176  predictor.newVector(psi);
177  END_CODE();
178  return res;
179  }
180 
181  private:
182  // Hide default constructor
186 
187  // Created and initialized here.
193 
194 
195  };
196 
197 
198 } // End namespace
199 
200 #endif
201 
Abstract interface for a Chronological Solution predictor.
virtual void newVector(const T &psi)=0
Even-odd preconditioned Clover-Dirac operator.
Even-odd preconditioned Clover-Dirac operator.
Support class for fermion actions and linear operators.
Definition: state.h:94
Class for counted reference semantics.
Definition: handle.h:33
Linear Operator.
Definition: linearop.h:27
M^dag.M linear operator.
Definition: lmdagm.h:22
const Subset & subset() const
Return the subset on which the operator acts.
MdagMSysSolverRichardsonClover(Handle< LinearOperator< T > > A_, Handle< FermState< T, Q, Q > > state_, const SysSolverRichardsonCloverParams &invParam_)
Constructor.
SystemSolverResults_t operator()(T &psi, const T &chi, AbsChronologicalPredictor4D< T > &predictor) const
Solve the linear system starting with a chrono guess.
SystemSolverResults_t operator()(T &psi, const T &chi) const
Solver the linear system.
SystemSolver disambiguator.
Periodic version of FermState.
static T & Instance()
Definition: singleton.h:432
Parameters for Clover fermion action.
int mu
Definition: cool.cc:24
Even-odd preconditioned Clover fermion linear operator.
Class for counted reference semantics.
Linear Operators.
M^dag*M composition of a linear operator.
Nd
Definition: meslate.cc:74
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
void InvMultiPrecRichardson(const SystemSolver< LatticeFermionF > &Dinv, const LinearOperator< LatticeFermionD > &D, const LatticeFermionD &b, LatticeFermionD &x, int MaxIter, Real RsdTarget, SystemSolverResults_t &res)
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
START_CODE()
Periodic ferm state and a creator.
Support class for fermion actions and linear operators.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Linear system solvers.
Disambiguator for MdagM system solvers.
Handle< LinearOperator< T > > MdagM
Factory for producing system solvers for MdagM*psi = chi.