fdjac1.h
Go to the documentation of this file.
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:01