CHROMA
syssolver_linop_rel_cg_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_linop_rel_cg_multiprec_h__
7 #define __syssolver_linop_rel_cg_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"
21 
22 #include <string>
23 
24 namespace Chroma
25 {
26 
27  //! Richardson system solver namespace
28  namespace LinOpSysSolverReliableCGCloverEnv
29  {
30  //! Register the syssolver
31  bool registerAll();
32  }
33 
34 
35 
36  //! Solve a system using Richardson iteration.
37  /*! \ingroup invert
38  *** WARNING THIS SOLVER WORKS FOR CLOVER FERMIONS ONLY ***
39  */
40 
41  class LinOpSysSolverReliableCGClover : public LinOpSystemSolver<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 SysSolverReliableCGCloverParams& 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  const Q& links = state_->getLinks();
72  for(int mu=0; mu < Nd; mu++) {
73  links_single[mu] = links[mu];
74  links_double[mu] = links[mu];
75  }
76 
77 
78  // Links single hold the possibly stouted links
79  // with gaugeBCs applied...
80  // Now I need to create a single prec state...
81  fstate_single = new PeriodicFermState<TF,QF,QF>(links_single);
82  fstate_double = new PeriodicFermState<TD,QD,QD>(links_double);
83 
84  // Make single precision M
87 
88 
89 
90  }
91 
92  //! Destructor is automatic
94 
95  //! Return the subset on which the operator acts
96  const Subset& subset() const {return A->subset();}
97 
98  //! Solver the linear system
99  /*!
100  * \param psi solution ( Modify )
101  * \param chi source ( Read )
102  * \return syssolver results
103  */
105  {
107 
108  START_CODE();
109  StopWatch swatch;
110  swatch.start();
111 
112  // T MdagChi;
113 
114  // This is a CGNE. So create new RHS
115  // (*A)(MdagChi, chi, MINUS);
116  // Handle< LinearOperator<T> > MM(new MdagMLinOp<T>(A));
117 
118  TD psi_d = psi;
119 
120  // Do this in floating...
121  T M_dag_chi;
122  (*A)(M_dag_chi, chi, MINUS);
123 
124  // Then convert to double
125  TD chi_d = M_dag_chi;
126 
127  res=InvCGReliable(*M_double,
128  *M_single,
129  chi_d,
130  psi_d,
132  invParam.Delta,
133  invParam.MaxIter);
134 
135  psi = psi_d;
136 
137  swatch.stop();
138 
139  {
140  T r;
141  r[A->subset()]=chi;
142  T tmp;
143  (*A)(tmp, psi, PLUS);
144  r[A->subset()] -= tmp;
145  res.resid = sqrt(norm2(r, A->subset()));
146  }
147  QDPIO::cout << "RELIABLE_CGNE_SOLVER: " << res.n_count << " iterations. Rsd = " << res.resid << " Relative Rsd = " << res.resid/sqrt(norm2(chi,A->subset())) << std::endl;
148 
149 
150  END_CODE();
151  return res;
152  }
153 
154 
155  private:
156  // Hide default constructor
160 
161  // Created and initialized here.
166 
167 
168  };
169 
170 
171 } // End namespace
172 
173 #endif
174 
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
Solve a system using Richardson iteration.
LinOpSysSolverReliableCGClover(Handle< LinearOperator< T > > A_, Handle< FermState< T, Q, Q > > state_, const SysSolverReliableCGCloverParams &invParam_)
Constructor.
SystemSolverResults_t operator()(T &psi, const T &chi) const
Solver the linear system.
Handle< FermState< TD, QD, QD > > fstate_double
Handle< FermState< TF, QF, QF > > fstate_single
const SysSolverReliableCGCloverParams invParam
const Subset & subset() const
Return the subset on which the operator acts.
SystemSolver disambiguator.
Linear Operator.
Definition: linearop.h:27
Periodic version of FermState.
Parameters for Clover fermion action.
int mu
Definition: cool.cc:24
Even-odd preconditioned Clover fermion linear operator.
SystemSolverResults_t InvCGReliable(const LinearOperator< LatticeFermionF > &A, const LatticeFermionF &chi, LatticeFermionF &psi, const Real &RsdCG, const Real &Delta, int MaxCG)
Bi-CG stabilized.
Definition: reliable_cg.cc:193
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
@ 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.
BiCGStab Solver with reliable updates.
Support class for fermion actions and linear operators.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Linear system solvers.
Factory for solving M*psi=chi where M is not hermitian or pos. def.