CHROMA
syssolver_mdagm_OPTeigcg.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /*! \file
3  * \brief Solve a M^dag*M*psi=chi linear system by EigCG
4  */
5 
6 #ifndef __syssolver_mdagm_OPTeigcg_h__
7 #define __syssolver_mdagm_OPTeigcg_h__
8 #include "chroma_config.h"
9 
10 #include "handle.h"
11 #include "syssolver.h"
12 #include "linearop.h"
13 #include "lmdagm.h"
14 #include "named_obj.h"
16 
20 
21 #include "util/info/unique_id.h"
22 
23 namespace Chroma
24 {
25 
26  //! Eigenstd::vector accelerated CG system solver namespace
27  namespace MdagMSysSolverOptEigCGEnv
28  {
29  //! Register the syssolver
30  bool registerAll();
31  }
32 
33 
34  //! Solve a M*psi=chi linear system by CG2 with eigenvectors
35  /*! \ingroup invert
36  */
37  template<typename T>
39  {
40  public:
41 
42  //! Write out an OptEigInfo Type
44  StopWatch swatch;
45  swatch.reset();
46  swatch.start();
47 
48  // A shorthand for the object
49  const LinAlg::OptEigInfo& obj =
51 
52  // File XML
53  XMLBufferWriter file_xml;
54  push(file_xml, "OptEigInfo");
55  write(file_xml, "id", uniqueId());
56  write(file_xml, "N", obj.N);
57  write(file_xml, "ncurEvals", obj.ncurEvals);
58  write(file_xml, "restartTol", obj.restartTol);
59  write(file_xml, "lde", obj.lde);
60  write(file_xml, "ldh", obj.evals.size());
61  pop(file_xml);
62 
63  // Open file
64  QDPFileWriter to(file_xml,
67  QDPIO_SERIAL,QDPIO_OPEN);
68 
69  for(int v(0);v<obj.ncurEvals;v++){
70  LatticeFermion lf ;
71  obj.CvToLatFerm(lf,subset(),v);
72  XMLBufferWriter record_xml;
73  push(record_xml, "EigenVector");
74  write(record_xml,"no",v);
75  pop(record_xml);
76  write(to, record_xml, lf);
77  }
78 
79  {
80  XMLBufferWriter record_xml;
81  push(record_xml, "EigenValues");
82  pop(record_xml);
83  write(to, record_xml, obj.evals);
84  }
85  {
86  XMLBufferWriter record_xml;
87  push(record_xml, "H");
88  pop(record_xml);
89  write(to, record_xml, obj.H);
90  }
91  {
92  XMLBufferWriter record_xml;
93  push(record_xml, "HU");
94  pop(record_xml);
95  write(to, record_xml, obj.HU);
96  }
97 
98  // Close
99  close(to);
100  swatch.stop();
101  QDPIO::cout<<" QIOWriteOptEvecs: Time to write evecs= "
102  << swatch.getTimeInSeconds() <<" secs "<<std::endl ;
103  }
104 
105  //-----------------------------------------------------------------------
106  //! Read a OptEigInfo Type
108  {
109 
110  StopWatch swatch;
111  swatch.reset();
112  swatch.start();
113 
114  // File XML
115  XMLReader file_xml;
116 
117  // Open file
118  QDPFileReader to(file_xml,invParam.file.file_name,QDPIO_SERIAL);
119 
120  // A shorthand for the object
123 
124  XMLReader record_xml;
125  //TheNamedObjMap::Instance().get(buffer_id).setRecordXML(record_xml);
126 
127  int ldh,lde,N ;
128  float restartTol ;
129  read(file_xml, "/OptEigInfo/ldh", ldh);
130  read(file_xml, "/OptEigInfo/lde", lde);
131  read(file_xml, "/OptEigInfo/N", N);
132  read(file_xml, "/OptEigInfo/ncurEvals", obj.ncurEvals);
133  read(file_xml, "/OptEigInfo/restartTol", restartTol);
134 
135  QDPIO::cout <<__func__ << " : Reading object with following properties"<<std::endl ;
136  QDPIO::cout <<__func__ << " : lde= " << lde << std::endl;
137  QDPIO::cout <<__func__ << " : ldh= " << ldh << std::endl;
138  QDPIO::cout <<__func__ << " : N = " << N << std::endl;
139 
140  QDPIO::cout << __func__ << " : ncurEvals = " << obj.ncurEvals<< std::endl;
141  QDPIO::cout << __func__ << " : restartTol = " << restartTol<< std::endl;
142 
143  if(obj.evals.size()< obj.ncurEvals){
144  QDPIO::cerr<<__func__<< " : ldh of the current object is not large enough to hold the vectors" ;
145  QDP_abort(1);
146  }
147 
148  for(int v(0);v<obj.ncurEvals;v++){
149  LatticeFermion lf ;
150  XMLReader record_xml;
151  read(to, record_xml, lf);
152  obj.CvToEigCGvec(lf, subset(), v) ;
153  }
154  //this is just fine
155  {
156  XMLReader record_xml;
157  multi1d<Real> evals(ldh) ;
158  multi1d<Complex> H(ldh*ldh) ;
159  multi1d<Complex> HU(ldh*ldh) ;
160  read(to, record_xml, evals);
161  read(to, record_xml, H);
162  read(to, record_xml, HU);
163  if(ldh<=obj.evals.size()){
164  for(int i(0);i<ldh;i++)
165  obj.evals[i] = evals[i] ;
166  for(int i(0);i<H.size();i++){
167  obj.H[i] = H[i] ;
168  obj.HU[i] = HU[i] ;
169  }
170  }
171  else{
172  QDPIO::cerr<<__func__<< " : ldh of the current object is not large enough to hold the Cholesky factors" ;
173  QDP_abort(1);
174  }
175  }
176 
177  // Close
178  close(to);
179 
180  swatch.stop();
181  QDPIO::cout<<" QIOReadOptEvecs: Time to read evecs= "
182  << swatch.getTimeInSeconds() <<" secs "<<std::endl ;
183  }
184 
185  //! Constructor
186  /*!
187  * \param M_ Linear operator ( Read )
188  * \param invParam_ inverter parameters ( Read )
189  */
191  const SysSolverOptEigCGParams& invParam_) :
192  MdagM(new MdagMLinOp<T>(A_)), A(A_), invParam(invParam_)
193  {
194 #ifndef QDP_IS_QDPJIT
195  numMatvecs = 0 ;
196  // NEED to grab the eignvectors from the named buffer here
198  {
202  int N = Layout::sitesOnNode()*Nc*Ns ;
203  int VectorSpaceSize = Nc*Ns*(A->subset()).numSiteTable();
204  EigInfo.init(invParam.Neig_max, N, VectorSpaceSize) ;
205  EigInfo.restartTol = invParam.restartTol.elem().elem().elem().elem();
206  if(invParam.file.read){
207  QDPIO::cout<<"MdagMSysSolverOptEigCG : reading evecs from disk"<<std::endl ;
208  QIOReadOptEvecs() ;
209  }
210  }
211 #endif
212  }
213 
214  //! Destructor is automatic
216  {
217  if(invParam.file.write){
218  QDPIO::cout<<"MdagMSysSolverOptEigCG : writing evecs to disk"<<std::endl ;
219  QIOWriteOptEvecs() ;
220  }
222  {
224  }
225  }
226 
227  //! Return the subset on which the operator acts
228  const Subset& subset() const {return A->subset();}
229 
230  //! Solver the linear system
231  /*!
232  * \param psi solution ( Modify )
233  * \param chi source ( Read )
234  * \return syssolver results
235  *
236  * Definitions supplied in the correspond .cc file
237  */
239 
240 
241  //! Solve the linear system starting with a chrono guess
242  /*!
243  * \param psi solution (Write)
244  * \param chi source (Read)
245  * \param predictor a chronological predictor (Read)
246  * \return syssolver results
247  */
248 
250  AbsChronologicalPredictor4D<T>& predictor) const
251  {
252 
253  START_CODE();
254 
255  // This solver uses InvCG2, so A is just the matrix.
256  // I need to predict with A^\dagger A
257  {
259  predictor(psi, (*MdagM), chi);
260  }
261  // Do solve
262  SystemSolverResults_t res=(*this)(psi,chi);
263 
264  // Store result
265  predictor.newVector(psi);
266  END_CODE();
267  return res;
268  }
269 
270 
271  private:
272 
273  // Hide default constructor
279  };
280 
281 } // End namespace
282 
283 
284 
285 #endif
286 
Abstract interface for a Chronological Solution predictor.
virtual void newVector(const T &psi)=0
Class for counted reference semantics.
Definition: handle.h:33
void init(int ldh, int lde_)
Definition: containers.h:134
Linear Operator.
Definition: linearop.h:27
M^dag.M linear operator.
Definition: lmdagm.h:22
Solve a M*psi=chi linear system by CG2 with eigenvectors.
SystemSolverResults_t operator()(T &psi, const T &chi) const
Solver the linear system.
void QIOReadOptEvecs()
Read a OptEigInfo Type
SystemSolverResults_t operator()(T &psi, const T &chi, AbsChronologicalPredictor4D< T > &predictor) const
Solve the linear system starting with a chrono guess.
MdagMSysSolverOptEigCG(Handle< LinearOperator< T > > A_, const SysSolverOptEigCGParams &invParam_)
Constructor.
const Subset & subset() const
Return the subset on which the operator acts.
void QIOWriteOptEvecs()
Write out an OptEigInfo Type
~MdagMSysSolverOptEigCG()
Destructor is automatic.
Handle< LinearOperator< T > > MdagM
Handle< LinearOperator< T > > A
SystemSolver disambiguator.
static T & Instance()
Definition: singleton.h:432
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
std::string uniqueId()
Return a unique id.
Definition: unique_id.cc:18
Class for counted reference semantics.
Handle< MapObject< int, EVPair< LatticeColorVector > > > obj
Linear Operators.
M^dag*M composition of a linear operator.
Named object support.
Named object function std::map.
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
push(xml_out,"Condensates")
LinOpSysSolverMGProtoClover::T T
int i
Definition: pbg5p_w.cc:55
multi1d< LatticeFermion > chi(Ncb)
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
START_CODE()
struct Chroma::SysSolverOptEigCGParams::File_t file
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Linear system solvers.
Solve a CG1 system.
Disambiguator for LinOp system solvers.
Generate a unique id.