lapack/experimental/xlamch.cpp

#include <limits>

#include "../../Eigen/Core"

extern "C" {
float slamch_(char* cmach);
double dlamch_(char* cmach);
float slamc3_(char* cmach);
double dlamc3_(char* cmach);
}

// Determines machine parameters.
template <typename T>
T xlamch(const char c) {
  if (c == 'E' || c == 'e') {
    // eps = relative machine precision. Looking at
    // https://gcc.gnu.org/onlinedocs/gfortran/EPSILON.html observe that it is
    // the smallest number E such that 1 + E > 1. It is NOT the same as
    // numeric_limits::epsilon, which is the difference between 1 and the
    // floating point number that succeeds it.
    return std::numeric_limits<T>::epsilon() / std::numeric_limits<T>::radix;
  }

  if (c == 'S' || c == 's') {
    // sfmin = safe minimum, such that 1 / sfmin does not overflow
    EIGEN_CONSTEXPR T tiny = std::numeric_limits<T>::min();
    EIGEN_CONSTEXPR T small = T(1) / std::numeric_limits<T>::max();
    EIGEN_CONSTEXPR T small_scaled =
        small * (T(1) + std::numeric_limits<T>::epsilon());
    return (small >= tiny) ? small_scaled : tiny;
  }

  if (c == 'B' || c == 'b') {
    // base =  base of the machine.
    return std::numeric_limits<T>::radix;
  }

  if (c == 'P' || c == 'p') {
    // prec = eps * base
    return std::numeric_limits<T>::epsilon();
  }

  if (c == 'N' || c == 'n') {
    // t = number of (base) digits in the mantissa
    return std::numeric_limits<T>::digits;
  }

  if (c == 'R' || c == 'r') {
    // rnd = 1.0 when rounding occurs in addition.
    return 1;
  }

  if (c == 'M' || c == 'm') {
    // emin = minimum exponent before (gradual) underflow
    return std::numeric_limits<T>::min_exponent;
  }

  if (c == 'U' || c == 'u') {
    // rmin =  underflow threshold - base**(emin-1)
    return std::numeric_limits<T>::min();
  }

  if (c == 'L' || c == 'l') {
    // emax = largest exponent before overflow
    return std::numeric_limits<T>::max_exponent;
  }
  if (c == 'O' || c == 'o') {
    // rmax  = overflow threshold  - (base**emax)*(1-eps)
    return std::numeric_limits<T>::max();
  }

  return 0;
}

float slamch_(char* cmach) { return xlamch<float>(*cmach); }
double dlamch_(char* cmach) { return xlamch<double>(*cmach); }

// xLAMC3  is intended to force  A  and  B  to be stored prior to doing
// the addition of  A  and  B ,  for use in situations where optimizers
// might hold one of these in a register.
float slamc3_(float* a, float* b) {
  volatile float retval = *a + *b;
  return retval;
}

double dlamc3_(double* a, double* b) {
  volatile double retval = *a + *b;
  return retval;
}