CHROMA
asqtad_cps_wrapper_qprop.cc
Go to the documentation of this file.
1 #include "chromabase.h"
2 
6 
7 
8 extern "C" {
9 
10 // double precision (recommended)
11 #define LEVEL3_PREC 2
12 
13 #if LEVEL3_PREC == 2
14 // double precision
15 #define QLA_Precision 'D'
16 #define QOP_Precision 2
17 #else
18 // single precision
19 #define QLA_Precision 'F'
20 #define QOP_Precision 1
21 #endif
22 
23 
24 #define QLA_Colors 3
25 
26 #include <qdp.h>
27 #include <qdp_types.h>
28 #include <qdp_common.h>
29 
30 #if LEVEL3_PREC == 2
31 #include <qdp_d3_generic.h>
32 #include <qdp_d3.h>
33 #else
34 #include <qdp_f3_generic.h>
35 #include <qdp_f3.h>
36 #endif
37 
38 #include <qop_qdp.h>
39 #include <qmp.h>
40 
41 #include <qla.h>
42 #include <qla_complex.h>
43 
44 
45 }
46 
47 
48 
49 
50 /**************************************/
51 /** level 3 interface code **/
52 
53 
54 //
55 // utility routines
56 //
57 
58 //
59 // from test_common.c
60 //
61 
62 QLA_Real
63 get_plaq(QDP_ColorMatrix *link[])
64 {
65  int mu, nu;
66  QLA_Real plaq;
67  QDP_ColorMatrix *temp1, *temp2, *temp3, *temp4;
68 
69 #ifdef LOCAL_SUM
70  QDP_Real *treal1, *treal2;
71  treal1 = QDP_create_R();
72  treal2 = QDP_create_R();
73  QDP_R_eq_zero(treal2, QDP_all);
74 #else
75  QLA_Real tplaq;
76  plaq = 0;
77 #endif
78 
79  temp1 = QDP_create_M();
80  temp2 = QDP_create_M();
81  temp3 = QDP_create_M();
82  temp4 = QDP_create_M();
83 
84  // http://usqcd.jlab.org/usqcd-docs/qdp/qdpc.html/Functions-involving-shifts.html#Functions-involving-shifts
85 
86  for(mu=0; mu<QDP_ndim()-1; ++mu) {
87  for(nu=mu+1; nu<QDP_ndim(); ++nu) {
88 
89  QDP_M_eq_sM(temp1, link[nu], QDP_neighbor[mu], QDP_forward, QDP_all);
90  QDP_M_eq_sM(temp2, link[mu], QDP_neighbor[nu], QDP_forward, QDP_all);
91 
92  QDP_M_eq_Ma_times_M(temp3, link[nu], link[mu], QDP_all);
93 
94  QDP_M_eq_M_times_M(temp4, temp3, temp1, QDP_all);
95  QDP_discard_M(temp1);
96 
97 #ifdef LOCAL_SUM
98  QDP_R_eq_re_M_dot_M(treal1, temp2, temp4, QDP_all);
99  QDP_discard_M(temp2);
100  QDP_R_peq_R(treal2, treal1, QDP_all);
101 #else
102  QDP_r_eq_re_M_dot_M(&tplaq, temp2, temp4, QDP_all);
103  QDP_discard_M(temp2);
104  plaq += tplaq;
105 #endif
106 
107  }
108  }
109 
110 #ifdef LOCAL_SUM
111  QDP_r_eq_sum_R(&plaq, treal2, QDP_all);
112  QDP_destroy_R(treal1);
113  QDP_destroy_R(treal2);
114 #endif
115 
116  QDP_destroy_M(temp1);
117  QDP_destroy_M(temp2);
118  QDP_destroy_M(temp3);
119  QDP_destroy_M(temp4);
120 
121  // normalisation of plaquette <p> = 1 for cold config.
122  return plaq/(3.0*0.5*QDP_ndim()*(QDP_ndim()-1)*QDP_volume());
123 }
124 
125 
126 
127 //
128 // Convert the gaugefield from chroma/qdp++ format
129 // into qdp format.
130 //
131 QDP_ColorMatrix ** QDP_create_gauge_from_chroma (
132  multi1d<LatticeColorMatrix> & u_chroma)
133 {
134  QDP_ColorMatrix **u;
135 
136  const int ndim = 4 ;
137  u = (QDP_ColorMatrix **) malloc(ndim*sizeof(QDP_ColorMatrix *));
138  for(int i=0; i<ndim; i++) u[i] = QDP_create_M();
139  // more work --- check error codes
140 
141  QLA_ColorMatrix *tmp;
142  tmp =
143  (QLA_ColorMatrix *) malloc(QDP_sites_on_node*sizeof(QLA_ColorMatrix));
144 
145  for(int i=0; i<ndim; i++)
146  {
147  for(int site=0 ; site < QDP_sites_on_node ; ++site)
148  {
149  for(int ic = 0 ; ic < 3 ; ++ic)
150  for(int jc = 0 ; jc < 3 ; ++jc)
151  {
152 
153  Real rrr = u_chroma[i].elem(site).elem().elem(ic,jc).real() ;
154  Real iii = u_chroma[i].elem(site).elem().elem(ic,jc).imag() ;
155 
156  QLA_Real zre = toFloat(rrr) ;
157  QLA_Real zim = toFloat(iii) ;
158  QLA_Complex z ;
159  QLA_C_eq_R_plus_i_R(&z,&zre,&zim) ;
160 
161  // QDPIO::cout << "" << ic << " , " << jc << " = " << zre << " " << zim << std::endl ;
162  QLA_M_eq_elem_C(&tmp[site],&z,ic,jc) ;
163  }
164 
165  // tmp[site]
166  }
167  // convert to qdp format
168  QDP_insert_M(u[i], tmp,QDP_all) ;
169  }
170 
171  free(tmp) ;
172 
173 
174  // convert to QOP format.
175  QLA_Real plaq;
176  plaq = get_plaq(u);
177  QDPIO::cout << "QDP_create_gauge_from_chroma:: plaquette = " << plaq << "\n";
178 
179  return u ;
180 }
181 
182 //
183 // Convert the colorstd::vector from chroma/qdp++ format
184 // into qdp format.
185 //
186 void convert_chroma_to_qdp(QDP_ColorVector *out,
187 const LatticeStaggeredFermion & in )
188 {
189  QLA_ColorVector *tmp;
190  tmp =
191  (QLA_ColorVector *) malloc(QDP_sites_on_node*sizeof(QLA_ColorVector));
192 
193  for(int site=0 ; site < QDP_sites_on_node ; ++site)
194  {
195 
196  for(int ic = 0 ; ic < 3 ; ++ic)
197  {
198  QLA_Real zre ;
199 
200  Real rrr = in.elem(site).elem(0).elem(ic).real() ;
201  Real iii = in.elem(site).elem(0).elem(ic).imag() ;
202 
203  zre = toFloat(rrr) ;
204  QLA_Real zim = toFloat(iii) ;
205  QLA_Complex z ;
206  QLA_C_eq_R_plus_i_R(&z,&zre,&zim) ;
207 
208  QLA_V_eq_elem_C(&tmp[site],&z,ic) ;
209 
210  // QDPIO::cout << "convert_chroma_to_qdp-DEBUG " << site << " " << ic << " " << zre << " " << zim << "\n" ;
211  }
212 
213  } // loop over sites
214 
215  // convert to qdp format
216  QDP_insert_V(out, tmp,QDP_all) ;
217 
218 
219  free(tmp) ;
220 
221 }
222 
223 
224 
225 //
226 // Convert the colorstd::vector in qdp into chroma/qdp++ format
227 //
228 //
229 void convert_qdp_to_chroma(LatticeStaggeredFermion & out,
230  QDP_ColorVector *in)
231 {
232  Complex zz ;
233  QLA_Complex z ;
234 
235  QLA_ColorVector *tmp;
236  tmp = QDP_expose_V(in) ;
237  // this is done at the pointer level
238 
239  for(int site=0 ; site < QDP_sites_on_node ; ++site)
240  {
241  for(int ic = 0 ; ic < 3 ; ++ic)
242  {
243  QLA_C_eq_elem_V(&z,&tmp[site],ic) ;
244 
245  QLA_Real zre ;
246  zre = QLA_real(z) ;
247 
248  QLA_Real zim ;
249  zim = QLA_imag(z) ;
250 
251  Real zre_chroma , zim_chroma ;
252  zre_chroma = (Real) zre ;
253  zim_chroma = (Real) zim ;
254 
255  // zre_chroma *= (4 * 0.03) ;
256  // zim_chroma *= (4 * 0.03) ;
257 
258  zre_chroma /= 2.0 ;
259  zim_chroma /= 2.0 ;
260 
261  out.elem(site).elem(0).elem(ic).real() = toFloat(zre_chroma) ;
262  out.elem(site).elem(0).elem(ic).imag() = toFloat(zim_chroma) ;
263 
264 
265  // QDPIO::cout << "convert_qdp_to_chroma-DEBUG " << site << " " << ic << " " << zre_chroma << " " << zim_chroma << "\n" ;
266 
267  } // loop over color
268 
269  } // loop over sites
270 
271 }
272 
273 
274 //
275 // set the coefficients from the MILC code.
276 // This routine came from
277 // generic_ks/quark_stuff.c
278 // generic_ks/imp_actions/asqtad_action.h
279 //
280 
281 void load_qop_asqtad_coeffs(QOP_asqtad_coeffs_t *c, QLA_Real weight)
282 {
283  QLA_Real ferm_epsilon;
284 #define MAX_BASIC_PATHS 6
285  static QLA_Real act_path_coeff[MAX_BASIC_PATHS] = {
286  ( 1.0/8.0)+(6.0/16.0)+(1.0/8.0), /* one link */
287  /*One link is 1/8 as in fat7 +3/8 for Lepage + 1/8 for Naik */
288  (-1.0/24.0), /* Naik */
289  (-1.0/8.0)*0.5, /* simple staple */
290  ( 1.0/8.0)*0.25*0.5, /* displace link in two directions */
291  (-1.0/8.0)*0.125*(1.0/6.0), /* displace link in three directions */
292  (-1.0/16 ), /* Correct O(a^2) errors */
293  };
294 
295 
296  /* Load path coefficients from table */
297  // act_path_coeff = get_quark_path_coeff();
298 
299  ferm_epsilon = weight;
300 
301  /* Path coefficients times fermion epsilon */
302 
303  c->one_link = act_path_coeff[0]*ferm_epsilon ;
304  c->naik = act_path_coeff[1]*ferm_epsilon ;
305  c->three_staple = act_path_coeff[2]*ferm_epsilon ;
306  c->five_staple = act_path_coeff[3]*ferm_epsilon ;
307  c->seven_staple = act_path_coeff[4]*ferm_epsilon ;
308  c->lepage = act_path_coeff[5]*ferm_epsilon ;
309 }
310 
311 
312 
313 /*** end of level 3 interface code ****/
314 /**************************************/
315 
316 
317 namespace Chroma
318 {
319 
320  //! Constructor
321  /*!
322  // Keeping the same interface as for the ordinary staggered
323  // qprop...
324  //
325  // But the M_ and A_ linop handles are no longer kept
326  // (are ignored) -- is there a nice way around this ?
327  // Perhaps not
328  */
329 
330  // internal variables to this file
331 static QDP_ColorVector *out ;
332 static QDP_ColorVector *in ;
333 
334 
335  static QOP_info_t info; // wot for ????
336  static QOP_FermionLinksAsqtad *fla ;
337 
338  // -----------------------------------------------------
339 
341  const EvenOddStaggeredTypeFermAct<LatticeStaggeredFermion,
342  multi1d<LatticeColorMatrix>,
343  multi1d<LatticeColorMatrix> >& S_,
345  const SysSolverCGParams& invParam_) :
346  // invParam(invParam_),Mass(S_.getQuarkMass()), state(state_)
347  invParam(invParam_),Mass(S_.getQuarkMass()),
348  state(state_.cast<AsqtadConnectStateBase>()) , M(S_.linOp(state_))
349  {
350  // Here is how to get at the gauge links: (Thanks Balint).
351  // gauge not needed
352  const multi1d<LatticeColorMatrix>& u = state->getLinks();
353 
354  const multi1d<LatticeColorMatrix>& u_fat = state->getFatLinks();
355  const multi1d<LatticeColorMatrix>& u_triple = state->getTripleLinks();
356 
357  multi1d<LatticeColorMatrix> u_chroma(Nd); // hack
358 
359  QDP_ColorMatrix **u_fat_qdp;
360  u_chroma = u_fat ;
361  u_fat_qdp = QDP_create_gauge_from_chroma(u_chroma) ;
362  // u_fat_qdp = QDP_create_gauge_from_chroma(u_fat) ;
363 
364  QDP_ColorMatrix **u_triple_qdp;
365  u_chroma = u_triple ;
366  u_triple_qdp = QDP_create_gauge_from_chroma(u_chroma) ;
367  // u_triple_qdp = QDP_create_gauge_from_chroma(u_triple) ;
368 
369  fla = QOP_asqtad_convert_L_from_qdp(u_fat_qdp,u_triple_qdp) ;
370 
371 #if 0
372  for(int i=0; i<Nd ; i++) QDP_destroy_M(u_fat_qdp[i]);
373  free(u_fat_qdp) ;
374 
375 
376  for(int i=0; i<Nd; i++) QDP_destroy_M(u_triple_qdp[i]);
377  free(u_triple_qdp) ;
378 
379 #endif
380 
381 
382 #if 0
383  multi1d<LatticeColorMatrix> u_with_phases(Nd);
384  state.getFermBC().modify(u_with_phases);
385 #endif
386 
387 #if 0
388  // add staggered phases to the chroma gauge field
389  multi1d<LatticeColorMatrix> u_chroma(Nd);
390  // alpha comes from the StagPhases:: namespace
391  for(int i = 0; i < Nd; i++) {
392  u_chroma[i] = u[i] ;
393  u_chroma[i] *= StagPhases::alpha(i);
394  }
395 
396 
397  // setup the qdp
398  QOP_asqtad_coeffs_t coeffs;
399  load_qop_asqtad_coeffs(&coeffs, 1.0) ;
400  Real u0 = 1.0 ;
401 
402  // ---- create the fat links ---------
403  QDP_ColorMatrix **uqdp ;
404  uqdp = QDP_create_gauge_from_chroma (u_chroma) ;
405 
406  QOP_GaugeField *gf;
407  gf = QOP_convert_G_from_qdp(uqdp);
408 
409  // create the fat links
410  fla = QOP_asqtad_create_L_from_G(&info, &coeffs, gf);
411 
412  // remove the space for the gauge configuration
413  // gf contains pointers to 4 components of uqdp,
414  // so uqdp[i] do not need to be destroyed as well
415  QOP_destroy_G(gf) ;
416  free(uqdp) ;
417 
418 #endif
419 
420 }
421 
422 
423  //! Destructor may not well be automatic
425  {
426  QOP_asqtad_destroy_L(fla) ;
427 
428  }
429 
431  AsqtadCPSWrapperQprop::operator() (LatticeStaggeredFermion& psi,
432  const LatticeStaggeredFermion& q_source) const
433  {
435 
436 
437  out = QDP_create_V();
438  in = QDP_create_V();
439 
440 
441  // convert "q_source" to in (qdp format)
442  convert_chroma_to_qdp(in,q_source) ;
444 
445  QOP_ColorVector *qopout, *qopin;
446  qopout = QOP_create_V_from_qdp(out);
447  qopin = QOP_create_V_from_qdp(in);
448 
449  QOP_invert_arg_t inv_arg;
450  QOP_resid_arg_t res_arg;
451 
452  res_arg.rsqmin = toFloat(invParam.RsdCG)*
453  toFloat(invParam.RsdCG);
454 
455 
456  if( invParam.MaxCGRestart > 0 )
457  inv_arg.max_restarts = invParam.MaxCGRestart ;
458  else
459  inv_arg.max_restarts = 1;
460 
461  inv_arg.restart = invParam.MaxCG ;
462  inv_arg.max_iter = inv_arg.max_restarts * invParam.MaxCG;
463  inv_arg.evenodd = QOP_EVENODD ;
464 
465  QLA_Real mass = toFloat(Mass) ;
466 
467  QDPIO::cout << "level3 asqtad inverter mass = " << Mass ;
468  QDPIO::cout << " iters = " << invParam.MaxCG << " restarts= " ;
469  QDPIO::cout << inv_arg.max_restarts << "\n" ;
470 
471  // invert
472  QOP_verbose(QOP_VERB_HI);
473  QOP_asqtad_invert(&info, fla, &inv_arg, &res_arg, mass, qopout, qopin);
474 
475  QDPIO::cout << "QOP Inversion results\n" ;
476  QDPIO::cout << "QOP iters = " << res_arg.final_iter << "\n" ;
477  QDPIO::cout << "QOP final restart = " << res_arg.final_restart << "\n" ;
478  QDPIO::cout << "QOP ||r||^2 = " << res_arg.final_rsq << "\n" ;
479  QDPIO::cout << "QOP ||r|| = " << sqrt(res_arg.final_rsq) << "\n" ;
480  res.n_count = res_arg.final_iter ;
481  res.resid = res_arg.final_rsq; // check norm convention
482 
483  // convert out to psi (qdp++ format)
484  QOP_extract_V_to_qdp(out,qopout) ;
486 
487 
488  // unscale the norm
489 #if 0
490  exit(0) ;
491  Real mass_scale = 4.0 * Mass ;
492  // psi /= 4.0 * Mass ;
493  T tt;
494  tt = psi ;
495  psi = tt / mass_scale ;
496 #endif
497 
498  // Compute residual
499  {
500  T r;
501  (*M)(r, psi, PLUS);
502  r -= q_source ;
503  res.resid = sqrt(norm2(r));
504  QDPIO::cout << "AsqtadCPSWrapperQprop: true residual: " << res.resid << std::endl;
505  QDPIO::cout << "AsqtadCPSWrapperQprop: || q_source ||: " << sqrt(norm2(q_source)) << std::endl;
506 
507  QDPIO::cout << "AsqtadCPSWrapperQprop: true residual/source: " << res.resid/ sqrt(norm2(q_source)) << std::endl;
508 
509  }
510 
511  QOP_destroy_V(qopout);
512  QOP_destroy_V(qopin);
513 
514  QDP_destroy_V(out);
515  QDP_destroy_V(in);
516 
517  return res;
518  }
519 
520 }
521 
522 
523 void setup_levelthree(multi1d<int> nrow, int *argc, char ***argv )
524 {
525  // QDP_start(argc, argv);
526  int tmp = ::QDP_initialize(argc,argv ) ;
527 
528  QDPIO::cout << "Setting up the level 3 code\n" ;
529  QDPIO::cout << "----------------------------\n";
530 
531  int ndim = 4;
532  int *lattice_size;
533  lattice_size = (int *) malloc(ndim*sizeof(int));
534  for(int i=0 ; i < ndim ; ++i)
535  lattice_size[i] = nrow[i] ;
536 
537  QDPIO::cout << "QDP/QOP lattice volume = " ;
538  QDPIO::cout << lattice_size[0] << " " << lattice_size[1] << " " ;
539  QDPIO::cout << lattice_size[2] << " " << lattice_size[3] << "\n" ;
540 
541 #if LEVEL3_PREC == 2
542  QDPIO::cout << "DOUBLE PRECISION level3\n" ;
543 #else
544  QDPIO::cout << "SINGLE PRECISION level3\n" ;
545 #endif
546 
547 
548 
549  QDP_set_latsize(ndim, lattice_size);
550  QDP_create_layout();
551 
552  /** set up the level 3 code **/
553  QOP_layout_t qoplayout;
554  qoplayout.latdim = ndim;
555  qoplayout.latsize = (int *) malloc(ndim*sizeof(int));
556  for(int i=0; i<ndim; i++) {
557  qoplayout.latsize[i] = lattice_size[i];
558  }
559  qoplayout.machdim = -1;
560  QOP_init(&qoplayout);
561 
562 
563  free(lattice_size) ;
564 
565 }
566 
567 
void load_qop_asqtad_coeffs(QOP_asqtad_coeffs_t *c, QLA_Real weight)
void convert_chroma_to_qdp(QDP_ColorVector *out, const LatticeStaggeredFermion &in)
void setup_levelthree(multi1d< int > nrow, int *argc, char ***argv)
QLA_Real get_plaq(QDP_ColorMatrix *link[])
QDP_ColorMatrix ** QDP_create_gauge_from_chroma(multi1d< LatticeColorMatrix > &u_chroma)
void convert_qdp_to_chroma(LatticeStaggeredFermion &out, QDP_ColorVector *in)
#define MAX_BASIC_PATHS
Propagator solver for an even-odd non-preconditioned fermion operator.
Primary include file for CHROMA library code.
Handle< AsqtadConnectStateBase > state
SystemSolverResults_t operator()(LatticeStaggeredFermion &psi, const LatticeStaggeredFermion &chi) const
Solver the linear system.
~AsqtadCPSWrapperQprop()
Destructor is automatic.
Basic "Connect State" for ASQTAD.
Definition: asqtad_state.h:24
Even-odd preconditioned Staggered-like fermion actions.
Definition: fermact.orig.h:714
Support class for fermion actions and linear operators.
Definition: state.h:94
Class for counted reference semantics.
Definition: handle.h:33
Real u0
int mu
Definition: cool.cc:24
int nu
Definition: cool.cc:25
unsigned i
Definition: ldumul_w.cc:34
Double tmp
Definition: meslate.cc:60
int z
Definition: meslate.cc:36
int c
Definition: meslate.cc:61
Double plaq
Definition: meslate.cc:58
Nd
Definition: meslate.cc:74
static const LatticeInteger & alpha(const int dim)
Definition: stag_phases_s.h:43
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
int i
Definition: pbg5p_w.cc:55
Double mass
Definition: pbg5p_w.cc:54
@ PLUS
Definition: chromabase.h:45
static QOP_FermionLinksAsqtad * fla
LatticeFermion psi
Definition: mespbg5p_w.cc:35
static QDP_ColorVector * out
Constructor.
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
static QOP_info_t info
static QDP_ColorVector * in
Double link
Definition: pade_trln_w.cc:146
Params for CG inverter.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Solve a CG1 system.