CHROMA
syssolver_mdagm_ibicgstab.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_ibicgstab_h__
7 #define __syssolver_mdagm_ibicgstab_h__
8 #include "chroma_config.h"
9 
10 #include "handle.h"
11 #include "syssolver.h"
12 #include "linearop.h"
16 #ifdef CHROMA_DO_ONE_CG_RESTART
18 #endif
19 
20 #include "lmdagm.h"
23 namespace Chroma
24 {
25 
26  //! IBiCGStab system solver namespace
27  namespace MdagMSysSolverIBiCGStabEnv
28  {
29  //! Register the syssolver
30  bool registerAll();
31  }
32 
33 
34  //! Solve a IBiCGStab system. Here, the operator is NOT assumed to be hermitian
35  /*! \ingroup invert
36  */
37  template<typename T>
39  {
40  public:
41  //! Constructor
42  /*!
43  * \param M_ Linear operator ( Read )
44  * \param invParam inverter parameters ( Read )
45  */
47  const SysSolverBiCGStabParams& invParam_) :
48  A(A_), invParam(invParam_)
49  {}
50 
51  //! Destructor is automatic
53 
54  //! Return the subset on which the operator acts
55  const Subset& subset() const {return A->subset();}
56 
57  //! Solver the linear system
58  /*!
59  * \param psi solution ( Modify )
60  * \param chi source ( Read )
61  * \return syssolver results
62  */
64  {
65  START_CODE();
66 
67 
68  StopWatch swatch;
69  SystemSolverResults_t res1,res2,res3; // initialized by a constructo
70  swatch.reset(); swatch.start();
71 
72  T Y;
74  (*A)(Y, psi, PLUS); // Y = M X
75 
77 
78  // Step 2: Solve M X = Y
80 
81  res3.n_count = 0;
82  // Potential safety polishup
83 #ifdef CHROMA_DO_ONE_CG_RESTART
84  // CG Polish - should be very quick
86 #endif
87  res3.n_count += res2.n_count + res1.n_count;
88 
89  { // Find true residuum
90  Y=zero;
91  T re=zero;
92  (*A)(Y, psi, PLUS);
93  (*A)(re,Y, MINUS);
94  re[A->subset()] -= chi;
95  res3.resid = sqrt(norm2(re,A->subset()));
96  }
97 
98  swatch.stop();
99  QDPIO::cout << "IBICGSTAB_SOLVER: " << res3.n_count
100  << " iterations. Rsd = " << res3.resid
101  << " Relative Rsd = " << res3.resid/sqrt(norm2(chi,A->subset())) << std::endl;
102 
103  double time = swatch.getTimeInSeconds();
104  QDPIO::cout << "IBICGSTAB_SOLVER_TIME: "<<time<< " sec" << std::endl;
105 
106 
107  END_CODE();
108 
109  return res3;
110 
111  }
112 
113  //! Solver the linear system
114  /*!
115  * \param psi solution ( Modify )
116  * \param chi source ( Read )
117  * \return syssolver results
118  */
120  AbsChronologicalPredictor4D<T>& predictor) const
121  {
122  START_CODE();
123 
124 
125  StopWatch swatch;
126  SystemSolverResults_t res1,res2,res3; // initialized by a constructo
127  swatch.reset(); swatch.start();
128 
129  T Y = psi;
130 
131  try {
132  // Get a two step solution plan
133  AbsTwoStepChronologicalPredictor4D<T>& two_step_predictor
134  = dynamic_cast<AbsTwoStepChronologicalPredictor4D<T>& >(predictor);
135 
136  // Hooray , we succeeded.
137  // Step 1: Solve M^\dagger Y = chi
138 
139  two_step_predictor.predictY(Y,*A,chi);
141  two_step_predictor.newYVector(Y);
142 
143  // Step 2: Solve M X = Y
145  two_step_predictor.predictX(psi,*MdagM, chi);
147  two_step_predictor.newXVector(psi);
148  }
149  catch(std::bad_cast) {
150 
151  // Boo Hiss, we can't downcast to a two step predictor.
152  // We rely on the fact that we can predict
153  // X ~ (M^\dagger M)^{-1} chi
154  // and then
155  // X ~ M^{-1} M^{-\dagger} chi
156  //
157  // Then MX ~ M^{-\dagger} chi ~ Y
158 
159  T Y ;
161  predictor(psi, (*MdagM), chi);
162  (*A)(Y, psi, PLUS); // Y = M X
164  // Step 2: Solve M X = Y
166 
167  predictor.newVector(psi);
168 
169  }
170 
171  res3.n_count = 0;
172  // Potential safety polishup
173 #ifdef CHROMA_DO_ONE_CG_RESTART
174  // CG Polish - should be very quick
176 #endif
177 
178  res3.n_count += res2.n_count + res1.n_count;
179 
180  { // Find true residuum
181  Y=zero;
182  T re=zero;
183  (*A)(Y, psi, PLUS);
184  (*A)(re,Y, MINUS);
185  re[A->subset()] -= chi;
186  res3.resid = sqrt(norm2(re,A->subset()));
187  }
188 
189  swatch.stop();
190  QDPIO::cout << "IBICGSTAB_SOLVER: " << res3.n_count
191  << " iterations. Rsd = " << res3.resid
192  << " Relative Rsd = " << res3.resid/sqrt(norm2(chi,A->subset())) << std::endl;
193 
194  double time = swatch.getTimeInSeconds();
195  QDPIO::cout << "IBICGSTAB_SOLVER_TIME: "<<time<< " sec" << std::endl;
196 
197 
198  END_CODE();
199 
200  return res3;
201  }
202 
203 
204 
205  private:
206  // Hide default constructor
208 
211  };
212 
213 
214 } // End namespace
215 
216 #endif
217 
Chronological predictor for HMC.
Abstract interface for a Chronological Solution predictor.
virtual void newVector(const T &psi)=0
Abstract interface for a Chronological Solution predictor.
virtual void predictY(T &Y, const T &chi, const Subset &s) const
virtual void newYVector(const T &Y)=0
virtual void predictX(T &X, const T &chi, const Subset &s) const
virtual void newXVector(const T &X)=0
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 IBiCGStab system. Here, the operator is NOT assumed to be hermitian.
const Subset & subset() const
Return the subset on which the operator acts.
SystemSolverResults_t operator()(T &psi, const T &chi) const
Solver the linear system.
~MdagMSysSolverIBiCGStab()
Destructor is automatic.
MdagMSysSolverIBiCGStab(Handle< LinearOperator< T > > A_, const SysSolverBiCGStabParams &invParam_)
Constructor.
SystemSolver disambiguator.
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.
Conjugate-Gradient algorithm for a generic Linear Operator.
Linear Operators.
M^dag*M composition of a linear operator.
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
SystemSolverResults_t InvIBiCGStab(const LinearOperator< LatticeFermionF > &A, const LatticeFermionF &chi, LatticeFermionF &psi, const Real &RsdBiCGStab, int MaxBiCGStab, enum PlusMinus isign)
LinOpSysSolverMGProtoClover::T T
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
START_CODE()
Double zero
Definition: invbicg.cc:106
Params for BiCGStab inverter.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Linear system solvers.
Solve a BICGSTAB system.
Disambiguator for MdagM system solvers.
Handle< LinearOperator< T > > MdagM
Zero initial guess predictor.