CHROMA
ilu2prec_s_cprec_t_wilson_linop_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Unpreconditioned Wilson linear operator
3  */
4 
5 #include "qdp_config.h"
6 #if QDP_NS == 4
7 #if QDP_ND == 4
8 #if QDP_NC == 3
9 
10 #include "chromabase.h"
13 
14 using namespace QDP::Hints;
15 
16 namespace Chroma
17 {
18 
19 
20  //! Creation routine with Anisotropy
21  /*!
22  * \param u_ gauge field (Read)
23  * \param Mass_ fermion kappa (Read)
24  * \param aniso anisotropy struct (Read)
25  */
26  void ILU2PrecSCprecTWilsonLinOp::create(Handle< FermState<T,P,Q> > fs_,
27  const Real& Mass_,
28  const AnisoParam_t& anisoParam_)
29  {
30 #ifndef QDP_IS_QDPJIT
31  START_CODE();
32 
33  // Check we are in 4D
34  if ( Nd != 4 ) {
35  QDPIO::cout << "This class (ILU2PrecSCprecTWilsonLinOp) only works in 4D" << std::endl;
36  QDP_abort(1);
37  }
38 
39  // Check Aniso Direction -- has to be 3.
40  if ( anisoParam_.t_dir != 3 ) {
41  QDPIO::cout << "This class (ILU2PrecSCprecTWilsonLinOp) is hardwired for t_dir=3"<< std::endl;
42  QDP_abort(1);
43  }
44 
45  // Check that the time direction is local:
46  const multi1d<int>& s_size = QDP::Layout::subgridLattSize(); // Local Lattice
47  const multi1d<int>& t_size = QDP::Layout::lattSize(); // Total Latt Size
48  if( t_size[3] != s_size[3] ) {
49  QDPIO::cout << "This class (ILU2PrecSCprecTWilsonLinOp) needs time to be local" << std::endl;
50  QDP_abort(1);
51  }
52 
53  // Store gauge state etc
54  fs = fs_;
55  aniso = anisoParam_;
56  Mass = Mass_;
57 
58  // Work out aniso factors
59  Real ff = where(aniso.anisoP, aniso.nu / aniso.xi_0, Real(1));
60  fact = 1 + (Nd-1)*ff + Mass;
61  invfact = Real(1)/fact;
62  u = fs->getLinks();
63 
64  // Incorporate into links
65  if (aniso.anisoP)
66  {
67  // Rescale the u fields by the anisotropy
68  for(int mu=0; mu < u.size(); ++mu)
69  {
70  if (mu != aniso.t_dir)
71  u[mu] *= ff;
72  }
73  }
74 
75  // These will be useful for controlling loops.
76  int t_index=3;
77  int Nx = QDP::Layout::subgridLattSize()[0];
78  int Ny = QDP::Layout::subgridLattSize()[1];
79  int Nz = QDP::Layout::subgridLattSize()[2];
80  int Nt = QDP::Layout::subgridLattSize()[3];
81 
82  // Resize lookup table: for spatial site z,y,x we will get back an array of length t giving the site indices
83  // for the t sites, with that spatial index.
84  //
85  int Nspaceby2 = Nx*Ny*Nz/2; // This always works because Nx has to be even...
86 
87  tsite.resize(2, Nspaceby2, Nt);
88 
89 
90  int ssite_even=0;
91  int ssite_odd=0;
92 
93  int color;
94  int site;
95 
96  // Loop over all spatial sites
97  for(int z=0; z < Nz; z++) {
98  for(int y=0; y < Ny; y++) {
99  for(int x=0; x < Nx; x++) {
100  // Assemble the site indices for all timeslices of this coordinate
101 
102  multi1d<int> coords(Nd);
103  coords[0]=x;
104  coords[1]=y;
105  coords[2]=z;
106  color = (x + y + z) & 1;
107 
108  if ( color == 0 ) {
109  for(int t=0; t < Nt; t++) {
110  coords[3] = t;
111  tsite(0,ssite_even, t) = QDP::Layout::linearSiteIndex(coords);
112  }
113  ssite_even++;
114 
115  }
116  else {
117  for(int t=0; t < Nt; t++) {
118  coords[3] = t;
119  tsite(1,ssite_odd, t) = QDP::Layout::linearSiteIndex(coords);
120  }
121  ssite_odd++;
122  }
123  }
124  }
125  }
126 
127  // Figure out whether we are Schroedinger in time:
128  const FermBC<T,P,Q>& fbc = getFermBC();
129  if( fbc.nontrivialP() ) {
130  try {
131  const SchrFermBC& schr_ref = dynamic_cast<const SchrFermBC&>(fbc);
132  schrTP = ( schr_ref.getDir() == tDir() );
133  }
134  catch(std::bad_cast) {
135  schrTP = false;
136  }
137 
138  }
139  else {
140  // fbc are not nontrivial
141  schrTP = false;
142  }
143 
144  logDetTSq = zero;
145 
146  if( !schrTP ) {
147  // P and P_mat dag are needed for the Woodbury
148  // (P_mat for inverting T, P_mat_dag for inverting T_dag - NB P_mat_dag != (P_mat)^dagger
149  P_mat.resize(2, Nspaceby2, Nt);
150  P_mat_dag.resize(2, Nspaceby2, Nt);
151 
152  Q_mat_inv.resize(2,Nspaceby2);
153  Q_mat_dag_inv.resize(2,Nspaceby2);
154 
155  for(int cb3=0; cb3 < 2; cb3++) {
156  for(int site=0; site < Nspaceby2; site++) {
157 
158 
159  Real minvfact = Real(-1)*invfact;
160 
161  // Compute P by backsubsitution - See Balint's notes eq 38-42
162  P_mat(cb3, site, Nt-1) = u[t_index].elem( tsite(cb3,site,Nt-1) );
163  P_mat(cb3, site, Nt-1) *= minvfact.elem();
164 
165 
166  for(int t=Nt-2; t >=0; t--) {
167  P_mat(cb3, site, t) = u[t_index].elem( tsite(cb3, site,t) ) * P_mat(cb3, site,t+1);
168  P_mat(cb3, site, t) *= invfact.elem();
169  }
170 
171  // Compute the dagger. Opposite order (forward sub) similara to eq 38-42
172  // in notes
173  P_mat_dag(cb3, site, 0) = adj( u[t_index].elem( tsite(cb3, site, Nt-1) ) );
174  P_mat_dag(cb3, site, 0) *= minvfact.elem();
175 
176  for(int t=1; t < Nt; t++) {
177  P_mat_dag(cb3, site,t) = adj( u[t_index].elem( tsite(cb3, site, t-1) ) )*P_mat_dag(cb3,site, t-1) ;
178  P_mat_dag(cb3, site,t) *= invfact.elem();
179  }
180 
181 
182  // Compute Q = (1 + W^\dag P)^{-1}, W = [1, 0, 0, 0...] => (1 + P_{0})^{-1} eq: 43
183  // NB: This is not necessarily SU(3) now, so we can't just take the dagger to get the inverse
184  Real one=Real(1);
185  CMat one_plus_P0 = one.elem() + P_mat(cb3, site,0);
186  CentralTPrecNoSpinUtils::invert3by3( Q_mat_inv(cb3, site), one_plus_P0 );
187 
188  // Compute Q_dag = 1 + W^\dag P^\dag, W = [ 0, ..., 1 ] = > Q_dag = P^\dag [Nt-1]
189  // Similar to eq 43
190  // NB: This is not necessarily SU(3) now, so we can't just take the dagger to get the inverse
191  CMat one_plus_P0_dag = one.elem() + P_mat_dag(cb3, site,Nt-1);
192  CentralTPrecNoSpinUtils::invert3by3( Q_mat_dag_inv(cb3, site), one_plus_P0_dag);
193 
194  CMat prod = one_plus_P0_dag*one_plus_P0;
195  logDetTSq += CentralTPrecNoSpinUtils::logDet(prod);
196  }
197  }
198  }
199  else {
200  logDetTSq = 0;
201  }
202  Dw3D.create(fs_, anisoParam_);
203  END_CODE();
204 #endif
205  }
206 
207 } // End Namespace Chroma
208 
209 #endif
210 #endif
211 #endif
Support for time preconditioning.
Primary include file for CHROMA library code.
#define END_CODE()
Definition: chromabase.h:65
#define START_CODE()
Definition: chromabase.h:64
int mu
Definition: cool.cc:24
Real Mass
int y
Definition: meslate.cc:35
int z
Definition: meslate.cc:36
int x
Definition: meslate.cc:34
int t
Definition: meslate.cc:37
Nd
Definition: meslate.cc:74
multi1d< int > coords(const int x, const int y, const int z, const int t)
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static int tDir()
Definition: coulgauge.cc:14
static multi1d< LatticeColorMatrix > u
Double one
Definition: invbicg.cc:105
Double zero
Definition: invbicg.cc:106