CHROMA
t_wlinvcg.cc
Go to the documentation of this file.
1 
2 #include <iostream>
3 #include <cstdio>
4 
5 #include "chroma.h"
6 #include "primitives.h" // GTF: for PLUS
7 
8 
9 
10 using namespace Chroma;
11 
12 
13 #include "chromabase.h"
15 
16 //! Conjugate-Gradient (CGNE) algorithm for a generic Linear Operator
17 /*! \ingroup invert
18  * This subroutine uses the Conjugate Gradient (CG) algorithm to find
19  * the solution of the set of linear equations
20  *
21  * Chi = A . Psi
22  *
23  * where A = M^dag . M
24  *
25  * Algorithm:
26 
27  * Psi[0] := initial guess; Linear interpolation (argument)
28  * r[0] := Chi - M^dag . M . Psi[0] ; Initial residual
29  * p[1] := r[0] ; Initial direction
30  * IF |r[0]| <= RsdCG |Chi| THEN RETURN; Converged?
31  * FOR k FROM 1 TO MaxCG DO CG iterations
32  * a[k] := |r[k-1]|**2 / <Mp[k],Mp[k]> ;
33  * Psi[k] += a[k] p[k] ; New solution std::vector
34  * r[k] -= a[k] M^dag . M . p[k] ; New residual
35  * IF |r[k]| <= RsdCG |Chi| THEN RETURN; Converged?
36  * b[k+1] := |r[k]|**2 / |r[k-1]|**2 ;
37  * p[k+1] := r[k] + b[k+1] p[k]; New direction
38  *
39  * Arguments:
40  *
41  * \param M Linear Operator (Read)
42  * \param chi Source (Read)
43  * \param psi Solution (Modify)
44  * \param RsdCG CG residual accuracy (Read)
45  * \param MaxCG Maximum CG iterations (Read)
46  * \param n_count Number of CG iteration (Write)
47  *
48  * Local Variables:
49  *
50  * p Direction std::vector
51  * r Residual std::vector
52  * cp | r[k] |**2
53  * c | r[k-1] |**2
54  * k CG iteration counter
55  * a a[k]
56  * b b[k+1]
57  * d < p[k], A.p[k] >
58  * Mp Temporary for M.p
59  *
60  * Subroutines:
61  * +
62  * A Apply matrix M or M to std::vector
63  *
64  * Operations:
65  *
66  * 2 A + 2 Nc Ns + N_Count ( 2 A + 10 Nc Ns )
67  */
68 
69 void WlInvCG2(const LinearOperator& M,
70  const LatticeFermion& chi,
71  LatticeFermion& psi,
72  const Real& RsdCG,
73  int MaxCG,
74  int& n_count)
75 {
76  const Subset& s = M.subset();
77 
78  LatticeFermion mp;
79  Real a;
80  Real b;
81  Double c;
82  Double d;
83 
84  Real rsd_sq = (RsdCG * RsdCG) * Real(norm2(chi,s));
85 
86  // +
87  // r[0] := Chi - A . Psi[0] where A = M . M
88 
89  // +
90  // r := [ Chi - M(u) . M(u) . psi ]
91  LatticeFermion r;
92  r[s] = chi - M(M(psi, PLUS), MINUS);
93 
94  // p[1] := r[0]
95  LatticeFermion p;
96  p[s] = r;
97 
98  // Cp = |r[0]|^2
99  Double cp = norm2(r, s); /* 2 Nc Ns flops */
100 
101  QDPIO::cout << "WlInvCG: k = 0 cp = " << cp << " rsd_sq = " << rsd_sq << std::endl;
102 
103  // IF |r[0]| <= RsdCG |Chi| THEN RETURN;
104  if ( toBool(cp <= rsd_sq) )
105  {
106  n_count = 0;
107  return;
108  }
109 
110  //
111  // FOR k FROM 1 TO MaxCG DO
112  //
113  for(int k = 1; k <= MaxCG; ++k)
114  {
115  // c = | r[k-1] |**2
116  c = cp;
117 
118  // a[k] := | r[k-1] |**2 / < p[k], Ap[k] > ;
119  // +
120  // First compute d = < p, A.p > = < p, M . M . p > = < M.p, M.p >
121  // Mp = M(u) * p
122  mp[s] = M(p, PLUS);
123 
124  // d = | mp | ** 2
125  d = norm2(mp, s); /* 2 Nc Ns flops */
126 
127  a = Real(c)/Real(d);
128 
129  // Psi[k] += a[k] p[k]
130  psi[s] += a * p; /* 2 Nc Ns flops */
131 
132  // r[k] -= a[k] A . p[k] ;
133  // + +
134  // r = r - M(u) . Mp = M . M . p = A . p
135  r[s] -= a * M(mp, MINUS);
136 
137  // IF |r[k]| <= RsdCG |Chi| THEN RETURN;
138 
139  // cp = | r[k] |**2
140  cp = norm2(r, s); /* 2 Nc Ns flops */
141 
142  QDPIO::cout << "WlInvCG: k = " << k << " cp = " << cp << std::endl;
143 
144  if ( toBool(cp <= rsd_sq) )
145  {
146  n_count = k;
147  return;
148  }
149 
150  // b[k+1] := |r[k]|**2 / |r[k-1]|**2
151  b = Real(cp) / Real(c);
152 
153  // p[k+1] := r[k] + b[k+1] p[k]
154  p[s] = r + b*p; /* Nc Ns flops */
155  }
156  n_count = MaxCG;
157  QDP_error_exit("too many CG iterations: count = %d", n_count);
158 }
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 int main(int argc, char **argv)
169 {
170  // Put the machine into a known state
171  Chroma::initialize(&argc, &argv);
172 
173  // Setup the layout
174  const int foo[] = {2,2,2,2};
175  multi1d<int> nrow(Nd);
176  nrow = foo; // Use only Nd elements
177  Layout::setLattSize(nrow);
178  Layout::create();
179 
180  XMLFileWriter xml("t_lwldslash.xml");
181  push(xml,"t_lwldslash");
182 
183  //! Test out dslash
184  multi1d<LatticeColorMatrix> u(Nd);
185  for(int m=0; m < u.size(); ++m)
186  gaussian(u[m]);
187 
188  LatticeFermion psi, chi;
189  random(psi);
190  chi = zero;
191 
192 
193  // Construct szin-like gauge field
194  multi3d<ColorMatrix> u_tmp(Nd,2,Layout::sitesOnNode()/2);
195  for(int m=0; m < u.size(); ++m)
196  {
197  multi1d<ColorMatrix> u_tt(Layout::sitesOnNode());
198  QDP_extract(u_tt, u[m], all);
199 
200  // Use this guy
201  for(int i=0; i < u_tt.size(); ++i)
202  {
203  int ii = i >> 1;
204  int cb = i / (Layout::sitesOnNode()/2);
205  u_tmp[m][cb][ii] = transpose(u_tt[ii]);
206  }
207  }
208 
209 
210 
211  pack_gauge_field(u_tmp);
212 
213  //! Create a linear operator
214  WilsonDslash D(u);
215  chi = D(psi, PLUS, 0);
216 
217  write(xml,"Nd", Nd);
218  write(xml,"Nc", Nc);
219  write(xml,"Ns", Ns);
220  write(xml,"nrow", nrow);
221  write(xml,"psi", psi);
222  write(xml,"chi", chi);
223 
224  //! Create and try a more sophisticated operator
225  Real Kappa = 0.1;
226  PreconditionedWilson M(u,Kappa);
227  LatticeFermion eta;
228  eta = M(psi, PLUS);
229 
230  write(xml,"eta", eta);
231 
232  pop(xml);
233 
234  // Time to bolt
236 
237  exit(0);
238 }
Primary include file for CHROMA in application codes.
Primary include file for CHROMA library code.
Linear Operator.
Definition: linearop.h:27
virtual const Subset & subset() const =0
Return the subset on which the operator acts.
General Wilson-Dirac dslash.
Definition: lwldslash_w.h:48
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void WlInvCG2(const LinearOperator &M, const LatticeFermion &chi, LatticeFermion &psi, const Real &RsdCG, int MaxCG, int &n_count)
Conjugate-Gradient (CGNE) algorithm for a generic Linear Operator.
Definition: t_wlinvcg.cc:69
Conjugate-Gradient algorithm for a generic Linear Operator.
static int m[4]
Definition: make_seeds.cc:16
Nd
Definition: meslate.cc:74
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > & RsdCG
Definition: pbg5p_w.cc:30
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
static multi1d< LatticeColorMatrix > u
int n_count
Definition: invbicg.cc:78
void transpose(multi2d< LatticeColorVector > &dist_rep, const multi2d< LatticeColorVector > &prop_rep)
Take transpose of a matrix in (explicit) spin space.
push(xml_out,"Condensates")
Double c
Definition: invbicg.cc:108
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
Definition: pbg5p_w.cc:32
Real rsd_sq
Definition: invbicg.cc:121
int i
Definition: pbg5p_w.cc:55
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
Double cp
Definition: invbicg.cc:107
multi1d< LatticeFermion > chi(Ncb)
Complex a
Definition: invbicg.cc:95
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
LatticeFermion eta
Definition: mespbg5p_w.cc:37
DComplex d
Definition: invbicg.cc:99
multi1d< LatticeFermion > mp(Ncb)
Complex b
Definition: invbicg.cc:96
int cb
Definition: invbicg.cc:120
Double zero
Definition: invbicg.cc:106
int k
Definition: invbicg.cc:119
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
int main(int argc, char **argv)
Definition: t_wlinvcg.cc:168