CHROMA
syssolver_mdagm_cg_lf_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_cg_lf_clover_h__
7 #define __syssolver_mdagm_cg_lf_clover_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 MdagMSysSolverCGLFCloverEnv
30  {
31  //! Register the syssolver
32  bool registerAll();
33  }
34 
35 
36 
37  //! Solve a system using CG iteration.
38  /*! \ingroup invert
39  */
40 
41  class MdagMSysSolverCGLFClover : 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 SysSolverCGCloverParams& 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 
70  const Q& links = state_->getLinks();
71  for(int mu=0; mu < Nd; mu++) {
72  links_single[mu] = links[mu];
73  }
74 
75 
76  // Links single hold the possibly stouted links
77  // with gaugeBCs applied...
78  // Now I need to create a single prec state...
79  fstate_single = new PeriodicFermState<TF,QF,QF>(links_single);
80 
81  // Make single precision M
83 
84  }
85 
86  //! Destructor is automatic
88 
89  //! Return the subset on which the operator acts
90  const Subset& subset() const {return A->subset();}
91 
92  //! Solver the linear system
93  /*!
94  * \param psi solution ( Modify )
95  * \param chi source ( Read )
96  * \return syssolver results
97  */
99  {
101 
102  START_CODE();
103  StopWatch swatch;
104  swatch.start();
105 
106  TF psi_f = psi;
107  TF chi_f = chi;
108 
109 
110  res = InvCG2( (*M_single),
111  chi_f,
112  psi_f,
113  invParam.RsdCG,
114  invParam.MaxCG);
115 
116 
117  psi = psi_f;
118 
119 
120  {
121  T r;
122  r[A->subset()]=chi;
123  T tmp, tmp1;
124  (*A)(tmp, psi, PLUS);
125  (*A)(tmp1, tmp, MINUS);
126  r[A->subset()] -= tmp1;
127  res.resid = sqrt(norm2(r, A->subset()));
128  }
129 
130  QDPIO::cout << "SINGLE_PREC_CLOVER_CG_SOLVER: " << res.n_count
131  << " iterations. Rsd = " << res.resid
132  << " Relative Rsd = " << res.resid/sqrt(norm2(chi,A->subset())) << std::endl;
133 
134  swatch.stop();
135 
136  double time = swatch.getTimeInSeconds();
137  QDPIO::cout << "SINGLE_PREC_CLOVER_CG_SOLVER_TIME: "<<time<< " sec" << std::endl;
138 
139  END_CODE();
140  return res;
141  }
142 
143  //! Solve the linear system starting with a chrono guess
144  /*!
145  * \param psi solution (Write)
146  * \param chi source (Read)
147  * \param predictor a chronological predictor (Read)
148  * \return syssolver results
149  */
150 
152  AbsChronologicalPredictor4D<T>& predictor) const
153  {
154 
155  START_CODE();
156 
157  // This solver uses InvCG2, so A is just the matrix.
158  // I need to predict with A^\dagger A
159  {
161  predictor(psi, (*MdagM), chi);
162  }
163  // Do solve
164  SystemSolverResults_t res=(*this)(psi,chi);
165 
166  // Store result
167  predictor.newVector(psi);
168  END_CODE();
169  return res;
170  }
171 
172  private:
173  // Hide default constructor
177 
178  // Created and initialized here.
181 
182 
183  };
184 
185 
186 } // End namespace
187 
188 #endif
189 
Abstract interface for a Chronological Solution predictor.
virtual void newVector(const T &psi)=0
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
Solve a system using CG iteration.
Handle< FermState< TF, QF, QF > > fstate_single
const Subset & subset() const
Return the subset on which the operator acts.
Handle< LinearOperator< TF > > M_single
MdagMSysSolverCGLFClover(Handle< LinearOperator< T > > A_, Handle< FermState< T, Q, Q > > state_, const SysSolverCGCloverParams &invParam_)
Constructor.
~MdagMSysSolverCGLFClover()
Destructor is automatic.
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.
Parameters for Clover fermion action.
int mu
Definition: cool.cc:24
Even-odd preconditioned Clover fermion linear operator.
SystemSolverResults_t InvCG2(const LinearOperator< LatticeFermionF > &M, const LatticeFermionF &chi, LatticeFermionF &psi, const Real &RsdCG, int MaxCG)
Conjugate-Gradient (CGNE) algorithm for a generic Linear Operator.
Definition: invcg2.cc:240
Class for counted reference semantics.
Conjugate-Gradient algorithm for a generic Linear Operator.
Linear Operators.
M^dag*M composition of a linear operator.
Nd
Definition: meslate.cc:74
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
@ 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.