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