Bridge++  Ver. 2.0.2
asolver_SAP_QWS.cpp
Go to the documentation of this file.
1 
10 //====================================================================
11 #include "asolver_SAP_QWS.h"
12 #include "lib_alt_QXS/Field/afield-inc.h"
13 namespace {
14 #ifndef AFILED_HAS_SUB
15  template<typename AFIELD>
16  inline void sub(AFIELD& v, const AFIELD& w)
17  {
18  typedef typename AFIELD::real_t real_t;
19  axpy(v, real_t(-1.0), w);
20  }
21 #endif
22 
23 #ifndef AFILED_HAS_ADD
24  template<typename AFIELD>
25  inline void add(AFIELD& v, const AFIELD& w)
26  {
27  typedef typename AFIELD::real_t real_t;
28  axpy(v, real_t(1.0), w);
29  }
30 #endif
31 }
32 
33 template<typename AFIELD>
34 const std::string ASolver_SAP_QWS<AFIELD>::class_name = "ASolver_SAP_QWS";
35 
36 #ifdef USE_QWSLIB
37 //====================================================================
38 template<typename AFIELD>
39 void ASolver_SAP_QWS<AFIELD>::init(void)
40 {
41  ThreadManager::assert_single_thread(class_name);
42 
43  int nin = m_fopr->field_nin();
44  int nvol = m_fopr->field_nvol();
45  int nex = m_fopr->field_nex();
46 
47  m_x.reset(nin, nvol, nex);
48 
49 #ifdef DEBUG
50  m_r.reset(nin, nvol, nex);
51  m_b.reset(nin, nvol, nex);
52 #endif
53 
54  { // working are for qws
55  int ret = 0;
56  void *p;
57  ret = posix_memalign(&p, CLS, sizeof(scs_t) * vols * 2);
58  m_s = (scs_t *)p;
59  ret = posix_memalign(&p, CLS, sizeof(scs_t) * vols * 2);
60  m_q = (scs_t *)p;
61  if (ret) {
62  vout.crucial(m_vl, "%s: allocation failed\n", class_name.c_str());
63  exit(EXIT_FAILURE);
64  }
65  }
66 
67  m_fopr->set_mode("D");
68 
69  m_nconv = -1;
70 }
71 
72 
73 //====================================================================
74 template<typename AFIELD>
75 void ASolver_SAP_QWS<AFIELD>::tidyup(void)
76 {
77  // free working area for qws
78  if (m_s) {
79  free(m_s);
80  }
81  if (m_q) {
82  free(m_q);
83  }
84  // ThreadManager::assert_single_thread(class_name);
85  //delete m_sap_minres;
86 }
87 
88 
89 //====================================================================
90 template<typename AFIELD>
91 void ASolver_SAP_QWS<AFIELD>::set_parameters(const Parameters& params)
92 {
93  const string str_vlevel = params.get_string("verbose_level");
94 
95  m_vl = vout.set_verbose_level(str_vlevel);
96 
97  //- fetch and check input parameters
98  int Niter, Nrestart;
99  double Stop_cond;
100 
101  int err = 0;
102  err += params.fetch_int("maximum_number_of_iteration", Niter);
103  err += params.fetch_int("maximum_number_of_restart", Nrestart);
104  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
105 
106  if (err) {
107  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
108  class_name.c_str());
109  exit(EXIT_FAILURE);
110  }
111 
112  int Niter2 = Niter * Nrestart;
113  set_parameters(Niter2, Stop_cond);
114 }
115 
116 
117 //====================================================================
118 template<typename AFIELD>
119 void ASolver_SAP_QWS<AFIELD>::set_parameters(const int Niter,
120  const real_t Stop_cond)
121 {
122  ThreadManager::assert_single_thread(class_name);
123 
124  m_Niter = Niter;
125  m_Stop_cond = Stop_cond;
126  std::string prec = "double";
127  if (sizeof(real_t) == 4) prec = "float";
128 
129  // N.B. Stop_cond is irrelevant
130  vout.general(m_vl, "%s:\n", class_name.c_str());
131  vout.general(m_vl, " Precision: %s\n", prec.c_str());
132  vout.general(m_vl, " Niter = %d\n", m_Niter);
133  vout.general(m_vl, " nm = %d\n", m_nm);
134  vout.general(m_vl, " Stop_cond is not used\n");
135 
136  // flop coung
137  int NPE = CommonParameters::NPE();
138  int Nx = CommonParameters::Nx();
139  int Ny = CommonParameters::Ny();
140  int Nz = CommonParameters::Nz();
141  int Nt = CommonParameters::Nt();
142  int Nx2 = Nx / 2;
143  int Nvol2 = Nx * Ny * Nz * Nt / 2;
144  double flop_clv = 576;
145  double flop_sap_hops
146  = 2 * ((Nx2 - 1) * Ny * Nz * Nt * 168 + Nx2 * (Ny - 1) * Nz * Nt * 168
147  + Nx2 * Ny * (Nz - 1) * Nt * 168 + Nx2 * Ny * Nz * (Nt - 1) * 156);
148  double flop_DEE_proc = (48 + flop_clv) * Nvol2 + flop_sap_hops; // per domain
149  double flop_DEO_proc // contains accum
150  = (2 * Nvol2 - 2 * (Nx2 - 2) * (Ny - 2) * (Nz - 2) * (Nt - 2)) * 624 // clov
151  + 2 * (Nx2 * Nz * Nt * 168 + Nx2 * Nz * Nt * 168 + Nx2 * Ny * Nt * 168 + Nx2 * Ny * Nz * 156);
152 
153  double flop_AEE_proc = m_nm * (flop_DEE_proc + 48 * Nvol2); // jinv_ddd
154  double flop
155  = m_Niter * (2 * flop_AEE_proc // jinv_ddd
156  + 2 * flop_DEE_proc // ddd_in
157  + 2 * flop_DEO_proc // ddd_out_pre, ddd_out_pos
158  + 2 * 2 * 12 * Nvol2); // accum_addsub
159  flop += 2 * flop_AEE_proc + flop_DEO_proc; // jinv_dd, ddd_out_pre, ddd_out_pos
160  flop += flop_clv * 2 * Nvol2; // clov_inv on the source
161  m_flop = flop * NPE;
162 }
163 
164 
165 //====================================================================
166 template<typename AFIELD>
167 void ASolver_SAP_QWS<AFIELD>::solve(AFIELD& x, const AFIELD& b,
168  int& Nconv, real_t& diff)
169 {
170  // vout.paranoiac(m_vl, "asolver_SAP_QWS: solve, started in %s\n", __func__);
171  vout.detailed(m_vl, "asolver_SAP_QWS: solve, started in %s\n", __func__);
172  using real_t = typename AFIELD::real_t;
173 
174 #ifdef USE_QWSLIB
175  m_fopr->convert(x, b);
176 
177 #ifdef DEBUG
178  {
179  double x2 = x.norm2();
180  vout.general("hoge: %s: after convert, |in|^2 = %15.7e\n", class_name.c_str(), x2);
181  }
182 #endif
183  int nsap = m_Niter;
184  int nm = m_nm;
185  m_fopr->mult_clv_inv(x);
186 #pragma omp barrier
187 #ifdef DEBUG
188  {
189  double x2 = x.norm2();
190  vout.general("hoge: %s: after mult_clv_inv, |in|^2 = %15.7e\n", class_name.c_str(), x2);
191  }
192 #endif
193 
194 #pragma omp barrier
195  prec_s_noprl_((scs_t *)m_x.ptr(0), (scs_t *)x.ptr(0), &nsap, &nm, m_s, m_q);
196 
197  diff = -1.0;
198 #ifdef DEBUG
199  {
200  double x2 = m_x.norm2();
201  m_fopr->convert(m_b, b);
202 
203  m_fopr->mult(m_r, m_x);
204  axpy(m_r, (real_t)-1.0, m_b);
205  double r2 = m_r.norm2();
206  diff = r2;
207  vout.general("hoge: %s: before reverse, |res|^2 = %15.7e\n", class_name.c_str(), r2);
208  vout.general("hoge: %s: |out|^2 = %15.7e\n", class_name.c_str(), x2);
209  }
210 #endif
211 
212  m_fopr->reverse(x, m_x);
213 
214 #ifdef DEBUG
215  {
216  double x2 = x.norm2();
217  vout.general("hoge: %s: after reverse, |out|^2 = %15.7e\n", class_name.c_str(), x2);
218  }
219 #endif
220 
221 
222  m_Nconv = m_Niter;
223 
224 #pragma omp barrier
225 #else
226  vout.crucial(m_vl, "%s: in solve(), USE_QWSLIB is not defined\n", class_name.c_str());
227 #endif
228 }
229 
230 
231 //====================================================================
232 template<typename AFIELD>
233 double ASolver_SAP_QWS<AFIELD>::flop_count()
234 {
235  return m_flop;
236 }
237 
238 
239 template<>
240 const std::string ASolver_SAP_QWS<AField<float, QXS> >::class_name
241  = "ASolver_SAP_QWS<AField<float,QXS> >";
242 
243 
244 template class ASolver_SAP_QWS<AField<float, QXS> >;
245 
246 #endif // USE_QWSLIB
247 //============================================================END=====
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
afield-inc.h
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
ASolver_SAP_QWS::flop_count
double flop_count()
returns the floating point operation count.
Parameters
Class for parameters.
Definition: parameters.h:46
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
asolver_SAP_QWS.h
SAP solver (qws version)
Field::real_t
double real_t
Definition: field.h:51
Field::norm2
double norm2() const
Definition: field.cpp:113
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
ASolver_SAP_QWS::init
void init(void)
ASolver_SAP_QWS::tidyup
void tidyup(void)
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
ASolver_SAP_QWS
Definition: asolver_SAP_QWS.h:33
ASolver_SAP_QWS::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
ASolver_SAP_QWS::solve
void solve(AFIELD &xq, const AFIELD &b, int &nconv, real_t &diff)
solver main.
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Field
Container of Field-type object.
Definition: field.h:46
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512