Magic numbers in Givens rotation regularization

[pelesh]

Givens rotation in (F)GMRES solvers is implemented as

        real_type gam = std::sqrt(Hii * Hii + Hii1 * Hii1);

        if(std::abs(gam - ZERO) <= EPSILON) {
          gam = EPSMAC;
        }

        /* next Given's rotation */
        h_c_[i] = Hii / gam;

where EPSILON and EPSMAC are defined as

  constexpr double EPSILON = 1.0e-18;
  constexpr double EPSMAC  = 1.0e-16;

It would be good to document these numbers to better understand why these specific values have been chosen.

Also, these values will probably not be appropriate if real_type is not chosen to be double floating precision. Perhaps a better idea would be to define EPSILON and EPSMAC in terms of std::limits::epsilon().