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