fdjac1.h
1 namespace Eigen {
2 
3 namespace internal {
4 
5 template<typename FunctorType, typename Scalar>
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