src/NonLinearOptimization/fdjac1.h

0001 namespace Eigen {
0002
0003 namespace internal {
0004
0005 template<typename FunctorType, typename Scalar>
0006 DenseIndex fdjac1(
0007         const FunctorType &Functor,
0008         Matrix< Scalar, Dynamic, 1 >  &x,
0009         Matrix< Scalar, Dynamic, 1 >  &fvec,
0010         Matrix< Scalar, Dynamic, Dynamic > &fjac,
0011         DenseIndex ml, DenseIndex mu,
0012         Scalar epsfcn)
0013 {
0014     using std::sqrt;
0015     using std::abs;
0016
0017     typedef DenseIndex Index;
0018
0019     /* Local variables */
0020     Scalar h;
0021     Index j, k;
0022     Scalar eps, temp;
0023     Index msum;
0024     int iflag;
0025     Index start, length;
0026
0027     /* Function Body */
0028     const Scalar epsmch = NumTraits<Scalar>::epsilon();
0029     const Index n = x.size();
0030     eigen_assert(fvec.size()==n);
0031     Matrix< Scalar, Dynamic, 1 >  wa1(n);
0032     Matrix< Scalar, Dynamic, 1 >  wa2(n);
0033
0034     eps = sqrt((std::max)(epsfcn,epsmch));
0035     msum = ml + mu + 1;
0036     if (msum >= n) {
0037         /* computation of dense approximate jacobian. */
0038         for (j = 0; j < n; ++j) {
0039             temp = x[j];
0040             h = eps * abs(temp);
0041             if (h == 0.)
0042                 h = eps;
0043             x[j] = temp + h;
0044             iflag = Functor(x, wa1);
0045             if (iflag < 0)
0046                 return iflag;
0047             x[j] = temp;
0048             fjac.col(j) = (wa1-fvec)/h;
0049         }
0050
0051     }else {
0052         /* computation of banded approximate jacobian. */
0053         for (k = 0; k < msum; ++k) {
0054             for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
0055                 wa2[j] = x[j];
0056                 h = eps * abs(wa2[j]);
0057                 if (h == 0.) h = eps;
0058                 x[j] = wa2[j] + h;
0059             }
0060             iflag = Functor(x, wa1);
0061             if (iflag < 0)
0062                 return iflag;
0063             for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
0064                 x[j] = wa2[j];
0065                 h = eps * abs(wa2[j]);
0066                 if (h == 0.) h = eps;
0067                 fjac.col(j).setZero();
0068                 start = std::max<Index>(0,j-mu);
0069                 length = (std::min)(n-1, j+ml) - start + 1;
0070                 fjac.col(j).segment(start, length) = ( wa1.segment(start, length)-fvec.segment(start, length))/h;
0071             }
0072         }
0073     }
0074     return 0;
0075 }
0076
0077 } // end namespace internal
0078
0079 } // end namespace Eigen