1. C++ implementations of ilaver to report the Eigen version number. Though you could make the argument that it should report the LAPACK test suite number that this passes.
2. xLAMCH implementation using <limits>. I can simplify the code for sfmin further because I know which branch of the if statement it takes.
3. Because xLAMCH and xLAMC2 live in the same fortran file, our overlay technique for testing LAPACK runs into problems because of linker issues. So we also split that function into its own file.
PiperOrigin-RevId: 381130439
diff --git a/lapack/experimental/ilaver.cpp b/lapack/experimental/ilaver.cpp
new file mode 100644
index 0000000..ddab5dd
--- /dev/null
+++ b/lapack/experimental/ilaver.cpp
@@ -0,0 +1,11 @@
+#include "../../Eigen/Core"
+
+extern "C" {
+void ilaver_(int* major, int* minor, int* patch);
+}
+
+void ilaver_(int* major, int* minor, int* patch) {
+ *major = EIGEN_WORLD_VERSION;
+ *minor = EIGEN_MAJOR_VERSION;
+ *patch = EIGEN_MINOR_VERSION;
+}
diff --git a/lapack/experimental/xlamch.cpp b/lapack/experimental/xlamch.cpp
new file mode 100644
index 0000000..df1527b
--- /dev/null
+++ b/lapack/experimental/xlamch.cpp
@@ -0,0 +1,89 @@
+#include <limits>
+
+#include "../../Eigen/src/Core/util/Macros.h"
+
+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;
+}
