28 bool init = Solver::Factory::Register(
"GMRES_m_Cmplx", create_object);
38 const string str_vlevel = params.
get_string(
"verbose_level");
48 err += params.
fetch_int(
"maximum_number_of_iteration", Niter);
49 err += params.
fetch_int(
"maximum_number_of_restart", Nrestart);
50 err += params.
fetch_double(
"convergence_criterion_squared", Stop_cond);
51 err += params.
fetch_int(
"number_of_orthonormal_vectors", N_M);
114 int& Nconv,
double& diff)
116 double bnorm2 = b.
norm2();
117 int bsize = b.
size();
123 bool is_converged =
false;
136 for (
int i_restart = 0; i_restart <
m_Nrestart; i_restart++) {
137 for (
int iter = 0; iter <
m_Niter; iter++) {
169 vout.
crucial(
m_vl,
" iter(final): %8d %22.15e\n", Nconv2, diff2 / bnorm2);
179 diff = sqrt(diff2 / bnorm2);
205 for (
int i = 0; i <
m_N_M + 1; ++i) {
206 m_v[i].reset(Nin, Nvol, Nex);
221 for (
int i = 0; i <
m_N_M + 1; ++i) {
252 for (
int j = 0; j <
m_N_M; ++j) {
255 for (
int i = 0; i < j + 1; ++i) {
263 for (
int i = 0; i < j + 1; ++i) {
268 double v_norm2 =
m_v[j + 1].norm2();
271 h[j1j] = cmplx(sqrt(v_norm2), 0.0);
273 scal(
m_v[j + 1], 1.0 / sqrt(v_norm2));
282 for (
int i = 0; i <
m_N_M; ++i) {
295 valarray<dcomplex>& h)
299 valarray<dcomplex> g(
m_N_M + 1);
304 for (
int i = 0; i <
m_N_M; ++i) {
306 double h_1_r = abs(h[ii]);
309 double h_2_r = abs(h[i1i]);
311 double denomi = sqrt(h_1_r * h_1_r + h_2_r * h_2_r);
313 dcomplex cs = h[ii] / denomi;
314 dcomplex sn = h[i1i] / denomi;
316 for (
int j = i; j <
m_N_M; ++j) {
320 dcomplex const_1_c = conj(cs) * h[ij] + sn * h[i1j];
321 dcomplex const_2_c = -sn * h[ij] + cs * h[i1j];
327 dcomplex const_1_c = conj(cs) * g[i] + sn * g[i + 1];
328 dcomplex const_2_c = -sn * g[i] + cs * g[i + 1];
331 g[i + 1] = const_2_c;
335 for (
int i = m_N_M - 1; i > -1; --i) {
336 for (
int j = i + 1; j <
m_N_M; ++j) {
338 g[i] -= h[ij] * y[j];
361 if (flop_fopr < eps) {
366 double flop_axpy =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
367 double flop_dotc =
static_cast<double>(Nin * Nex * 4) * (Nvol * NPE);
368 double flop_norm =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
369 double flop_scal =
static_cast<double>(Nin * Nex * 2) * (Nvol * NPE);
374 for (
int j = 0; j <
m_N_M; ++j) {
375 for (
int i = 0; i < j + 1; ++i) {
380 double flop_init = flop_fopr + flop_axpy + flop_norm;
381 double flop_step = m_N_M * flop_fopr
382 + N_M_part * flop_dotc
383 + (N_M_part +
m_N_M) * flop_axpy
385 double flop_true_residual = flop_fopr + flop_axpy + flop_norm;
387 double flop = flop_norm + flop_init + flop_step * N_iter + flop_true_residual
void scal(Field &x, const double a)
scal(x, a): x = a * x
void detailed(const char *format,...)
static double epsilon_criterion()
void general(const char *format,...)
static const std::string class_name
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)
GMRES(m) algorithm with complex variables.
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)
int index_ij(int i, int j)
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.
void reset_field(const Field &)
virtual void mult(Field &, const Field &)=0
multiplies fermion operator to a given field (2nd argument)
void solve(Field &solution, const Field &source, int &Nconv, double &diff)
int non_negative(const int v)
void min_J(std::valarray< dcomplex > &y, std::valarray< dcomplex > &h)
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
void set_parameters_GMRES_m(const int N_M)
Base class of fermion operator family.
string get_string(const string &key) const
Bridge::VerboseLevel m_vl
static VerboseLevel set_verbose_level(const std::string &str)
void solve_step(const Field &, double &)
void set_parameters(const Parameters ¶ms)