1 files changed, 157 insertions, 0 deletions
diff --git a/ast/cminpack/enorm.c b/ast/cminpack/enorm.c
new file mode 100644
index 0000000..ad10824
--- /dev/null
+++ b/ast/cminpack/enorm.c
@@ -0,0 +1,157 @@
+#include "cminpack.h"
+#include <math.h>
+#include "cminpackP.h"
+
+/*
+  About the values for rdwarf and rgiant.
+
+  The original values, both in signe-precision FORTRAN source code and in double-precision code were:
+#define rdwarf 3.834e-20
+#define rgiant 1.304e19
+  See for example:
+    http://www.netlib.org/slatec/src/denorm.f
+    http://www.netlib.org/slatec/src/enorm.f
+  However, rdwarf is smaller than sqrt(FLT_MIN) = 1.0842021724855044e-19, so that rdwarf**2 will
+  underflow. This contradicts the constraints expressed in the comments below.
+
+  We changed these constants to be sqrt(dpmpar(2))*0.9 and sqrt(dpmpar(3))*0.9, as proposed by the
+  implementation found in MPFIT http://cow.physics.wisc.edu/~craigm/idl/fitting.html
+*/
+
+#define double_dwarf (1.4916681462400413e-154*0.9)
+#define double_giant (1.3407807929942596e+154*0.9)
+#define float_dwarf (1.0842021724855044e-19f*0.9f)
+#define float_giant (1.8446743523953730e+19f*0.9f)
+#define half_dwarf (2.4414062505039999e-4f*0.9f)
+#define half_giant (255.93749236874225497222f*0.9f)
+
+#define dwarf(type) _dwarf(type)
+#define _dwarf(type) type ## _dwarf
+#define giant(type) _giant(type)
+#define _giant(type) type ## _giant
+
+#define rdwarf dwarf(real)
+#define rgiant giant(real)
+
+__cminpack_attr__
+real __cminpack_func__(enorm)(int n, const real *x)
+{
+#ifdef USE_CBLAS
+    return cblas_dnrm2(n, x, 1);
+#else /* !USE_CBLAS */
+    /* System generated locals */
+    real ret_val, d1;
+
+    /* Local variables */
+    int i;
+    real s1, s2, s3, xabs, x1max, x3max, agiant, floatn;
+
+/*     ********** */
+
+/*     function enorm */
+
+/*     given an n-vector x, this function calculates the */
+/*     euclidean norm of x. */
+
+/*     the euclidean norm is computed by accumulating the sum of */
+/*     squares in three different sums. the sums of squares for the */
+/*     small and large components are scaled so that no overflows */
+/*     occur. non-destructive underflows are permitted. underflows */
+/*     and overflows do not occur in the computation of the unscaled */
+/*     sum of squares for the intermediate components. */
+/*     the definitions of small, intermediate and large components */
+/*     depend on two constants, rdwarf and rgiant. the main */
+/*     restrictions on these constants are that rdwarf**2 not */
+/*     underflow and rgiant**2 not overflow. the constants */
+/*     given here are suitable for every known computer. */
+
+/*     the function statement is */
+
+/*       double precision function enorm(n,x) */
+
+/*     where */
+
+/*       n is a positive integer input variable. */
+
+/*       x is an input array of length n. */
+
+/*     subprograms called */
+
+/*       fortran-supplied ... dabs,dsqrt */
+
+/*     argonne national laboratory. minpack project. march 1980. */
+/*     burton s. garbow, kenneth e. hillstrom, jorge j. more */
+
+/*     ********** */
+
+    s1 = 0.;
+    s2 = 0.;
+    s3 = 0.;
+    x1max = 0.;
+    x3max = 0.;
+    floatn = (real) (n);
+    agiant = rgiant / floatn;
+    for (i = 0; i < n; ++i) {
+	xabs = fabs(x[i]);
+	if (xabs <= rdwarf || xabs >= agiant) {
+            if (xabs > rdwarf) {
+
+/*              sum for large components. */
+
+                if (xabs > x1max) {
+                    /* Computing 2nd power */
+                    d1 = x1max / xabs;
+                    s1 = 1. + s1 * (d1 * d1);
+                    x1max = xabs;
+                } else {
+                    /* Computing 2nd power */
+                    d1 = xabs / x1max;
+                    s1 += d1 * d1;
+                }
+            } else {
+
+/*              sum for small components. */
+
+                if (xabs > x3max) {
+                    /* Computing 2nd power */
+                    d1 = x3max / xabs;
+                    s3 = 1. + s3 * (d1 * d1);
+                    x3max = xabs;
+                } else {
+                    if (xabs != 0.) {
+                        /* Computing 2nd power */
+                        d1 = xabs / x3max;
+                        s3 += d1 * d1;
+                    }
+                }
+            }
+	} else {
+
+/*           sum for intermediate components. */
+
+            /* Computing 2nd power */
+            s2 += xabs * xabs;
+        }
+    }
+
+/*     calculation of norm. */
+
+    if (s1 != 0.) {
+        ret_val = x1max * sqrt(s1 + (s2 / x1max) / x1max);
+    } else {
+        if (s2 != 0.) {
+            if (s2 >= x3max) {
+                ret_val = sqrt(s2 * (1. + (x3max / s2) * (x3max * s3)));
+            } else {
+                ret_val = sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+            }
+        } else {
+            ret_val = x3max * sqrt(s3);
+        }
+    }
+    return ret_val;
+
+/*     last card of function enorm. */
+#endif /* !USE_CBLAS */
+} /* enorm_ */
+