![]() |
Eigen-unsupported
3.3.3
|
00001 namespace Eigen { 00002 00003 namespace internal { 00004 00005 template<typename FunctorType, typename Scalar> 00006 DenseIndex fdjac1( 00007 const FunctorType &Functor, 00008 Matrix< Scalar, Dynamic, 1 > &x, 00009 Matrix< Scalar, Dynamic, 1 > &fvec, 00010 Matrix< Scalar, Dynamic, Dynamic > &fjac, 00011 DenseIndex ml, DenseIndex mu, 00012 Scalar epsfcn) 00013 { 00014 using std::sqrt; 00015 using std::abs; 00016 00017 typedef DenseIndex Index; 00018 00019 /* Local variables */ 00020 Scalar h; 00021 Index j, k; 00022 Scalar eps, temp; 00023 Index msum; 00024 int iflag; 00025 Index start, length; 00026 00027 /* Function Body */ 00028 const Scalar epsmch = NumTraits<Scalar>::epsilon(); 00029 const Index n = x.size(); 00030 eigen_assert(fvec.size()==n); 00031 Matrix< Scalar, Dynamic, 1 > wa1(n); 00032 Matrix< Scalar, Dynamic, 1 > wa2(n); 00033 00034 eps = sqrt((std::max)(epsfcn,epsmch)); 00035 msum = ml + mu + 1; 00036 if (msum >= n) { 00037 /* computation of dense approximate jacobian. */ 00038 for (j = 0; j < n; ++j) { 00039 temp = x[j]; 00040 h = eps * abs(temp); 00041 if (h == 0.) 00042 h = eps; 00043 x[j] = temp + h; 00044 iflag = Functor(x, wa1); 00045 if (iflag < 0) 00046 return iflag; 00047 x[j] = temp; 00048 fjac.col(j) = (wa1-fvec)/h; 00049 } 00050 00051 }else { 00052 /* computation of banded approximate jacobian. */ 00053 for (k = 0; k < msum; ++k) { 00054 for (j = k; (msum<0) ? (j>n): (j<n); j += msum) { 00055 wa2[j] = x[j]; 00056 h = eps * abs(wa2[j]); 00057 if (h == 0.) h = eps; 00058 x[j] = wa2[j] + h; 00059 } 00060 iflag = Functor(x, wa1); 00061 if (iflag < 0) 00062 return iflag; 00063 for (j = k; (msum<0) ? (j>n): (j<n); j += msum) { 00064 x[j] = wa2[j]; 00065 h = eps * abs(wa2[j]); 00066 if (h == 0.) h = eps; 00067 fjac.col(j).setZero(); 00068 start = std::max<Index>(0,j-mu); 00069 length = (std::min)(n-1, j+ml) - start + 1; 00070 fjac.col(j).segment(start, length) = ( wa1.segment(start, length)-fvec.segment(start, length))/h; 00071 } 00072 } 00073 } 00074 return 0; 00075 } 00076 00077 } // end namespace internal 00078 00079 } // end namespace Eigen