diff options
Diffstat (limited to 'ast/cminpack/enorm.c')
| -rw-r--r-- | ast/cminpack/enorm.c | 157 |
1 files changed, 0 insertions, 157 deletions
diff --git a/ast/cminpack/enorm.c b/ast/cminpack/enorm.c deleted file mode 100644 index ad10824..0000000 --- a/ast/cminpack/enorm.c +++ /dev/null @@ -1,157 +0,0 @@ -#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_ */ - |