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


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:24