1 #include "cminpack.h"
2 #include <math.h>
3 #include "cminpackP.h"
4
5 #define log10e 0.43429448190325182765
6 #define factor 100.
7
8 /* Table of constant values */
9
10 __cminpack_attr__
__cminpack_func__(chkder)11 void __cminpack_func__(chkder)(int m, int n, const real *x,
12 real *fvec, real *fjac, int ldfjac, real *xp,
13 real *fvecp, int mode, real *err)
14 {
15 /* Local variables */
16 int i, j;
17 real eps, epsf, temp, epsmch;
18 real epslog;
19
20 /* ********** */
21
22 /* subroutine chkder */
23
24 /* this subroutine checks the gradients of m nonlinear functions */
25 /* in n variables, evaluated at a point x, for consistency with */
26 /* the functions themselves. the user must call chkder twice, */
27 /* first with mode = 1 and then with mode = 2. */
28
29 /* mode = 1. on input, x must contain the point of evaluation. */
30 /* on output, xp is set to a neighboring point. */
31
32 /* mode = 2. on input, fvec must contain the functions and the */
33 /* rows of fjac must contain the gradients */
34 /* of the respective functions each evaluated */
35 /* at x, and fvecp must contain the functions */
36 /* evaluated at xp. */
37 /* on output, err contains measures of correctness of */
38 /* the respective gradients. */
39
40 /* the subroutine does not perform reliably if cancellation or */
41 /* rounding errors cause a severe loss of significance in the */
42 /* evaluation of a function. therefore, none of the components */
43 /* of x should be unusually small (in particular, zero) or any */
44 /* other value which may cause loss of significance. */
45
46 /* the subroutine statement is */
47
48 /* subroutine chkder(m,n,x,fvec,fjac,ldfjac,xp,fvecp,mode,err) */
49
50 /* where */
51
52 /* m is a positive integer input variable set to the number */
53 /* of functions. */
54
55 /* n is a positive integer input variable set to the number */
56 /* of variables. */
57
58 /* x is an input array of length n. */
59
60 /* fvec is an array of length m. on input when mode = 2, */
61 /* fvec must contain the functions evaluated at x. */
62
63 /* fjac is an m by n array. on input when mode = 2, */
64 /* the rows of fjac must contain the gradients of */
65 /* the respective functions evaluated at x. */
66
67 /* ldfjac is a positive integer input parameter not less than m */
68 /* which specifies the leading dimension of the array fjac. */
69
70 /* xp is an array of length n. on output when mode = 1, */
71 /* xp is set to a neighboring point of x. */
72
73 /* fvecp is an array of length m. on input when mode = 2, */
74 /* fvecp must contain the functions evaluated at xp. */
75
76 /* mode is an integer input variable set to 1 on the first call */
77 /* and 2 on the second. other values of mode are equivalent */
78 /* to mode = 1. */
79
80 /* err is an array of length m. on output when mode = 2, */
81 /* err contains measures of correctness of the respective */
82 /* gradients. if there is no severe loss of significance, */
83 /* then if err(i) is 1.0 the i-th gradient is correct, */
84 /* while if err(i) is 0.0 the i-th gradient is incorrect. */
85 /* for values of err between 0.0 and 1.0, the categorization */
86 /* is less certain. in general, a value of err(i) greater */
87 /* than 0.5 indicates that the i-th gradient is probably */
88 /* correct, while a value of err(i) less than 0.5 indicates */
89 /* that the i-th gradient is probably incorrect. */
90
91 /* subprograms called */
92
93 /* minpack supplied ... dpmpar */
94
95 /* fortran supplied ... dabs,dlog10,dsqrt */
96
97 /* argonne national laboratory. minpack project. march 1980. */
98 /* burton s. garbow, kenneth e. hillstrom, jorge j. more */
99
100 /* ********** */
101
102 /* epsmch is the machine precision. */
103
104 epsmch = __cminpack_func__(dpmpar)(1);
105
106 eps = sqrt(epsmch);
107
108 if (mode != 2) {
109
110 /* mode = 1. */
111
112 for (j = 0; j < n; ++j) {
113 temp = eps * fabs(x[j]);
114 if (temp == 0.) {
115 temp = eps;
116 }
117 xp[j] = x[j] + temp;
118 }
119 return;
120 }
121
122 /* mode = 2. */
123
124 epsf = factor * epsmch;
125 epslog = log10e * log(eps);
126 for (i = 0; i < m; ++i) {
127 err[i] = 0.;
128 }
129 for (j = 0; j < n; ++j) {
130 temp = fabs(x[j]);
131 if (temp == 0.) {
132 temp = 1.;
133 }
134 for (i = 0; i < m; ++i) {
135 err[i] += temp * fjac[i + j * ldfjac];
136 }
137 }
138 for (i = 0; i < m; ++i) {
139 temp = 1.;
140 if (fvec[i] != 0. && fvecp[i] != 0. &&
141 fabs(fvecp[i] - fvec[i]) >= epsf * fabs(fvec[i]))
142 {
143 temp = eps * fabs((fvecp[i] - fvec[i]) / eps - err[i])
144 / (fabs(fvec[i]) +
145 fabs(fvecp[i]));
146 }
147 err[i] = 1.;
148 if (temp > epsmch && temp < eps) {
149 err[i] = (log10e * log(temp) - epslog) / epslog;
150 }
151 if (temp >= eps) {
152 err[i] = 0.;
153 }
154 }
155
156 /* last card of subroutine chkder. */
157
158 } /* chkder_ */
159
160