25 bool init = Solver::Factory::Register(
"BiCGStab_L_Cmplx", create_object);
35 const string str_vlevel = params.
get_string(
"verbose_level");
45 err += params.
fetch_int(
"maximum_number_of_iteration", Niter);
46 err += params.
fetch_int(
"maximum_number_of_restart", Nrestart);
47 err += params.
fetch_double(
"convergence_criterion_squared", Stop_cond);
48 err += params.
fetch_int(
"number_of_orthonormal_vectors", N_L);
111 int& Nconv,
double& diff)
113 double bnorm2 = b.
norm2();
114 int bsize = b.
size();
120 bool is_converged =
false;
133 for (
int i_restart = 0; i_restart <
m_Nrestart; i_restart++) {
134 for (
int iter = 0; iter <
m_Niter; iter++) {
166 vout.
crucial(
m_vl,
" iter(final): %8d %22.15e\n", Nconv2, diff2 / bnorm2);
176 diff = sqrt(diff2 / bnorm2);
203 for (
int i = 0; i <
m_N_L + 1; ++i) {
204 m_u[i].reset(Nin, Nvol, Nex);
205 m_r[i].reset(Nin, Nvol, Nex);
219 for (
int i = 0; i <
m_N_L + 1; ++i) {
250 for (
int j = 0; j <
m_N_L; ++j) {
254 dcomplex beta = alpha_prev2 * (rho / rho_prev2);
258 for (
int i = 0; i < j + 1; ++i) {
265 alpha_prev2 = rho_prev2 / conj(gamma);
267 for (
int i = 0; i < j + 1; ++i) {
277 std::vector<double> sigma(m_N_L + 1);
278 std::vector<dcomplex> gamma_prime(m_N_L + 1);
281 std::vector<dcomplex> tau(m_N_L * (m_N_L + 1));
284 for (
int j = 1; j < m_N_L + 1; ++j) {
285 for (
int i = 1; i < j; ++i) {
289 tau[ij] = conj(r_ji) / sigma[i];
293 sigma[j] =
m_r[j].norm2();
296 gamma_prime[j] = conj(r_0j) / sigma[j];
300 std::vector<dcomplex> gamma(m_N_L + 1);
305 for (
int j = m_N_L - 1; j > 0; --j) {
306 c_tmp = cmplx(0.0, 0.0);
308 for (
int i = j + 1; i < m_N_L + 1; ++i) {
310 c_tmp += tau[ji] * gamma[i];
313 gamma[j] = gamma_prime[j] - c_tmp;
318 std::vector<dcomplex> gamma_double_prime(m_N_L);
320 for (
int j = 1; j <
m_N_L; ++j) {
321 c_tmp = cmplx(0.0, 0.0);
323 for (
int i = j + 1; i <
m_N_L; ++i) {
325 c_tmp += tau[ji] * gamma[i + 1];
328 gamma_double_prime[j] = gamma[j + 1] + c_tmp;
333 axpy(
m_r[0], -gamma_prime[m_N_L],
m_r[m_N_L]);
336 for (
int j = 1; j <
m_N_L; ++j) {
370 if (flop_fopr < eps) {
375 double flop_axpy =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
376 double flop_dotc =
static_cast<double>(Nin * Nex * 4) * (Nvol * NPE);
377 double flop_norm =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
382 for (
int j = 0; j <
m_N_L; ++j) {
383 for (
int i = 0; i < j + 1; ++i) {
388 double flop_init = flop_fopr + flop_axpy + flop_norm;
389 double flop_step_BiCG_part = 2 * m_N_L * flop_fopr
390 + 2 * m_N_L * flop_dotc
391 + (m_N_L + 2 * N_L_part) * flop_axpy;
392 double flop_step_L_part = (N_L_part +
m_N_L) * flop_dotc
393 + (N_L_part + 3 * m_N_L) * flop_axpy
394 + (m_N_L + 1) * flop_norm;
395 double flop_step = flop_step_BiCG_part + flop_step_L_part;
396 double flop_true_residual = flop_fopr + flop_axpy + flop_norm;
398 double flop = flop_norm + flop_init + flop_step * N_iter + flop_true_residual
static const std::string class_name
void detailed(const char *format,...)
int index_ij(int i, int j)
static double epsilon_criterion()
void general(const char *format,...)
void set_parameters_L(const int N_L)
void solve_step(double &)
Container of Field-type object.
int fetch_double(const string &key, double &value) const
void solve_init(const Field &, double &)
void copy(Field &y, const Field &x)
copy(y, x): y = x
int square_non_zero(const double v)
void set_parameters(const Parameters ¶ms)
dcomplex dotc(const Field &y, const Field &x)
int fetch_int(const string &key, int &value) const
virtual double flop_count()
returns the flops per site.
void reset(const int Nin, const int Nvol, const int Nex, const element_type cmpl=COMPLEX)
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
void paranoiac(const char *format,...)
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
void crucial(const char *format,...)
Base class for linear solver class family.
virtual void mult(Field &, const Field &)=0
multiplies fermion operator to a given field (2nd argument)
int non_negative(const int v)
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Base class of fermion operator family.
string get_string(const string &key) const
void reset_field(const Field &)
Bridge::VerboseLevel m_vl
static VerboseLevel set_verbose_level(const std::string &str)
void solve(Field &solution, const Field &source, int &Nconv, double &diff)