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
|
#include "cminpack.h"
#include <math.h>
#include "cminpackP.h"
__cminpack_attr__
int __cminpack_func__(fdjac2)(__cminpack_decl_fcn_mn__ void *p, int m, int n, real *x,
const real *fvec, real *fjac, int ldfjac,
real epsfcn, real *wa)
{
/* Local variables */
real h;
int i, j;
real eps, temp, epsmch;
int iflag;
/* ********** */
/* subroutine fdjac2 */
/* this subroutine computes a forward-difference approximation */
/* to the m by n jacobian matrix associated with a specified */
/* problem of m functions in n variables. */
/* the subroutine statement is */
/* subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa) */
/* where */
/* fcn is the name of the user-supplied subroutine which */
/* calculates the functions. fcn must be declared */
/* in an external statement in the user calling */
/* program, and should be written as follows. */
/* subroutine fcn(m,n,x,fvec,iflag) */
/* integer m,n,iflag */
/* double precision x(n),fvec(m) */
/* ---------- */
/* calculate the functions at x and */
/* return this vector in fvec. */
/* ---------- */
/* return */
/* end */
/* the value of iflag should not be changed by fcn unless */
/* the user wants to terminate execution of fdjac2. */
/* in this case set iflag to a negative integer. */
/* m is a positive integer input variable set to the number */
/* of functions. */
/* n is a positive integer input variable set to the number */
/* of variables. n must not exceed m. */
/* x is an input array of length n. */
/* fvec is an input array of length m which must contain the */
/* functions evaluated at x. */
/* fjac is an output m by n array which contains the */
/* approximation to the jacobian matrix evaluated at x. */
/* ldfjac is a positive integer input variable not less than m */
/* which specifies the leading dimension of the array fjac. */
/* iflag is an integer variable which can be used to terminate */
/* the execution of fdjac2. see description of fcn. */
/* epsfcn is an input variable used in determining a suitable */
/* step length for the forward-difference approximation. this */
/* approximation assumes that the relative errors in the */
/* functions are of the order of epsfcn. if epsfcn is less */
/* than the machine precision, it is assumed that the relative */
/* errors in the functions are of the order of the machine */
/* precision. */
/* wa is a work array of length m. */
/* subprograms called */
/* user-supplied ...... fcn */
/* minpack-supplied ... dpmpar */
/* fortran-supplied ... dabs,dmax1,dsqrt */
/* argonne national laboratory. minpack project. march 1980. */
/* burton s. garbow, kenneth e. hillstrom, jorge j. more */
/* ********** */
/* epsmch is the machine precision. */
epsmch = __cminpack_func__(dpmpar)(1);
eps = sqrt((max(epsfcn,epsmch)));
for (j = 0; j < n; ++j) {
temp = x[j];
h = eps * fabs(temp);
if (h == 0.) {
h = eps;
}
x[j] = temp + h;
/* the last parameter of fcn_mn() is set to 2 to differentiate
calls made to compute the function from calls made to compute
the Jacobian (see fcn() in examples/lmfdrv.c, and how njev
is used to compute the number of Jacobian evaluations) */
iflag = fcn_mn(p, m, n, x, wa, 2);
if (iflag < 0) {
return iflag;
}
x[j] = temp;
for (i = 0; i < m; ++i) {
fjac[i + j * ldfjac] = (wa[i] - fvec[i]) / h;
}
}
return 0;
/* last card of subroutine fdjac2. */
} /* fdjac2_ */
|
