lapack/experimental/xlamch.cpp

#include <limits>

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
    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
    T sfmin = std::numeric_limits<T>::min();
    const T small = T(1) / std::numeric_limits<T>::max();
    if (small >= sfmin) {
      sfmin = small * (T(1) + std::numeric_limits<T>::epsilon() /
                                  std::numeric_limits<T>::radix);
    }
    return sfmin;
  }

  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;
}