Move work in progress Levenberg Marquardt module in unsupported

11 files changed, 2511 insertions, 1 deletions

diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt
index e961e72c5..e06f1238b 100644
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@@ -1,6 +1,6 @@
 set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials
                   FFT NonLinearOptimization SparseExtra IterativeSolvers
-                  NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct Splines
+                  NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct Splines LevenbergMarquardt
    )
 
 install(FILES
diff --git a/unsupported/Eigen/LevenbergMarquardt b/unsupported/Eigen/LevenbergMarquardt
new file mode 100644
index 000000000..5ce1a23f6
--- /dev/null
+++ b/unsupported/Eigen/LevenbergMarquardt
@@ -0,0 +1,46 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_LEVENBERGMARQUARDT_MODULE
+#define EIGEN_LEVENBERGMARQUARDT_MODULE
+
+// #include <vector>
+
+#include <Eigen/Core>
+#include <Eigen/Jacobi>
+#include <Eigen/QR>
+#include <unsupported/Eigen/NumericalDiff> 
+#ifdef EIGEN_SPQR_SUPPORT
+#include <Eigen/SPQRSupport>
+#endif
+
+/** \ingroup NonLinearOptimization_modules
+  * \defgroup LevenbergMarquardt_Module Levenberg-Marquardt optimization module
+  *
+  * \code
+  * #include </Eigen/LevenbergMarquardt>
+  * \endcode
+  *
+  * 
+  */
+
+#include "Eigen/SparseCore"
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+
+#include "src/LevenbergMarquardt/LMqrsolv.h"
+#include "src/LevenbergMarquardt/LMcovar.h"
+#include "src/LevenbergMarquardt/LMpar.h"
+
+#endif
+
+#include "src/LevenbergMarquardt/LevenbergMarquardt.h"
+#include "src/LevenbergMarquardt/LMonestep.h"
+
+
+#endif // EIGEN_LEVENBERGMARQUARDT_MODULE
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt b/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt
new file mode 100644
index 000000000..8513803ce
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_LevenbergMarquardt_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_LevenbergMarquardt_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LevenbergMarquardt COMPONENT Devel
+  )
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/CopyrightMINPACK.txt b/unsupported/Eigen/src/LevenbergMarquardt/CopyrightMINPACK.txt
new file mode 100644
index 000000000..ae7984dae
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/CopyrightMINPACK.txt
@@ -0,0 +1,52 @@
+Minpack Copyright Notice (1999) University of Chicago.  All rights reserved
+
+Redistribution and use in source and binary forms, with or
+without modification, are permitted provided that the
+following conditions are met:
+
+1. Redistributions of source code must retain the above
+copyright notice, this list of conditions and the following
+disclaimer.
+
+2. Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following
+disclaimer in the documentation and/or other materials
+provided with the distribution.
+
+3. The end-user documentation included with the
+redistribution, if any, must include the following
+acknowledgment:
+
+   "This product includes software developed by the
+   University of Chicago, as Operator of Argonne National
+   Laboratory.
+
+Alternately, this acknowledgment may appear in the software
+itself, if and wherever such third-party acknowledgments
+normally appear.
+
+4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
+WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
+UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
+THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
+OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
+OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
+OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
+USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
+THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
+DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
+UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
+BE CORRECTED.
+
+5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
+HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
+ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
+INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
+ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
+PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
+SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
+(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
+EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
+POSSIBILITY OF SUCH LOSS OR DAMAGES.
+
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMcovar.h b/unsupported/Eigen/src/LevenbergMarquardt/LMcovar.h
new file mode 100644
index 000000000..3210587e4
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LMcovar.h
@@ -0,0 +1,85 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This code initially comes from MINPACK whose original authors are:
+// Copyright Jorge More - Argonne National Laboratory
+// Copyright Burt Garbow - Argonne National Laboratory
+// Copyright Ken Hillstrom - Argonne National Laboratory
+//
+// This Source Code Form is subject to the terms of the Minpack license
+// (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
+
+#ifndef EIGEN_LMCOVAR_H
+#define EIGEN_LMCOVAR_H
+
+namespace Eigen { 
+
+namespace internal {
+
+template <typename Scalar>
+void covar(
+        Matrix< Scalar, Dynamic, Dynamic > &r,
+        const VectorXi& ipvt,
+        Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) )
+{
+    using std::abs;
+    typedef DenseIndex Index;
+    /* Local variables */
+    Index i, j, k, l, ii, jj;
+    bool sing;
+    Scalar temp;
+
+    /* Function Body */
+    const Index n = r.cols();
+    const Scalar tolr = tol * abs(r(0,0));
+    Matrix< Scalar, Dynamic, 1 > wa(n);
+    assert(ipvt.size()==n);
+
+    /* form the inverse of r in the full upper triangle of r. */
+    l = -1;
+    for (k = 0; k < n; ++k)
+        if (abs(r(k,k)) > tolr) {
+            r(k,k) = 1. / r(k,k);
+            for (j = 0; j <= k-1; ++j) {
+                temp = r(k,k) * r(j,k);
+                r(j,k) = 0.;
+                r.col(k).head(j+1) -= r.col(j).head(j+1) * temp;
+            }
+            l = k;
+        }
+
+    /* form the full upper triangle of the inverse of (r transpose)*r */
+    /* in the full upper triangle of r. */
+    for (k = 0; k <= l; ++k) {
+        for (j = 0; j <= k-1; ++j)
+            r.col(j).head(j+1) += r.col(k).head(j+1) * r(j,k);
+        r.col(k).head(k+1) *= r(k,k);
+    }
+
+    /* form the full lower triangle of the covariance matrix */
+    /* in the strict lower triangle of r and in wa. */
+    for (j = 0; j < n; ++j) {
+        jj = ipvt[j];
+        sing = j > l;
+        for (i = 0; i <= j; ++i) {
+            if (sing)
+                r(i,j) = 0.;
+            ii = ipvt[i];
+            if (ii > jj)
+                r(ii,jj) = r(i,j);
+            if (ii < jj)
+                r(jj,ii) = r(i,j);
+        }
+        wa[jj] = r(j,j);
+    }
+
+    /* symmetrize the covariance matrix in r. */
+    r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose();
+    r.diagonal() = wa;
+}
+
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_LMCOVAR_H
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h b/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h
new file mode 100644
index 000000000..aa0c02668
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LMonestep.h
@@ -0,0 +1,179 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
+//
+// This code initially comes from MINPACK whose original authors are:
+// Copyright Jorge More - Argonne National Laboratory
+// Copyright Burt Garbow - Argonne National Laboratory
+// Copyright Ken Hillstrom - Argonne National Laboratory
+//
+// This Source Code Form is subject to the terms of the Minpack license
+// (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
+
+#ifndef EIGEN_LMONESTEP_H
+#define EIGEN_LMONESTEP_H
+
+namespace Eigen {
+
+template<typename FunctorType>
+LevenbergMarquardtSpace::Status
+LevenbergMarquardt<FunctorType>::minimizeOneStep(FVectorType  &x)
+{
+  typedef typename FunctorType::JacobianType JacobianType; 
+  using std::abs;
+  using std::sqrt;
+  RealScalar temp, temp1,temp2; 
+  RealScalar ratio; 
+  RealScalar pnorm, xnorm, fnorm1, actred, dirder, prered;
+  assert(x.size()==n); // check the caller is not cheating us
+
+  temp = 0.0; xnorm = 0.0;
+  /* calculate the jacobian matrix. */
+  Index df_ret = m_functor.df(x, m_fjac);
+  if (df_ret<0)
+      return LevenbergMarquardtSpace::UserAsked;
+  if (df_ret>0)
+      // numerical diff, we evaluated the function df_ret times
+      m_nfev += df_ret;
+  else m_njev++;
+
+  /* compute the qr factorization of the jacobian. */
+  for (int j = 0; j < x.size(); ++j)
+    m_wa2(j) = m_fjac.col(j).norm(); 
+  //FIXME Implement bluenorm for sparse vectors
+//     m_wa2 = m_fjac.colwise().blueNorm();
+  QRSolver qrfac(m_fjac); //FIXME Check if the QR decomposition succeed
+  // Make a copy of the first factor with the associated permutation
+  JacobianType rfactor;
+  rfactor = qrfac.matrixQR();
+  m_permutation = (qrfac.colsPermutation());
+
+  /* on the first iteration and if external scaling is not used, scale according */
+  /* to the norms of the columns of the initial jacobian. */
+  if (m_iter == 1) {
+      if (!m_useExternalScaling)
+          for (Index j = 0; j < n; ++j)
+              m_diag[j] = (m_wa2[j]==0.)? 1. : m_wa2[j];
+
+      /* on the first iteration, calculate the norm of the scaled x */
+      /* and initialize the step bound m_delta. */
+      xnorm = m_diag.cwiseProduct(x).stableNorm();
+      m_delta = m_factor * xnorm;
+      if (m_delta == 0.)
+          m_delta = m_factor;
+  }
+
+  /* form (q transpose)*m_fvec and store the first n components in */
+  /* m_qtf. */
+  m_wa4 = m_fvec;
+  m_wa4 = qrfac.matrixQ().adjoint() * m_fvec; 
+  m_qtf = m_wa4.head(n);
+
+  /* compute the norm of the scaled gradient. */
+  m_gnorm = 0.;
+  if (m_fnorm != 0.)
+      for (Index j = 0; j < n; ++j)
+          if (m_wa2[m_permutation.indices()[j]] != 0.)
+              m_gnorm = (std::max)(m_gnorm, abs( rfactor.col(j).head(j+1).dot(m_qtf.head(j+1)/m_fnorm) / m_wa2[m_permutation.indices()[j]]));
+
+  /* test for convergence of the gradient norm. */
+  if (m_gnorm <= m_gtol)
+      return LevenbergMarquardtSpace::CosinusTooSmall;
+
+  /* rescale if necessary. */
+  if (!m_useExternalScaling)
+      m_diag = m_diag.cwiseMax(m_wa2);
+
+  do {
+    /* determine the levenberg-marquardt parameter. */
+    internal::lmpar2(qrfac, m_diag, m_qtf, m_delta, m_par, m_wa1);
+
+    /* store the direction p and x + p. calculate the norm of p. */
+    m_wa1 = -m_wa1;
+    m_wa2 = x + m_wa1;
+    pnorm = m_diag.cwiseProduct(m_wa1).stableNorm();
+
+    /* on the first iteration, adjust the initial step bound. */
+    if (m_iter == 1)
+        m_delta = (std::min)(m_delta,pnorm);
+
+    /* evaluate the function at x + p and calculate its norm. */
+    if ( m_functor(m_wa2, m_wa4) < 0)
+        return LevenbergMarquardtSpace::UserAsked;
+    ++m_nfev;
+    fnorm1 = m_wa4.stableNorm();
+
+    /* compute the scaled actual reduction. */
+    actred = -1.;
+    if (Scalar(.1) * fnorm1 < m_fnorm)
+        actred = 1. - internal::abs2(fnorm1 / m_fnorm);
+
+    /* compute the scaled predicted reduction and */
+    /* the scaled directional derivative. */
+    m_wa3 = rfactor.template triangularView<Upper>() * (m_permutation.inverse() *m_wa1);
+    temp1 = internal::abs2(m_wa3.stableNorm() / m_fnorm);
+    temp2 = internal::abs2(sqrt(m_par) * pnorm / m_fnorm);
+    prered = temp1 + temp2 / Scalar(.5);
+    dirder = -(temp1 + temp2);
+
+    /* compute the ratio of the actual to the predicted */
+    /* reduction. */
+    ratio = 0.;
+    if (prered != 0.)
+        ratio = actred / prered;
+
+    /* update the step bound. */
+    if (ratio <= Scalar(.25)) {
+        if (actred >= 0.)
+            temp = RealScalar(.5);
+        if (actred < 0.)
+            temp = RealScalar(.5) * dirder / (dirder + RealScalar(.5) * actred);
+        if (RealScalar(.1) * fnorm1 >= m_fnorm || temp < RealScalar(.1))
+            temp = Scalar(.1);
+        /* Computing MIN */
+        m_delta = temp * (std::min)(m_delta, pnorm / RealScalar(.1));
+        m_par /= temp;
+    } else if (!(m_par != 0. && ratio < RealScalar(.75))) {
+        m_delta = pnorm / RealScalar(.5);
+        m_par = RealScalar(.5) * m_par;
+    }
+
+    /* test for successful iteration. */
+    if (ratio >= RealScalar(1e-4)) {
+        /* successful iteration. update x, m_fvec, and their norms. */
+        x = m_wa2;
+        m_wa2 = m_diag.cwiseProduct(x);
+        m_fvec = m_wa4;
+        xnorm = m_wa2.stableNorm();
+        m_fnorm = fnorm1;
+        ++m_iter;
+    }
+
+    /* tests for convergence. */
+    if (abs(actred) <= m_ftol && prered <= m_ftol && Scalar(.5) * ratio <= 1. && m_delta <= m_xtol * xnorm)
+        return LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall;
+    if (abs(actred) <= m_ftol && prered <= m_ftol && Scalar(.5) * ratio <= 1.)
+        return LevenbergMarquardtSpace::RelativeReductionTooSmall;
+    if (m_delta <= m_xtol * xnorm)
+        return LevenbergMarquardtSpace::RelativeErrorTooSmall;
+
+    /* tests for termination and stringent tolerances. */
+    if (m_nfev >= m_maxfev)
+        return LevenbergMarquardtSpace::TooManyFunctionEvaluation;
+    if (abs(actred) <= NumTraits<Scalar>::epsilon() && prered <= NumTraits<Scalar>::epsilon() && Scalar(.5) * ratio <= 1.)
+        return LevenbergMarquardtSpace::FtolTooSmall;
+    if (m_delta <= NumTraits<Scalar>::epsilon() * xnorm)
+        return LevenbergMarquardtSpace::XtolTooSmall;
+    if (m_gnorm <= NumTraits<Scalar>::epsilon())
+        return LevenbergMarquardtSpace::GtolTooSmall;
+
+  } while (ratio < Scalar(1e-4));
+
+  return LevenbergMarquardtSpace::Running;
+}
+
+  
+} // end namespace Eigen
+
+#endif // EIGEN_LMONESTEP_H
\ No newline at end of file
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
new file mode 100644
index 000000000..dc60ca05a
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LMpar.h
@@ -0,0 +1,160 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This code initially comes from MINPACK whose original authors are:
+// Copyright Jorge More - Argonne National Laboratory
+// Copyright Burt Garbow - Argonne National Laboratory
+// Copyright Ken Hillstrom - Argonne National Laboratory
+//
+// This Source Code Form is subject to the terms of the Minpack license
+// (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
+
+#ifndef EIGEN_LMPAR_H
+#define EIGEN_LMPAR_H
+
+namespace Eigen {
+
+namespace internal {
+  
+  template <typename QRSolver, typename VectorType>
+    void lmpar2(
+    const QRSolver &qr,
+    const VectorType  &diag,
+    const VectorType  &qtb,
+    typename VectorType::Scalar m_delta,
+    typename VectorType::Scalar &par,
+    VectorType  &x)
+
+  {
+    using std::sqrt;
+    using std::abs;
+    typedef typename QRSolver::MatrixType MatrixType;
+    typedef typename QRSolver::Scalar Scalar;
+    typedef typename QRSolver::Index Index;
+
+    /* Local variables */
+    Index j;
+    Scalar fp;
+    Scalar parc, parl;
+    Index iter;
+    Scalar temp, paru;
+    Scalar gnorm;
+    Scalar dxnorm;
+
+
+    /* Function Body */
+    const Scalar dwarf = (std::numeric_limits<Scalar>::min)();
+    const Index n = qr.matrixQR().cols();
+    assert(n==diag.size());
+    assert(n==qtb.size());
+
+    VectorType  wa1, wa2;
+
+    /* compute and store in x the gauss-newton direction. if the */
+    /* jacobian is rank-deficient, obtain a least squares solution. */
+
+    //    const Index rank = qr.nonzeroPivots(); // exactly double(0.)
+    const Index rank = qr.rank(); // use a threshold
+    wa1 = qtb;
+    wa1.tail(n-rank).setZero();
+    //FIXME There is no solve in place for sparse triangularView
+    //qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solveInPlace(wa1.head(rank));
+    wa1.head(rank) = qr.matrixQR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(qtb.head(rank));
+
+    x = qr.colsPermutation()*wa1;
+
+    /* initialize the iteration counter. */
+    /* evaluate the function at the origin, and test */
+    /* for acceptance of the gauss-newton direction. */
+    iter = 0;
+    wa2 = diag.cwiseProduct(x);
+    dxnorm = wa2.blueNorm();
+    fp = dxnorm - m_delta;
+    if (fp <= Scalar(0.1) * m_delta) {
+      par = 0;
+      return;
+    }
+
+    /* if the jacobian is not rank deficient, the newton */
+    /* step provides a lower bound, parl, for the zero of */
+    /* the function. otherwise set this bound to zero. */
+    parl = 0.;
+    if (rank==n) {
+      wa1 = qr.colsPermutation().inverse() *  diag.cwiseProduct(wa2)/dxnorm;
+      qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1);
+      temp = wa1.blueNorm();
+      parl = fp / m_delta / temp / temp;
+    }
+
+    /* calculate an upper bound, paru, for the zero of the function. */
+    for (j = 0; j < n; ++j)
+      wa1[j] = qr.matrixQR().col(j).head(j+1).dot(qtb.head(j+1)) / diag[qr.colsPermutation().indices()(j)];
+
+    gnorm = wa1.stableNorm();
+    paru = gnorm / m_delta;
+    if (paru == 0.)
+      paru = dwarf / (std::min)(m_delta,Scalar(0.1));
+
+    /* if the input par lies outside of the interval (parl,paru), */
+    /* set par to the closer endpoint. */
+    par = (std::max)(par,parl);
+    par = (std::min)(par,paru);
+    if (par == 0.)
+      par = gnorm / dxnorm;
+
+    /* beginning of an iteration. */
+    MatrixType s;
+    s = qr.matrixQR();
+    while (true) {
+      ++iter;
+
+      /* evaluate the function at the current value of par. */
+      if (par == 0.)
+        par = (std::max)(dwarf,Scalar(.001) * paru); /* Computing MAX */
+      wa1 = sqrt(par)* diag;
+
+      VectorType sdiag(n);
+      lmqrsolv(s, qr.colsPermutation(), wa1, qtb, x, sdiag);
+
+      wa2 = diag.cwiseProduct(x);
+      dxnorm = wa2.blueNorm();
+      temp = fp;
+      fp = dxnorm - m_delta;
+
+      /* if the function is small enough, accept the current value */
+      /* of par. also test for the exceptional cases where parl */
+      /* is zero or the number of iterations has reached 10. */
+      if (abs(fp) <= Scalar(0.1) * m_delta || (parl == 0. && fp <= temp && temp < 0.) || iter == 10)
+        break;
+
+      /* compute the newton correction. */
+      wa1 = qr.colsPermutation().inverse() * diag.cwiseProduct(wa2/dxnorm);
+      // we could almost use this here, but the diagonal is outside qr, in sdiag[]
+      // qr.matrixQR().topLeftCorner(n, n).transpose().template triangularView<Lower>().solveInPlace(wa1);
+      for (j = 0; j < n; ++j) {
+        wa1[j] /= sdiag[j];
+        temp = wa1[j];
+        for (Index i = j+1; i < n; ++i)
+          wa1[i] -= s.coeff(i,j) * temp;
+      }
+      temp = wa1.blueNorm();
+      parc = fp / m_delta / temp / temp;
+
+      /* depending on the sign of the function, update parl or paru. */
+      if (fp > 0.)
+        parl = (std::max)(parl,par);
+      if (fp < 0.)
+        paru = (std::min)(paru,par);
+
+      /* compute an improved estimate for par. */
+      par = (std::max)(parl,par+parc);
+    }
+    if (iter == 0)
+      par = 0.;
+    return;
+  }
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_LMPAR_H
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h b/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h
new file mode 100644
index 000000000..ed6f97fe8
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LMqrsolv.h
@@ -0,0 +1,189 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
+// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
+//
+// This code initially comes from MINPACK whose original authors are:
+// Copyright Jorge More - Argonne National Laboratory
+// Copyright Burt Garbow - Argonne National Laboratory
+// Copyright Ken Hillstrom - Argonne National Laboratory
+//
+// This Source Code Form is subject to the terms of the Minpack license
+// (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
+
+#ifndef EIGEN_LMQRSOLV_H
+#define EIGEN_LMQRSOLV_H
+
+namespace Eigen { 
+
+namespace internal {
+
+template <typename Scalar,int SizeAtCompileTime, typename Index>
+void lmqrsolv(
+  Matrix<Scalar,SizeAtCompileTime,SizeAtCompileTime> &s,
+  const PermutationMatrix<Dynamic,Dynamic,Index> &iPerm,
+  const Matrix<Scalar,Dynamic,1> &diag,
+  const Matrix<Scalar,Dynamic,1> &qtb,
+  Matrix<Scalar,Dynamic,1> &x,
+  Matrix<Scalar,Dynamic,1> &sdiag)
+{
+
+    /* Local variables */
+    Index i, j, k, l;
+    Scalar temp;
+    Index n = s.cols();
+    Matrix<Scalar,Dynamic,1>  wa(n);
+    JacobiRotation<Scalar> givens;
+
+    /* Function Body */
+    // the following will only change the lower triangular part of s, including
+    // the diagonal, though the diagonal is restored afterward
+
+    /*     copy r and (q transpose)*b to preserve input and initialize s. */
+    /*     in particular, save the diagonal elements of r in x. */
+    x = s.diagonal();
+    wa = qtb;
+    
+   
+    s.topLeftCorner(n,n).template triangularView<StrictlyLower>() = s.topLeftCorner(n,n).transpose();
+    /*     eliminate the diagonal matrix d using a givens rotation. */
+    for (j = 0; j < n; ++j) {
+
+        /*        prepare the row of d to be eliminated, locating the */
+        /*        diagonal element using p from the qr factorization. */
+        l = iPerm.indices()(j);
+        if (diag[l] == 0.)
+            break;
+        sdiag.tail(n-j).setZero();
+        sdiag[j] = diag[l];
+
+        /*        the transformations to eliminate the row of d */
+        /*        modify only a single element of (q transpose)*b */
+        /*        beyond the first n, which is initially zero. */
+        Scalar qtbpj = 0.;
+        for (k = j; k < n; ++k) {
+            /*           determine a givens rotation which eliminates the */
+            /*           appropriate element in the current row of d. */
+            givens.makeGivens(-s(k,k), sdiag[k]);
+
+            /*           compute the modified diagonal element of r and */
+            /*           the modified element of ((q transpose)*b,0). */
+            s(k,k) = givens.c() * s(k,k) + givens.s() * sdiag[k];
+            temp = givens.c() * wa[k] + givens.s() * qtbpj;
+            qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj;
+            wa[k] = temp;
+
+            /*           accumulate the tranformation in the row of s. */
+            for (i = k+1; i<n; ++i) {
+                temp = givens.c() * s(i,k) + givens.s() * sdiag[i];
+                sdiag[i] = -givens.s() * s(i,k) + givens.c() * sdiag[i];
+                s(i,k) = temp;
+            }
+        }
+    }
+  
+    /*     solve the triangular system for z. if the system is */
+    /*     singular, then obtain a least squares solution. */
+    Index nsing;
+    for(nsing=0; nsing<n && sdiag[nsing]!=0; nsing++) {}
+
+    wa.tail(n-nsing).setZero();
+    s.topLeftCorner(nsing, nsing).transpose().template triangularView<Upper>().solveInPlace(wa.head(nsing));
+  
+    // restore
+    sdiag = s.diagonal();
+    s.diagonal() = x;
+
+    /* permute the components of z back to components of x. */
+    x = iPerm * wa; 
+}
+
+template <typename Scalar, int _Options, typename Index>
+void lmqrsolv(
+  SparseMatrix<Scalar,_Options,Index> &s,
+  const PermutationMatrix<Dynamic,Dynamic> &iPerm,
+  const Matrix<Scalar,Dynamic,1> &diag,
+  const Matrix<Scalar,Dynamic,1> &qtb,
+  Matrix<Scalar,Dynamic,1> &x,
+  Matrix<Scalar,Dynamic,1> &sdiag)
+{
+  /* Local variables */
+  typedef SparseMatrix<Scalar,RowMajor,Index> FactorType;
+    Index i, j, k, l;
+    Scalar temp;
+    Index n = s.cols();
+    Matrix<Scalar,Dynamic,1>  wa(n);
+    JacobiRotation<Scalar> givens;
+
+    /* Function Body */
+    // the following will only change the lower triangular part of s, including
+    // the diagonal, though the diagonal is restored afterward
+
+    /*     copy r and (q transpose)*b to preserve input and initialize R. */
+    wa = qtb;
+    FactorType R(s);
+    // Eliminate the diagonal matrix d using a givens rotation
+    for (j = 0; j < n; ++j)
+    {
+      // Prepare the row of d to be eliminated, locating the 
+      // diagonal element using p from the qr factorization
+      l = iPerm.indices()(j);
+      if (diag(l) == Scalar(0)) 
+        break; 
+      sdiag.tail(n-j).setZero();
+      sdiag[j] = diag[l];
+      // the transformations to eliminate the row of d
+      // modify only a single element of (q transpose)*b
+      // beyond the first n, which is initially zero. 
+      
+      Scalar qtbpj = 0; 
+      // Browse the nonzero elements of row j of the upper triangular s
+      for (k = j; k < n; ++k)
+      {
+        typename FactorType::InnerIterator itk(R,k);
+        for (; itk; ++itk){
+          if (itk.index() < k) continue;
+          else break;
+        }
+        //At this point, we have the diagonal element R(k,k)
+        // Determine a givens rotation which eliminates 
+        // the appropriate element in the current row of d
+        givens.makeGivens(-itk.value(), sdiag(k));
+        
+        // Compute the modified diagonal element of r and 
+        // the modified element of ((q transpose)*b,0).
+        itk.valueRef() = givens.c() * itk.value() + givens.s() * sdiag(k);
+        temp = givens.c() * wa(k) + givens.s() * qtbpj; 
+        qtbpj = -givens.s() * wa(k) + givens.c() * qtbpj;
+        wa(k) = temp;
+        
+        // Accumulate the transformation in the remaining k row/column of R
+        for (++itk; itk; ++itk)
+        {
+          i = itk.index();
+          temp = givens.c() *  itk.value() + givens.s() * sdiag(i);
+          sdiag(i) = -givens.s() * itk.value() + givens.c() * sdiag(i);
+          itk.valueRef() = temp;
+        }
+      }
+    }
+    
+    // Solve the triangular system for z. If the system is 
+    // singular, then obtain a least squares solution
+    Index nsing;
+    for(nsing = 0; nsing<n && sdiag(nsing) !=0; nsing++) {}
+    
+    wa.tail(n-nsing).setZero();
+//     x = wa; 
+    wa.head(nsing) = R.topLeftCorner(nsing,nsing).template triangularView<Upper>().solve/*InPlace*/(wa.head(nsing));
+    
+    sdiag = R.diagonal();
+    // Permute the components of z back to components of x
+    x = iPerm * wa; 
+}
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_LMQRSOLV_H
diff --git a/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
new file mode 100644
index 000000000..e45e73ab5
--- /dev/null
+++ b/unsupported/Eigen/src/LevenbergMarquardt/LevenbergMarquardt.h
@@ -0,0 +1,344 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
+// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
+//
+// The algorithm of this class initially comes from MINPACK whose original authors are:
+// Copyright Jorge More - Argonne National Laboratory
+// Copyright Burt Garbow - Argonne National Laboratory
+// Copyright Ken Hillstrom - Argonne National Laboratory
+//
+// This Source Code Form is subject to the terms of the Minpack license
+// (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_LEVENBERGMARQUARDT_H
+#define EIGEN_LEVENBERGMARQUARDT_H
+
+
+namespace Eigen {
+namespace LevenbergMarquardtSpace {
+    enum Status {
+        NotStarted = -2,
+        Running = -1,
+        ImproperInputParameters = 0,
+        RelativeReductionTooSmall = 1,
+        RelativeErrorTooSmall = 2,
+        RelativeErrorAndReductionTooSmall = 3,
+        CosinusTooSmall = 4,
+        TooManyFunctionEvaluation = 5,
+        FtolTooSmall = 6,
+        XtolTooSmall = 7,
+        GtolTooSmall = 8,
+        UserAsked = 9
+    };
+}
+
+template <typename _Scalar, int NX=Dynamic, int NY=Dynamic>
+struct DenseFunctor
+{
+  typedef _Scalar Scalar;
+  enum {
+    InputsAtCompileTime = NX,
+    ValuesAtCompileTime = NY
+  };
+  typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
+  typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
+  typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
+  typedef ColPivHouseholderQR<JacobianType> QRSolver;
+  const int m_inputs, m_values;
+
+  DenseFunctor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
+  DenseFunctor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
+
+  int inputs() const { return m_inputs; }
+  int values() const { return m_values; }
+
+  //int operator()(const InputType &x, ValueType& fvec) { }
+  // should be defined in derived classes
+  
+  //int df(const InputType &x, JacobianType& fjac) { }
+  // should be defined in derived classes
+};
+
+#ifdef EIGEN_SPQR_SUPPORT
+template <typename _Scalar, typename _Index>
+struct SparseFunctor
+{
+  typedef _Scalar Scalar;
+  typedef _Index Index;
+  typedef Matrix<Scalar,Dynamic,1> InputType;
+  typedef Matrix<Scalar,Dynamic,1> ValueType;
+  typedef SparseMatrix<Scalar, ColMajor, Index> JacobianType;
+  typedef SPQR<JacobianType> QRSolver;
+  
+  SparseFunctor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
+
+  int inputs() const { return m_inputs; }
+  int values() const { return m_values; }
+  
+  const int m_inputs, m_values;
+  //int operator()(const InputType &x, ValueType& fvec) { }
+  // to be defined in the functor
+  
+  //int df(const InputType &x, JacobianType& fjac) { }
+  // to be defined in the functor if no automatic differentiation
+  
+};
+#endif
+namespace internal {
+template <typename QRSolver, typename VectorType>
+void lmpar2(const QRSolver &qr, const VectorType  &diag, const VectorType  &qtb,
+	    typename VectorType::Scalar m_delta, typename VectorType::Scalar &par,
+	    VectorType  &x);
+    }
+/**
+  * \ingroup NonLinearOptimization_Module
+  * \brief Performs non linear optimization over a non-linear function,
+  * using a variant of the Levenberg Marquardt algorithm.
+  *
+  * Check wikipedia for more information.
+  * http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
+  */
+template<typename _FunctorType>
+class LevenbergMarquardt
+{
+  public:
+    typedef _FunctorType FunctorType;
+    typedef typename FunctorType::QRSolver QRSolver;
+    typedef typename FunctorType::JacobianType JacobianType;
+    typedef typename JacobianType::Scalar Scalar;
+    typedef typename JacobianType::RealScalar RealScalar; 
+    typedef typename JacobianType::Index Index;
+    typedef typename QRSolver::Index PermIndex;
+    typedef Matrix<Scalar,Dynamic,1> FVectorType;
+    typedef PermutationMatrix<Dynamic,Dynamic> PermutationType;
+  public:
+    LevenbergMarquardt(FunctorType& functor) 
+    : m_functor(functor),m_nfev(0),m_njev(0),m_fnorm(0.0),m_gnorm(0)
+    {
+      resetParameters();
+      m_useExternalScaling=false; 
+    }
+    
+    LevenbergMarquardtSpace::Status minimize(FVectorType &x);
+    LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
+    LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
+    LevenbergMarquardtSpace::Status lmder1(
+      FVectorType  &x, 
+      const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon())
+    );
+    static LevenbergMarquardtSpace::Status lmdif1(
+            FunctorType &functor,
+            FVectorType  &x,
+            Index *nfev,
+            const Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon())
+            );
+    
+    /** Sets the default parameters */
+    void resetParameters() 
+    { 
+      m_factor = 100.; 
+      m_maxfev = 400; 
+      m_ftol = std::sqrt(NumTraits<RealScalar>::epsilon());
+      m_xtol = std::sqrt(NumTraits<RealScalar>::epsilon());
+      m_gtol = 0. ; 
+      m_epsfcn = 0. ;
+    }
+    
+    /** Sets the tolerance for the norm of the solution vector*/
+    void setXtol(RealScalar xtol) { m_xtol = xtol; }
+    
+    /** Sets the tolerance for the norm of the vector function*/
+    void setFtol(RealScalar ftol) { m_ftol = ftol; }
+    
+    /** Sets the tolerance for the norm of the gradient of the error vector*/
+    void setGtol(RealScalar gtol) { m_gtol = gtol; }
+    
+    /** Sets the step bound for the diagonal shift */
+    void setFactor(RealScalar factor) { m_factor = factor; }    
+    
+    /** Sets the error precision  */
+    void setEpsilon (RealScalar epsfcn) { m_epsfcn = epsfcn; }
+    
+    /** Sets the maximum number of function evaluation */
+    void setMaxfev(Index maxfev) {m_maxfev = maxfev; }
+    
+    /** Use an external Scaling. If set to true, pass a nonzero diagonal to diag() */
+    void setExternalScaling(bool value) {m_useExternalScaling  = value; }
+    
+    /** Get a reference to the diagonal of the jacobian */
+    FVectorType& diag() {return m_diag; }
+    
+    /** Number of iterations performed */
+    Index iterations() { return m_iter; }
+    
+    /** Number of functions evaluation */
+    Index nfev() { return m_nfev; }
+    
+    /** Number of jacobian evaluation */
+    Index njev() { return m_njev; }
+    
+    /** Norm of current vector function */
+    RealScalar fnorm() {return m_fnorm; }
+    
+    /** Norm of the gradient of the error */
+    RealScalar gnorm() {return m_gnorm; }
+    
+    /** the LevenbergMarquardt parameter */
+    RealScalar lm_param(void) { return m_par; }
+    
+    /** reference to the  current vector function 
+     */
+    FVectorType& fvec() {return m_fvec; }
+    
+    /** reference to the matrix where the current Jacobian matrix is stored
+     */
+    JacobianType& fjac() {return m_fjac; }
+    
+    /** the permutation used
+     */
+    PermutationType permutation() {return m_permutation; }
+    
+  private:
+    JacobianType m_fjac; 
+    FunctorType &m_functor;
+    FVectorType m_fvec, m_qtf, m_diag; 
+    Index n;
+    Index m; 
+    Index m_nfev;
+    Index m_njev; 
+    RealScalar m_fnorm; // Norm of the current vector function
+    RealScalar m_gnorm; //Norm of the gradient of the error 
+    RealScalar m_factor; //
+    Index m_maxfev; // Maximum number of function evaluation
+    RealScalar m_ftol; //Tolerance in the norm of the vector function
+    RealScalar m_xtol; // 
+    RealScalar m_gtol; //tolerance of the norm of the error gradient
+    RealScalar m_epsfcn; //
+    Index m_iter; // Number of iterations performed
+    RealScalar m_delta;
+    bool m_useExternalScaling;
+    PermutationType m_permutation;
+    FVectorType m_wa1, m_wa2, m_wa3, m_wa4; //Temporary vectors
+    RealScalar m_par;
+};
+
+template<typename FunctorType>
+LevenbergMarquardtSpace::Status
+LevenbergMarquardt<FunctorType>::minimize(FVectorType  &x)
+{
+    LevenbergMarquardtSpace::Status status = minimizeInit(x);
+    if (status==LevenbergMarquardtSpace::ImproperInputParameters)
+        return status;
+    do {
+//       std::cout << " uv " << x.transpose() << "\n";
+        status = minimizeOneStep(x);
+    } while (status==LevenbergMarquardtSpace::Running);
+    return status;
+}
+
+template<typename FunctorType>
+LevenbergMarquardtSpace::Status
+LevenbergMarquardt<FunctorType>::minimizeInit(FVectorType  &x)
+{
+    n = x.size();
+    m = m_functor.values();
+
+    m_wa1.resize(n); m_wa2.resize(n); m_wa3.resize(n);
+    m_wa4.resize(m);
+    m_fvec.resize(m);
+    //FIXME Sparse Case : Allocate space for the jacobian
+    m_fjac.resize(m, n);
+//     m_fjac.reserve(VectorXi::Constant(n,5)); // FIXME Find a better alternative
+    if (!m_useExternalScaling)
+        m_diag.resize(n);
+    assert( (!m_useExternalScaling || m_diag.size()==n) || "When m_useExternalScaling is set, the caller must provide a valid 'm_diag'");
+    m_qtf.resize(n);
+
+    /* Function Body */
+    m_nfev = 0;
+    m_njev = 0;
+
+    /*     check the input parameters for errors. */
+    if (n <= 0 || m < n || m_ftol < 0. || m_xtol < 0. || m_gtol < 0. || m_maxfev <= 0 || m_factor <= 0.)
+        return LevenbergMarquardtSpace::ImproperInputParameters;
+
+    if (m_useExternalScaling)
+        for (Index j = 0; j < n; ++j)
+            if (m_diag[j] <= 0.)
+                return LevenbergMarquardtSpace::ImproperInputParameters;
+
+    /*     evaluate the function at the starting point */
+    /*     and calculate its norm. */
+    m_nfev = 1;
+    if ( m_functor(x, m_fvec) < 0)
+        return LevenbergMarquardtSpace::UserAsked;
+    m_fnorm = m_fvec.stableNorm();
+
+    /*     initialize levenberg-marquardt parameter and iteration counter. */
+    m_par = 0.;
+    m_iter = 1;
+
+    return LevenbergMarquardtSpace::NotStarted;
+}
+
+template<typename FunctorType>
+LevenbergMarquardtSpace::Status
+LevenbergMarquardt<FunctorType>::lmder1(
+        FVectorType  &x,
+        const Scalar tol
+        )
+{
+    n = x.size();
+    m = m_functor.values();
+
+    /* check the input parameters for errors. */
+    if (n <= 0 || m < n || tol < 0.)
+        return LevenbergMarquardtSpace::ImproperInputParameters;
+
+    resetParameters();
+    m_ftol = tol;
+    m_xtol = tol;
+    m_maxfev = 100*(n+1);
+
+    return minimize(x);
+}
+
+
+template<typename FunctorType>
+LevenbergMarquardtSpace::Status
+LevenbergMarquardt<FunctorType>::lmdif1(
+        FunctorType &functor,
+        FVectorType  &x,
+        Index *nfev,
+        const Scalar tol
+        )
+{
+    Index n = x.size();
+    Index m = functor.values();
+
+    /* check the input parameters for errors. */
+    if (n <= 0 || m < n || tol < 0.)
+        return LevenbergMarquardtSpace::ImproperInputParameters;
+
+    NumericalDiff<FunctorType> numDiff(functor);
+    // embedded LevenbergMarquardt
+    LevenbergMarquardt<NumericalDiff<FunctorType> > lm(numDiff);
+    lm.setFtol(tol);
+    lm.setXtol(tol);
+    lm.setMaxfev(200*(n+1));
+
+    LevenbergMarquardtSpace::Status info = LevenbergMarquardtSpace::Status(lm.minimize(x));
+    if (nfev)
+        * nfev = lm.nfev();
+    return info;
+}
+
+} // end namespace Eigen
+
+#endif // EIGEN_LEVENBERGMARQUARDT_H
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index cc690160d..0e0c8a6bf 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -86,3 +86,4 @@ ei_add_test(kronecker_product)
 ei_add_test(splines)
 ei_add_test(gmres)
 ei_add_test(minres)
+ei_add_test(levenberg_marquardt)
diff --git a/unsupported/test/levenberg_marquardt.cpp b/unsupported/test/levenberg_marquardt.cpp
new file mode 100644
index 000000000..c7061f017
--- /dev/null
+++ b/unsupported/test/levenberg_marquardt.cpp
@@ -0,0 +1,1448 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
+// Copyright (C) 2012 desire Nuentsa <desire.nuentsa_wakam@inria.fr
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+
+#include <stdio.h>
+
+#include "main.h"
+#include <Eigen/LevenbergMarquardt>
+
+// This disables some useless Warnings on MSVC.
+// It is intended to be done for this test only.
+#include <Eigen/src/Core/util/DisableStupidWarnings.h>
+
+using std::sqrt;
+
+struct lmder_functor : DenseFunctor<double>
+{
+    lmder_functor(void): DenseFunctor<double>(3,15) {}
+    int operator()(const VectorXd &x, VectorXd &fvec) const
+    {
+        double tmp1, tmp2, tmp3;
+        static const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
+            3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
+
+        for (int i = 0; i < values(); i++)
+        {
+            tmp1 = i+1;
+            tmp2 = 16 - i - 1;
+            tmp3 = (i>=8)? tmp2 : tmp1;
+            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
+        }
+        return 0;
+    }
+
+    int df(const VectorXd &x, MatrixXd &fjac) const
+    {
+        double tmp1, tmp2, tmp3, tmp4;
+        for (int i = 0; i < values(); i++)
+        {
+            tmp1 = i+1;
+            tmp2 = 16 - i - 1;
+            tmp3 = (i>=8)? tmp2 : tmp1;
+            tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
+            fjac(i,0) = -1;
+            fjac(i,1) = tmp1*tmp2/tmp4;
+            fjac(i,2) = tmp1*tmp3/tmp4;
+        }
+        return 0;
+    }
+};
+
+void testLmder1()
+{
+  int n=3, info;
+
+  VectorXd x;
+
+  /* the following starting values provide a rough fit. */
+  x.setConstant(n, 1.);
+
+  // do the computation
+  lmder_functor functor;
+  LevenbergMarquardt<lmder_functor> lm(functor);
+  info = lm.lmder1(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 6);
+  VERIFY_IS_EQUAL(lm.njev(), 5);
+
+  // check norm
+  VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596);
+
+  // check x
+  VectorXd x_ref(n);
+  x_ref << 0.08241058, 1.133037, 2.343695;
+  VERIFY_IS_APPROX(x, x_ref);
+}
+
+void testLmder()
+{
+  const int m=15, n=3;
+  int info;
+  double fnorm, covfac;
+  VectorXd x;
+
+  /* the following starting values provide a rough fit. */
+  x.setConstant(n, 1.);
+
+  // do the computation
+  lmder_functor functor;
+  LevenbergMarquardt<lmder_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return values
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 6);
+  VERIFY_IS_EQUAL(lm.njev(), 5);
+
+  // check norm
+  fnorm = lm.fvec().blueNorm();
+  VERIFY_IS_APPROX(fnorm, 0.09063596);
+
+  // check x
+  VectorXd x_ref(n);
+  x_ref << 0.08241058, 1.133037, 2.343695;
+  VERIFY_IS_APPROX(x, x_ref);
+
+  // check covariance
+  covfac = fnorm*fnorm/(m-n);
+  internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm
+
+  MatrixXd cov_ref(n,n);
+  cov_ref <<
+      0.0001531202,   0.002869941,  -0.002656662,
+      0.002869941,    0.09480935,   -0.09098995,
+      -0.002656662,   -0.09098995,    0.08778727;
+
+//  std::cout << fjac*covfac << std::endl;
+
+  MatrixXd cov;
+  cov =  covfac*lm.fjac().topLeftCorner<n,n>();
+  VERIFY_IS_APPROX( cov, cov_ref);
+  // TODO: why isn't this allowed ? :
+  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
+}
+
+struct lmdif_functor : DenseFunctor<double>
+{
+    lmdif_functor(void) : DenseFunctor<double>(3,15) {}
+    int operator()(const VectorXd &x, VectorXd &fvec) const
+    {
+        int i;
+        double tmp1,tmp2,tmp3;
+        static const double y[15]={1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1,
+            3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0};
+
+        assert(x.size()==3);
+        assert(fvec.size()==15);
+        for (i=0; i<15; i++)
+        {
+            tmp1 = i+1;
+            tmp2 = 15 - i;
+            tmp3 = tmp1;
+
+            if (i >= 8) tmp3 = tmp2;
+            fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
+        }
+        return 0;
+    }
+};
+
+void testLmdif1()
+{
+  const int n=3;
+  int info;
+
+  VectorXd x(n), fvec(15);
+
+  /* the following starting values provide a rough fit. */
+  x.setConstant(n, 1.);
+
+  // do the computation
+  lmdif_functor functor;
+  DenseIndex nfev;
+  info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(nfev, 26);
+
+  // check norm
+  functor(x, fvec);
+  VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596);
+
+  // check x
+  VectorXd x_ref(n);
+  x_ref << 0.0824106, 1.1330366, 2.3436947;
+  VERIFY_IS_APPROX(x, x_ref);
+
+}
+
+void testLmdif()
+{
+  const int m=15, n=3;
+  int info;
+  double fnorm, covfac;
+  VectorXd x(n);
+
+  /* the following starting values provide a rough fit. */
+  x.setConstant(n, 1.);
+
+  // do the computation
+  lmdif_functor functor;
+  NumericalDiff<lmdif_functor> numDiff(functor);
+  LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff);
+  info = lm.minimize(x);
+
+  // check return values
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 26);
+
+  // check norm
+  fnorm = lm.fvec().blueNorm();
+  VERIFY_IS_APPROX(fnorm, 0.09063596);
+
+  // check x
+  VectorXd x_ref(n);
+  x_ref << 0.08241058, 1.133037, 2.343695;
+  VERIFY_IS_APPROX(x, x_ref);
+
+  // check covariance
+  covfac = fnorm*fnorm/(m-n);
+  internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm
+
+  MatrixXd cov_ref(n,n);
+  cov_ref <<
+      0.0001531202,   0.002869942,  -0.002656662,
+      0.002869942,    0.09480937,   -0.09098997,
+      -0.002656662,   -0.09098997,    0.08778729;
+
+//  std::cout << fjac*covfac << std::endl;
+
+  MatrixXd cov;
+  cov =  covfac*lm.fjac().topLeftCorner<n,n>();
+  VERIFY_IS_APPROX( cov, cov_ref);
+  // TODO: why isn't this allowed ? :
+  // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref);
+}
+
+struct chwirut2_functor : DenseFunctor<double>
+{
+    chwirut2_functor(void) : DenseFunctor<double>(3,54) {}
+    static const double m_x[54];
+    static const double m_y[54];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        int i;
+
+        assert(b.size()==3);
+        assert(fvec.size()==54);
+        for(i=0; i<54; i++) {
+            double x = m_x[i];
+            fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i];
+        }
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==3);
+        assert(fjac.rows()==54);
+        assert(fjac.cols()==3);
+        for(int i=0; i<54; i++) {
+            double x = m_x[i];
+            double factor = 1./(b[1]+b[2]*x);
+            double e = exp(-b[0]*x);
+            fjac(i,0) = -x*e*factor;
+            fjac(i,1) = -e*factor*factor;
+            fjac(i,2) = -x*e*factor*factor;
+        }
+        return 0;
+    }
+};
+const double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0};
+const double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0  };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
+void testNistChwirut2(void)
+{
+  const int n=3;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 0.1, 0.01, 0.02;
+  // do the computation
+  chwirut2_functor functor;
+  LevenbergMarquardt<chwirut2_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 10);
+  VERIFY_IS_EQUAL(lm.njev(), 8);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
+  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
+  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
+
+  /*
+   * Second try
+   */
+  x<< 0.15, 0.008, 0.010;
+  // do the computation
+  lm.resetParameters();
+  lm.setFtol(1.E6*NumTraits<double>::epsilon());
+  lm.setXtol(1.E6*NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 7);
+  VERIFY_IS_EQUAL(lm.njev(), 6);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.6657666537E-01);
+  VERIFY_IS_APPROX(x[1], 5.1653291286E-03);
+  VERIFY_IS_APPROX(x[2], 1.2150007096E-02);
+}
+
+
+struct misra1a_functor : DenseFunctor<double>
+{
+    misra1a_functor(void) : DenseFunctor<double>(2,14) {}
+    static const double m_x[14];
+    static const double m_y[14];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==2);
+        assert(fvec.size()==14);
+        for(int i=0; i<14; i++) {
+            fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ;
+        }
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==2);
+        assert(fjac.rows()==14);
+        assert(fjac.cols()==2);
+        for(int i=0; i<14; i++) {
+            fjac(i,0) = (1.-exp(-b[1]*m_x[i]));
+            fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i]));
+        }
+        return 0;
+    }
+};
+const double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
+const double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
+
+// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
+void testNistMisra1a(void)
+{
+  const int n=2;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 500., 0.0001;
+  // do the computation
+  misra1a_functor functor;
+  LevenbergMarquardt<misra1a_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 19);
+  VERIFY_IS_EQUAL(lm.njev(), 15);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
+  // check x
+  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
+  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
+
+  /*
+   * Second try
+   */
+  x<< 250., 0.0005;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 5);
+  VERIFY_IS_EQUAL(lm.njev(), 4);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01);
+  // check x
+  VERIFY_IS_APPROX(x[0], 2.3894212918E+02);
+  VERIFY_IS_APPROX(x[1], 5.5015643181E-04);
+}
+
+struct hahn1_functor : DenseFunctor<double>
+{
+    hahn1_functor(void) : DenseFunctor<double>(7,236) {}
+    static const double m_x[236];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0  , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0  , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0  , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 ,
+        16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0  , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0   , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 
+12.596E0 , 
+13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0  , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0  , 20.935E0 , 21.035E0 , 20.93E0  , 21.074E0 , 21.085E0 , 20.935E0 };
+
+        //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
+
+        assert(b.size()==7);
+        assert(fvec.size()==236);
+        for(int i=0; i<236; i++) {
+            double x=m_x[i], xx=x*x, xxx=xx*x;
+            fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - m_y[i];
+        }
+        return 0;
+    }
+
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==7);
+        assert(fjac.rows()==236);
+        assert(fjac.cols()==7);
+        for(int i=0; i<236; i++) {
+            double x=m_x[i], xx=x*x, xxx=xx*x;
+            double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
+            fjac(i,0) = 1.*fact;
+            fjac(i,1) = x*fact;
+            fjac(i,2) = xx*fact;
+            fjac(i,3) = xxx*fact;
+            fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
+            fjac(i,4) = x*fact;
+            fjac(i,5) = xx*fact;
+            fjac(i,6) = xxx*fact;
+        }
+        return 0;
+    }
+};
+const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0  , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 ,
+282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 ,
+141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0  , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0  , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0};
+
+// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml
+void testNistHahn1(void)
+{
+  const int  n=7;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 10., -1., .05, -.00001, -.05, .001, -.000001;
+  // do the computation
+  hahn1_functor functor;
+  LevenbergMarquardt<hahn1_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 11);
+  VERIFY_IS_EQUAL(lm.njev(), 10);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.0776351733E+00);
+  VERIFY_IS_APPROX(x[1],-1.2269296921E-01);
+  VERIFY_IS_APPROX(x[2], 4.0863750610E-03);
+  VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06
+  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
+  VERIFY_IS_APPROX(x[5], 2.4053735503E-04);
+  VERIFY_IS_APPROX(x[6],-1.2314450199E-07);
+
+  /*
+   * Second try
+   */
+  x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 11);
+  VERIFY_IS_EQUAL(lm.njev(), 10);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.077640); // should be :  1.0776351733E+00
+  VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01
+  VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03
+  VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06
+  VERIFY_IS_APPROX(x[4],-5.7609940901E-03);
+  VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04
+  VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07
+
+}
+
+struct misra1d_functor : DenseFunctor<double>
+{
+    misra1d_functor(void) : DenseFunctor<double>(2,14) {}
+    static const double x[14];
+    static const double y[14];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==2);
+        assert(fvec.size()==14);
+        for(int i=0; i<14; i++) {
+            fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i];
+        }
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==2);
+        assert(fjac.rows()==14);
+        assert(fjac.cols()==2);
+        for(int i=0; i<14; i++) {
+            double den = 1.+b[1]*x[i];
+            fjac(i,0) = b[1]*x[i] / den;
+            fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den;
+        }
+        return 0;
+    }
+};
+const double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0};
+const double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0};
+
+// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml
+void testNistMisra1d(void)
+{
+  const int n=2;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 500., 0.0001;
+  // do the computation
+  misra1d_functor functor;
+  LevenbergMarquardt<misra1d_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 9);
+  VERIFY_IS_EQUAL(lm.njev(), 7);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
+  // check x
+  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
+  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
+
+  /*
+   * Second try
+   */
+  x<< 450., 0.0003;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 4);
+  VERIFY_IS_EQUAL(lm.njev(), 3);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02);
+  // check x
+  VERIFY_IS_APPROX(x[0], 4.3736970754E+02);
+  VERIFY_IS_APPROX(x[1], 3.0227324449E-04);
+}
+
+
+struct lanczos1_functor : DenseFunctor<double>
+{
+    lanczos1_functor(void) : DenseFunctor<double>(6,24) {}
+    static const double x[24];
+    static const double y[24];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==6);
+        assert(fvec.size()==24);
+        for(int i=0; i<24; i++)
+            fvec[i] = b[0]*exp(-b[1]*x[i]) + b[2]*exp(-b[3]*x[i]) + b[4]*exp(-b[5]*x[i])  - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==6);
+        assert(fjac.rows()==24);
+        assert(fjac.cols()==6);
+        for(int i=0; i<24; i++) {
+            fjac(i,0) = exp(-b[1]*x[i]);
+            fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]);
+            fjac(i,2) = exp(-b[3]*x[i]);
+            fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]);
+            fjac(i,4) = exp(-b[5]*x[i]);
+            fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]);
+        }
+        return 0;
+    }
+};
+const double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 };
+const double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml
+void testNistLanczos1(void)
+{
+  const int n=6;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6;
+  // do the computation
+  lanczos1_functor functor;
+  LevenbergMarquardt<lanczos1_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 2);
+  VERIFY_IS_EQUAL(lm.nfev(), 79);
+  VERIFY_IS_EQUAL(lm.njev(), 72);
+  // check norm^2
+//   VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.430899764097e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
+  // check x
+  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
+  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
+  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
+  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
+  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
+  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
+
+  /*
+   * Second try
+   */
+  x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 2);
+  VERIFY_IS_EQUAL(lm.nfev(), 9);
+  VERIFY_IS_EQUAL(lm.njev(), 8);
+  // check norm^2
+//   VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.428595533845e-25);  // should be 1.4307867721E-25, but nist results are on 128-bit floats
+  // check x
+  VERIFY_IS_APPROX(x[0], 9.5100000027E-02);
+  VERIFY_IS_APPROX(x[1], 1.0000000001E+00);
+  VERIFY_IS_APPROX(x[2], 8.6070000013E-01);
+  VERIFY_IS_APPROX(x[3], 3.0000000002E+00);
+  VERIFY_IS_APPROX(x[4], 1.5575999998E+00);
+  VERIFY_IS_APPROX(x[5], 5.0000000001E+00);
+
+}
+
+struct rat42_functor : DenseFunctor<double>
+{
+    rat42_functor(void) : DenseFunctor<double>(3,9) {}
+    static const double x[9];
+    static const double y[9];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==3);
+        assert(fvec.size()==9);
+        for(int i=0; i<9; i++) {
+            fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i];
+        }
+        return 0;
+    }
+
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==3);
+        assert(fjac.rows()==9);
+        assert(fjac.cols()==3);
+        for(int i=0; i<9; i++) {
+            double e = exp(b[1]-b[2]*x[i]);
+            fjac(i,0) = 1./(1.+e);
+            fjac(i,1) = -b[0]*e/(1.+e)/(1.+e);
+            fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e);
+        }
+        return 0;
+    }
+};
+const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 };
+const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml
+void testNistRat42(void)
+{
+  const int n=3;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 100., 1., 0.1;
+  // do the computation
+  rat42_functor functor;
+  LevenbergMarquardt<rat42_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 10);
+  VERIFY_IS_EQUAL(lm.njev(), 8);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
+  // check x
+  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
+  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
+  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
+
+  /*
+   * Second try
+   */
+  x<< 75., 2.5, 0.07;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 6);
+  VERIFY_IS_EQUAL(lm.njev(), 5);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00);
+  // check x
+  VERIFY_IS_APPROX(x[0], 7.2462237576E+01);
+  VERIFY_IS_APPROX(x[1], 2.6180768402E+00);
+  VERIFY_IS_APPROX(x[2], 6.7359200066E-02);
+}
+
+struct MGH10_functor : DenseFunctor<double>
+{
+    MGH10_functor(void) : DenseFunctor<double>(3,16) {}
+    static const double x[16];
+    static const double y[16];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==3);
+        assert(fvec.size()==16);
+        for(int i=0; i<16; i++)
+            fvec[i] =  b[0] * exp(b[1]/(x[i]+b[2])) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==3);
+        assert(fjac.rows()==16);
+        assert(fjac.cols()==3);
+        for(int i=0; i<16; i++) {
+            double factor = 1./(x[i]+b[2]);
+            double e = exp(b[1]*factor);
+            fjac(i,0) = e;
+            fjac(i,1) = b[0]*factor*e;
+            fjac(i,2) = -b[1]*b[0]*factor*factor*e;
+        }
+        return 0;
+    }
+};
+const double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 };
+const double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml
+void testNistMGH10(void)
+{
+  const int n=3;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 2., 400000., 25000.;
+  // do the computation
+  MGH10_functor functor;
+  LevenbergMarquardt<MGH10_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1); 
+  VERIFY_IS_EQUAL(lm.nfev(), 284 ); 
+  VERIFY_IS_EQUAL(lm.njev(), 249 ); 
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
+  // check x
+  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
+  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
+  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
+
+  /*
+   * Second try
+   */
+  x<< 0.02, 4000., 250.;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 126);
+  VERIFY_IS_EQUAL(lm.njev(), 116);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01);
+  // check x
+  VERIFY_IS_APPROX(x[0], 5.6096364710E-03);
+  VERIFY_IS_APPROX(x[1], 6.1813463463E+03);
+  VERIFY_IS_APPROX(x[2], 3.4522363462E+02);
+}
+
+
+struct BoxBOD_functor : DenseFunctor<double>
+{
+    BoxBOD_functor(void) : DenseFunctor<double>(2,6) {}
+    static const double x[6];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        static const double y[6] = { 109., 149., 149., 191., 213., 224. };
+        assert(b.size()==2);
+        assert(fvec.size()==6);
+        for(int i=0; i<6; i++)
+            fvec[i] =  b[0]*(1.-exp(-b[1]*x[i])) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==2);
+        assert(fjac.rows()==6);
+        assert(fjac.cols()==2);
+        for(int i=0; i<6; i++) {
+            double e = exp(-b[1]*x[i]);
+            fjac(i,0) = 1.-e;
+            fjac(i,1) = b[0]*x[i]*e;
+        }
+        return 0;
+    }
+};
+const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml
+void testNistBoxBOD(void)
+{
+  const int n=2;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 1., 1.;
+  // do the computation
+  BoxBOD_functor functor;
+  LevenbergMarquardt<BoxBOD_functor> lm(functor);
+  lm.setFtol(1.E6*NumTraits<double>::epsilon());
+  lm.setXtol(1.E6*NumTraits<double>::epsilon());
+  lm.setFactor(10);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 31);
+  VERIFY_IS_EQUAL(lm.njev(), 25);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
+  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
+
+  /*
+   * Second try
+   */
+  x<< 100., 0.75;
+  // do the computation
+  lm.resetParameters();
+  lm.setFtol(NumTraits<double>::epsilon());
+  lm.setXtol( NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1); 
+  VERIFY_IS_EQUAL(lm.nfev(), 15 ); 
+  VERIFY_IS_EQUAL(lm.njev(), 14 ); 
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 2.1380940889E+02);
+  VERIFY_IS_APPROX(x[1], 5.4723748542E-01);
+}
+
+struct MGH17_functor : DenseFunctor<double>
+{
+    MGH17_functor(void) : DenseFunctor<double>(5,33) {}
+    static const double x[33];
+    static const double y[33];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==5);
+        assert(fvec.size()==33);
+        for(int i=0; i<33; i++)
+            fvec[i] =  b[0] + b[1]*exp(-b[3]*x[i]) +  b[2]*exp(-b[4]*x[i]) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==5);
+        assert(fjac.rows()==33);
+        assert(fjac.cols()==5);
+        for(int i=0; i<33; i++) {
+            fjac(i,0) = 1.;
+            fjac(i,1) = exp(-b[3]*x[i]);
+            fjac(i,2) = exp(-b[4]*x[i]);
+            fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]);
+            fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]);
+        }
+        return 0;
+    }
+};
+const double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 };
+const double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml
+void testNistMGH17(void)
+{
+  const int n=5;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 50., 150., -100., 1., 2.;
+  // do the computation
+  MGH17_functor functor;
+  LevenbergMarquardt<MGH17_functor> lm(functor);
+  lm.setFtol(NumTraits<double>::epsilon());
+  lm.setXtol(NumTraits<double>::epsilon());
+  lm.setMaxfev(1000);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 2); 
+//   VERIFY_IS_EQUAL(lm.nfev(), 602 ); 
+  VERIFY_IS_EQUAL(lm.njev(), 545 ); 
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
+  // check x
+  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
+  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
+  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
+  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
+  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
+
+  /*
+   * Second try
+   */
+  x<< 0.5  ,1.5  ,-1   ,0.01 ,0.02;
+  // do the computation
+  lm.resetParameters();
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 18);
+  VERIFY_IS_EQUAL(lm.njev(), 15);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05);
+  // check x
+  VERIFY_IS_APPROX(x[0], 3.7541005211E-01);
+  VERIFY_IS_APPROX(x[1], 1.9358469127E+00);
+  VERIFY_IS_APPROX(x[2], -1.4646871366E+00);
+  VERIFY_IS_APPROX(x[3], 1.2867534640E-02);
+  VERIFY_IS_APPROX(x[4], 2.2122699662E-02);
+}
+
+struct MGH09_functor : DenseFunctor<double>
+{
+    MGH09_functor(void) : DenseFunctor<double>(4,11) {}
+    static const double _x[11];
+    static const double y[11];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==4);
+        assert(fvec.size()==11);
+        for(int i=0; i<11; i++) {
+            double x = _x[i], xx=x*x;
+            fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i];
+        }
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==4);
+        assert(fjac.rows()==11);
+        assert(fjac.cols()==4);
+        for(int i=0; i<11; i++) {
+            double x = _x[i], xx=x*x;
+            double factor = 1./(xx+x*b[2]+b[3]);
+            fjac(i,0) = (xx+x*b[1]) * factor;
+            fjac(i,1) = b[0]*x* factor;
+            fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor;
+            fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor;
+        }
+        return 0;
+    }
+};
+const double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01,  1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 };
+const double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml
+void testNistMGH09(void)
+{
+  const int n=4;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 25., 39, 41.5, 39.;
+  // do the computation
+  MGH09_functor functor;
+  LevenbergMarquardt<MGH09_functor> lm(functor);
+  lm.setMaxfev(1000);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1); 
+  VERIFY_IS_EQUAL(lm.nfev(), 490 ); 
+  VERIFY_IS_EQUAL(lm.njev(), 376 ); 
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
+  // check x
+  VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01
+  VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01
+  VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01
+  VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01
+
+  /*
+   * Second try
+   */
+  x<< 0.25, 0.39, 0.415, 0.39;
+  // do the computation
+  lm.resetParameters();
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 18);
+  VERIFY_IS_EQUAL(lm.njev(), 16);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04);
+  // check x
+  VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01
+  VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01
+  VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01
+  VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01
+}
+
+
+
+struct Bennett5_functor : DenseFunctor<double>
+{
+    Bennett5_functor(void) : DenseFunctor<double>(3,154) {}
+    static const double x[154];
+    static const double y[154];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==3);
+        assert(fvec.size()==154);
+        for(int i=0; i<154; i++)
+            fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==3);
+        assert(fjac.rows()==154);
+        assert(fjac.cols()==3);
+        for(int i=0; i<154; i++) {
+            double e = pow(b[1]+x[i],-1./b[2]);
+            fjac(i,0) = e;
+            fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]);
+            fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2];
+        }
+        return 0;
+    }
+};
+const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0,
+11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 };
+const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0
+,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,
+-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml
+void testNistBennett5(void)
+{
+  const int  n=3;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< -2000., 50., 0.8;
+  // do the computation
+  Bennett5_functor functor;
+  LevenbergMarquardt<Bennett5_functor> lm(functor);
+  lm.setMaxfev(1000);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 758);
+  VERIFY_IS_EQUAL(lm.njev(), 744);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
+  // check x
+  VERIFY_IS_APPROX(x[0], -2.5235058043E+03);
+  VERIFY_IS_APPROX(x[1], 4.6736564644E+01);
+  VERIFY_IS_APPROX(x[2], 9.3218483193E-01);
+  /*
+   * Second try
+   */
+  x<< -1500., 45., 0.85;
+  // do the computation
+  lm.resetParameters();
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 203);
+  VERIFY_IS_EQUAL(lm.njev(), 192);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04);
+  // check x
+  VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03
+  VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01);
+  VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01);
+}
+
+struct thurber_functor : DenseFunctor<double>
+{
+    thurber_functor(void) : DenseFunctor<double>(7,37) {}
+    static const double _x[37];
+    static const double _y[37];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        //        int called=0; printf("call hahn1_functor with  iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++;
+        assert(b.size()==7);
+        assert(fvec.size()==37);
+        for(int i=0; i<37; i++) {
+            double x=_x[i], xx=x*x, xxx=xx*x;
+            fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i];
+        }
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==7);
+        assert(fjac.rows()==37);
+        assert(fjac.cols()==7);
+        for(int i=0; i<37; i++) {
+            double x=_x[i], xx=x*x, xxx=xx*x;
+            double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx);
+            fjac(i,0) = 1.*fact;
+            fjac(i,1) = x*fact;
+            fjac(i,2) = xx*fact;
+            fjac(i,3) = xxx*fact;
+            fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact;
+            fjac(i,4) = x*fact;
+            fjac(i,5) = xx*fact;
+            fjac(i,6) = xxx*fact;
+        }
+        return 0;
+    }
+};
+const double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 };
+const double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0};
+
+// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml
+void testNistThurber(void)
+{
+  const int n=7;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ;
+  // do the computation
+  thurber_functor functor;
+  LevenbergMarquardt<thurber_functor> lm(functor);
+  lm.setFtol(1.E4*NumTraits<double>::epsilon());
+  lm.setXtol(1.E4*NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 39);
+  VERIFY_IS_EQUAL(lm.njev(), 36);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
+  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
+  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
+  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
+  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
+  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
+  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
+
+  /*
+   * Second try
+   */
+  x<< 1300 ,1500 ,500  ,75   ,1    ,0.4  ,0.05  ;
+  // do the computation
+  lm.resetParameters();
+  lm.setFtol(1.E4*NumTraits<double>::epsilon());
+  lm.setXtol(1.E4*NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 29);
+  VERIFY_IS_EQUAL(lm.njev(), 28);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.2881396800E+03);
+  VERIFY_IS_APPROX(x[1], 1.4910792535E+03);
+  VERIFY_IS_APPROX(x[2], 5.8323836877E+02);
+  VERIFY_IS_APPROX(x[3], 7.5416644291E+01);
+  VERIFY_IS_APPROX(x[4], 9.6629502864E-01);
+  VERIFY_IS_APPROX(x[5], 3.9797285797E-01);
+  VERIFY_IS_APPROX(x[6], 4.9727297349E-02);
+}
+
+struct rat43_functor : DenseFunctor<double>
+{
+    rat43_functor(void) : DenseFunctor<double>(4,15) {}
+    static const double x[15];
+    static const double y[15];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==4);
+        assert(fvec.size()==15);
+        for(int i=0; i<15; i++)
+            fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==4);
+        assert(fjac.rows()==15);
+        assert(fjac.cols()==4);
+        for(int i=0; i<15; i++) {
+            double e = exp(b[1]-b[2]*x[i]);
+            double power = -1./b[3];
+            fjac(i,0) = pow(1.+e, power);
+            fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.);
+            fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.);
+            fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power);
+        }
+        return 0;
+    }
+};
+const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. };
+const double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml
+void testNistRat43(void)
+{
+  const int n=4;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 100., 10., 1., 1.;
+  // do the computation
+  rat43_functor functor;
+  LevenbergMarquardt<rat43_functor> lm(functor);
+  lm.setFtol(1.E6*NumTraits<double>::epsilon());
+  lm.setXtol(1.E6*NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 27);
+  VERIFY_IS_EQUAL(lm.njev(), 20);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
+  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
+  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
+  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
+
+  /*
+   * Second try
+   */
+  x<< 700., 5., 0.75, 1.3;
+  // do the computation
+  lm.resetParameters();
+  lm.setFtol(1.E5*NumTraits<double>::epsilon());
+  lm.setXtol(1.E5*NumTraits<double>::epsilon());
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 9);
+  VERIFY_IS_EQUAL(lm.njev(), 8);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 6.9964151270E+02);
+  VERIFY_IS_APPROX(x[1], 5.2771253025E+00);
+  VERIFY_IS_APPROX(x[2], 7.5962938329E-01);
+  VERIFY_IS_APPROX(x[3], 1.2792483859E+00);
+}
+
+
+
+struct eckerle4_functor : DenseFunctor<double>
+{
+    eckerle4_functor(void) : DenseFunctor<double>(3,35) {}
+    static const double x[35];
+    static const double y[35];
+    int operator()(const VectorXd &b, VectorXd &fvec)
+    {
+        assert(b.size()==3);
+        assert(fvec.size()==35);
+        for(int i=0; i<35; i++)
+            fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i];
+        return 0;
+    }
+    int df(const VectorXd &b, MatrixXd &fjac)
+    {
+        assert(b.size()==3);
+        assert(fjac.rows()==35);
+        assert(fjac.cols()==3);
+        for(int i=0; i<35; i++) {
+            double b12 = b[1]*b[1];
+            double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12);
+            fjac(i,0) = e / b[1];
+            fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12;
+            fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12;
+        }
+        return 0;
+    }
+};
+const double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0};
+const double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 };
+
+// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml
+void testNistEckerle4(void)
+{
+  const int n=3;
+  int info;
+
+  VectorXd x(n);
+
+  /*
+   * First try
+   */
+  x<< 1., 10., 500.;
+  // do the computation
+  eckerle4_functor functor;
+  LevenbergMarquardt<eckerle4_functor> lm(functor);
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 18);
+  VERIFY_IS_EQUAL(lm.njev(), 15);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.5543827178);
+  VERIFY_IS_APPROX(x[1], 4.0888321754);
+  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
+
+  /*
+   * Second try
+   */
+  x<< 1.5, 5., 450.;
+  // do the computation
+  info = lm.minimize(x);
+
+  // check return value
+  VERIFY_IS_EQUAL(info, 1);
+  VERIFY_IS_EQUAL(lm.nfev(), 7);
+  VERIFY_IS_EQUAL(lm.njev(), 6);
+  // check norm^2
+  VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03);
+  // check x
+  VERIFY_IS_APPROX(x[0], 1.5543827178);
+  VERIFY_IS_APPROX(x[1], 4.0888321754);
+  VERIFY_IS_APPROX(x[2], 4.5154121844E+02);
+}
+
+void test_levenberg_marquardt()
+{
+    // Tests using the examples provided by (c)minpack
+    CALL_SUBTEST(testLmder1());
+    CALL_SUBTEST(testLmder());
+    CALL_SUBTEST(testLmdif1());
+//     CALL_SUBTEST(testLmstr1());
+//     CALL_SUBTEST(testLmstr());
+    CALL_SUBTEST(testLmdif());
+
+    // NIST tests, level of difficulty = "Lower"
+    CALL_SUBTEST(testNistMisra1a());
+    CALL_SUBTEST(testNistChwirut2());
+
+    // NIST tests, level of difficulty = "Average"
+    CALL_SUBTEST(testNistHahn1());
+    CALL_SUBTEST(testNistMisra1d());
+    CALL_SUBTEST(testNistMGH17());
+    CALL_SUBTEST(testNistLanczos1());
+
+//     // NIST tests, level of difficulty = "Higher"
+    CALL_SUBTEST(testNistRat42());
+    CALL_SUBTEST(testNistMGH10());
+    CALL_SUBTEST(testNistBoxBOD());
+//     CALL_SUBTEST(testNistMGH09());
+    CALL_SUBTEST(testNistBennett5());
+    CALL_SUBTEST(testNistThurber());
+    CALL_SUBTEST(testNistRat43());
+    CALL_SUBTEST(testNistEckerle4());
+}