1 namespace Eigen {
2
3 namespace internal {
4
5 template<typename FunctorType, typename Scalar>
fdjac1(const FunctorType & Functor,Matrix<Scalar,Dynamic,1> & x,Matrix<Scalar,Dynamic,1> & fvec,Matrix<Scalar,Dynamic,Dynamic> & fjac,DenseIndex ml,DenseIndex mu,Scalar epsfcn)6 DenseIndex fdjac1(
7 const FunctorType &Functor,
8 Matrix< Scalar, Dynamic, 1 > &x,
9 Matrix< Scalar, Dynamic, 1 > &fvec,
10 Matrix< Scalar, Dynamic, Dynamic > &fjac,
11 DenseIndex ml, DenseIndex mu,
12 Scalar epsfcn)
13 {
14 typedef DenseIndex Index;
15
16 /* Local variables */
17 Scalar h;
18 Index j, k;
19 Scalar eps, temp;
20 Index msum;
21 int iflag;
22 Index start, length;
23
24 /* Function Body */
25 const Scalar epsmch = NumTraits<Scalar>::epsilon();
26 const Index n = x.size();
27 assert(fvec.size()==n);
28 Matrix< Scalar, Dynamic, 1 > wa1(n);
29 Matrix< Scalar, Dynamic, 1 > wa2(n);
30
31 eps = sqrt((std::max)(epsfcn,epsmch));
32 msum = ml + mu + 1;
33 if (msum >= n) {
34 /* computation of dense approximate jacobian. */
35 for (j = 0; j < n; ++j) {
36 temp = x[j];
37 h = eps * abs(temp);
38 if (h == 0.)
39 h = eps;
40 x[j] = temp + h;
41 iflag = Functor(x, wa1);
42 if (iflag < 0)
43 return iflag;
44 x[j] = temp;
45 fjac.col(j) = (wa1-fvec)/h;
46 }
47
48 }else {
49 /* computation of banded approximate jacobian. */
50 for (k = 0; k < msum; ++k) {
51 for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
52 wa2[j] = x[j];
53 h = eps * abs(wa2[j]);
54 if (h == 0.) h = eps;
55 x[j] = wa2[j] + h;
56 }
57 iflag = Functor(x, wa1);
58 if (iflag < 0)
59 return iflag;
60 for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
61 x[j] = wa2[j];
62 h = eps * abs(wa2[j]);
63 if (h == 0.) h = eps;
64 fjac.col(j).setZero();
65 start = std::max<Index>(0,j-mu);
66 length = (std::min)(n-1, j+ml) - start + 1;
67 fjac.col(j).segment(start, length) = ( wa1.segment(start, length)-fvec.segment(start, length))/h;
68 }
69 }
70 }
71 return 0;
72 }
73
74 } // end namespace internal
75
76 } // end namespace Eigen
77