25 bool init = Solver::Factory::Register(
"BiCGStab_IDS_L_Cmplx", create_object);
35 const string str_vlevel = params.
get_string(
"verbose_level");
46 err += params.
fetch_int(
"maximum_number_of_iteration", Niter);
47 err += params.
fetch_int(
"maximum_number_of_restart", Nrestart);
48 err += params.
fetch_double(
"convergence_criterion_squared", Stop_cond);
49 err += params.
fetch_int(
"number_of_orthonormal_vectors", N_L);
50 err += params.
fetch_double(
"tolerance_for_DynamicSelection_of_L", Tol_L);
117 int& Nconv,
double& diff)
119 double bnorm2 = b.
norm2();
120 int bsize = b.
size();
126 bool is_converged =
false;
139 for (
int i_restart = 0; i_restart <
m_Nrestart; i_restart++) {
140 for (
int iter = 0; iter <
m_Niter; iter++) {
169 if (is_converged)
break;
175 vout.
crucial(
m_vl,
" iter(final): %8d %22.15e\n", Nconv2, diff2 / bnorm2);
185 diff = sqrt(diff2 / bnorm2);
212 for (
int i = 0; i <
m_N_L + 1; ++i) {
213 m_u[i].reset(Nin, Nvol, Nex);
214 m_r[i].reset(Nin, Nvol, Nex);
228 for (
int i = 0; i <
m_N_L + 1; ++i) {
262 dcomplex c_Rayleigh_prev = cmplx(0.0, 0.0);
264 bool is_converged_L =
false;
268 if (!is_converged_L) {
272 dcomplex beta = alpha_prev2 * (rho / rho_prev2);
276 for (
int i = 0; i < j + 1; ++i) {
284 alpha_prev2 = rho_prev2 / conj(gamma);
286 for (
int i = 0; i < j + 1; ++i) {
298 double r_tmp =
m_r[j].norm2();
299 dcomplex c_Rayleigh =
dotc(
m_r[j],
m_r[j + 1]) / r_tmp;
301 dcomplex c_E = (c_Rayleigh - c_Rayleigh_prev) / c_Rayleigh;
304 c_Rayleigh_prev = c_Rayleigh;
314 is_converged_L =
true;
320 std::vector<double> sigma(
m_N_L + 1);
321 std::vector<dcomplex> gamma_prime(
m_N_L + 1);
324 std::vector<dcomplex> tau(
m_N_L * (
m_N_L + 1));
327 for (
int j = 1; j < N_L_tmp + 1; ++j) {
328 for (
int i = 1; i < j; ++i) {
332 tau[ij] = conj(r_ji) / sigma[i];
336 sigma[j] =
m_r[j].norm2();
339 gamma_prime[j] = conj(r_0j) / sigma[j];
343 std::vector<dcomplex> gamma(
m_N_L + 1);
346 gamma[N_L_tmp] = gamma_prime[N_L_tmp];
348 for (
int j = N_L_tmp - 1; j > 0; --j) {
349 c_tmp = cmplx(0.0, 0.0);
351 for (
int i = j + 1; i < N_L_tmp + 1; ++i) {
353 c_tmp += tau[ji] * gamma[i];
356 gamma[j] = gamma_prime[j] - c_tmp;
361 std::vector<dcomplex> gamma_double_prime(
m_N_L);
363 for (
int j = 1; j < N_L_tmp; ++j) {
364 c_tmp = cmplx(0.0, 0.0);
366 for (
int i = j + 1; i < N_L_tmp; ++i) {
368 c_tmp += tau[ji] * gamma[i + 1];
371 gamma_double_prime[j] = gamma[j + 1] + c_tmp;
376 axpy(
m_r[0], -gamma_prime[N_L_tmp],
m_r[N_L_tmp]);
379 for (
int j = 1; j < N_L_tmp; ++j) {
402 m_N_L_prev = N_L_tmp;
425 if (flop_fopr < eps) {
430 double flop_axpy =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
431 double flop_dotc =
static_cast<double>(Nin * Nex * 4) * (Nvol * NPE);
432 double flop_norm =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
436 double flop_init = flop_fopr + flop_axpy + flop_norm;
437 double flop_step_BiCG_part = 2 * N_L_prev_total * flop_fopr
438 + 3 * N_L_prev_total * flop_dotc
440 + N_L_prev_total * flop_norm;
444 double flop_step = flop_step_BiCG_part + flop_step_L_part;
445 double flop_true_residual = flop_fopr + flop_axpy + flop_norm;
447 double flop = flop_norm + flop_init + flop_step + flop_true_residual
int index_ij(int i, int j)
void detailed(const char *format,...)
static double epsilon_criterion()
void set_parameters_DS_L(const int N_L, const double Tol_L)
void general(const char *format,...)
Container of Field-type object.
int fetch_double(const string &key, double &value) const
void copy(Field &y, const Field &x)
copy(y, x): y = x
BiCGStab(IDS_L) algorithm.
void solve_init(const Field &, double &)
int square_non_zero(const double v)
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
static const std::string class_name
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)
void solve_step(double &)
void reset_field(const Field &)
int non_negative(const int v)
void set_parameters(const Parameters ¶ms)
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 solve(Field &solution, const Field &source, int &Nconv, double &diff)
Bridge::VerboseLevel m_vl
static VerboseLevel set_verbose_level(const std::string &str)