1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
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_ */
|
