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 77
